{"paper_id":"e74fc1ce-e0d6-425d-96cd-5492b039533c","body_text":"The B-cell lymphoma 2\n(BCL-2) family plays a central role in the\nregulation of apoptosis, and dysregulation of this pathway is a hallmark\nof many cancers. \n − \n \n \n \n \n  Overexpression of antiapoptotic proteins such as BCL-2 and BCL-xL\nenables malignant cells to evade programmed cell death and contributes\nto therapeutic resistance. \n , − \n \n \n  As a result, BCL-2 has become a clinically validated target for\nanticancer drug discovery, most notably through the development of\nBH3-mimetic inhibitors such as Venetoclax. \n , , − \n \n  Although Venetoclax has\ndemonstrated major therapeutic benefit, particularly in hematologic\nmalignancies, the continued need for compounds with improved selectivity,\nbroader applicability, and favorable safety profiles has sustained\ninterest in the discovery of novel BCL-2-targeting ligands. \n , , −\nDrug repurposing offers an attractive strategy\nfor accelerating\nthis process because approved drugs already possess substantial pharmacological\nand safety information. \n , \n  In the context of BCL-2,\nrepurposing may enable rapid identification of chemically diverse\nscaffolds capable of modulating apoptosis-related signaling. However,\nsuccessful repurposing remains challenging because accurate prioritization\ndepends on more than ligand similarity or static binding scores. In\nparticular, the binding pose of the ligand within the BCL-2 binding\ngroove, the induced structural adaptation of the target protein, and\nthe resulting residue-level interaction network can all influence\nwhether a candidate behaves as a productive binder. Thus, computational\nworkflows that can capture both plausible bound structures and their\ndynamic consequences are especially valuable for repurposing studies\ntargeting BCL-2.\nTraditional target-based screening methods,\nincluding molecular\ndocking and molecular dynamics (MD) simulations, remain widely used\nfor protein–ligand modeling, but they carry important limitations.\nDocking methods rely on simplified sampling schemes and approximate\nscoring functions, which can inadequately represent structural cooperativity\nbetween ligand binding and receptor rearrangement. Although MD simulations\nprovide full time-dependent trajectory information, their practical\napplication in large-scale screening is constrained by substantial\ncomputational cost and the need for extensive sampling to approach\nconvergence. Consequently, standard MD-based workflows are most informative\nwhen applied to a focused set of preselected candidates rather than\nto large compound libraries. By contrast, recent deep generative models\nprovide an opportunity to predict protein–ligand complexes\nin a more state-aware manner. Such as NeuralPlexer is a multiscale\ngenerative framework that predicts protein–ligand complex structures\ndirectly from protein sequence and ligand graph representations. \n − \n \n  By combining residue-level contact prediction with progressive all-atom\nrefinement and diffusion-based ligand coordinate generation, NeuralPlexer\ncan model bound complexes while incorporating key biophysical constraints\nsuch as bond geometry and steric compatibility. \n − \n \n  These features\nmake it particularly attractive for large-scale generation of ligand-specific\ncomplex conformations prior to downstream energetic analysis. NeuralPlexer\ncaptures the structural cooperativity between ligand binding and protein\nconformation, overcoming a major limitation of conventional docking\nmethods that often treat the receptor as rigid. Using a diffusion-based\ngenerative process, it iteratively samples diverse ligand poses while\nenforcing key biophysical constraints, including bond lengths, bond\nangles, and steric compatibility. By incorporating these priors directly\ninto model generation, NeuralPlexer produces binding-site conformations\nthat are both geometrically plausible and energetically favorable.\nIts performance has been extensively benchmarked in tasks such as\nflexible binding-site structure recovery and blind protein–ligand\nbinding. \n − \n \n  In these evaluations, NeuralPlexer has consistently\noutperformed state-of-the-art approaches, including AlphaFold2, in\nboth global structure accuracy and the prediction of ligand-induced\nconformational changes. \n −\nEven when plausible bound\ncomplexes have been generated and MD\nsimulations have been performed, prioritization reduced to a single\nend point conformation, such as the lowest-energy docked pose or a\ntrajectory end-state, may fail to capture the time-dependent fluctuations\nand ligand-induced conformational changes that govern binding stability.\nApproaches that extract and integrate interaction patterns from MD\ntrajectories across many frames, therefore, provide a level of mechanistic\ndiscrimination that single-structure or single-point energetic scoring\nalone may not achieve. Neural relational inference (NRI) is a graph-based\nvariational framework that learns latent interaction maps from dynamical\nsystems and has recently been applied to protein motions. \n , \n  In its original form, NRI was used primarily to study residue–residue\ncommunication in proteins. Here, first time in literature, NRI was\nincorporated to go beyond per-residue contact frequency analysis.\nBy treating each residue as a node in a latent interaction graph and\nlearning edge types from trajectory data, NRI can detect noncovalent\ninteraction channels, including allosteric communication pathways,\nthat remain invisible to conventional distance-based contact maps.\nThis is particularly relevant for BCL-2, where induced-fit rearrangements\nof the hydrophobic groove upon ligand binding are mechanistically\nimportant and not fully captured by pairwise distance cutoffs alone.\nSuch an approach could help distinguish compounds that merely occupy\nthe binding pocket from those that reproduce a more functionally relevant\ndynamic interaction architecture.\nThus, in this study, we developed\na multimodal repurposing workflow\ncentered on diffusion-based generative model complex generation and\nan extended NRI analysis of protein–ligand MD trajectories\nto identify FDA-approved compounds with potential BCL-2 inhibitory\nactivity. A library of 3094 approved drugs was screened by NeuralPlexer\nto generate ligand-specific BCL-2 complex conformations. This process\nyielded around 1300 successful binding poses, representing individual\n3D conformations of BCL-2 for each bound ligand. Subsequently, we\ncalculated the interaction energies (i.e., score-in-place) for these\nprotein–ligand poses, which were then filtered and prioritized\nusing docking refinement, anticancer Quantitative Structure–Activity\nRelationships (QSAR) classification, MD simulations, and Molecular\nMechanics/Generalized Born Surface Area (MM/GBSA) binding free-energy\ncalculations. In parallel, a Chemprop-based pIC 50  prediction\nmodel  and GNINA docking  were used as an orthogonal ligand-centric prioritization\narm. We then applied the extended NRI framework to characterize residue-ligand\nand residue–residue dynamic couplings and to compare candidate\ninteraction signatures against Venetoclax. This integrated strategy\nprioritized the most compelling hit compounds, and subsequent Time-Resolved\nFluorescence Resonance Energy Transfer (TR-FRET) and cell-viability\nassays provided initial validation of selected hit compounds as potential\nBCL-2 inhibitors. To the best of our knowledge, this work represents\nthe first application of diffusion-based generative modeling to large-scale\nscreening of an FDA-approved drug library against BCL-2. Collectively,\nour findings illustrate how generative complex prediction and graph-based\ndynamic inference can be combined to support structure-guided drug\nrepurposing for BCL-2-driven malignancies.\n\nIn this study, we performed\nan  in silico  analysis to identify FDA-approved compounds\nthat show promise in targeting the BCL-2 protein, an important apoptosis\nregulator and a desirable target for cancer treatment. Utilizing the\ncutting-edge protein structure prediction model NeuralPlexer, we produced\na library of 3094 FDA-approved compounds as well as precise and sensitive\nconformations of the BCL-2 protein. NeuralPlexer derived 8 conformations\nfor every protein–ligand interaction using the whole sequence\nof the BCL-2 protein and then ranked them according to confidence\nscores. Therefore, we were able to focus on the most promising binding\npositions for additional examinations. The protein–ligand complexes\nobtained for each BCL-2/ligand complex were re-scored by Glide SP\n(Standard Precision) with the “refine only” option,\nwhich optimizes the ligand conformation without requiring any sampling\nin docking. This method enabled efficient evaluation of the binding\nof FDA-approved drugs within the BCL-2 binding site. The initial screening\nlibrary contained 3094 FDA-approved compounds obtained from DrugBank.\nFor each compound, NeuralPlexer generated up to eight BCL-2–ligand\ncomplex conformations, and the highest-confidence poses were evaluated\nfor structural acceptability based on pose geometry and model confidence.\nThis first filtering step retained 1294 structurally acceptable successful\nBCL-2 ligand complexes for downstream rescoring. These complexes were\nthen subjected to Glide SP scoring in “refine only”\nmode, which locally optimized the NeuralPlexer-generated ligand pose\nwithout performing full redocking. To ensure that Glide refinement\ndid not substantially alter the original NeuralPlexer binding mode,\nonly complexes with a ligand RMSD <2 Å between the NeuralPlexer-generated\nand Glide-refined poses were retained, yielding 463 compounds. These\n463 compounds were ranked according to their Glide SP docking scores,\nand the 20 compounds with the most favorable scores were selected\nfor further filtering. Finally, the top-20 compounds were evaluated\nusing the MetaCore anticancer binary QSAR model, and compounds with\nanticancer activity scores >0.5 were advanced, resulting in 10\nfinal\ncandidates for all-atom MD simulations, MM/GBSA binding free-energy\ncalculations, and NRI-based trajectory analysis. Thus, the final compound\nset was obtained through a sequential, stage-gated prioritization\nstrategy rather than by a single docking or scoring criterion.\nImportantly, compound prioritization in this workflow was implemented\nas a sequential, stage-gated decision process rather than as a single-score\nranking procedure. NeuralPlexer was first used to generate ligand-specific\nBCL-2 complex conformations, and only structurally acceptable poses\nwere retained for downstream assessment. Glide SP “refine only”\nscoring was then applied to evaluate whether the NeuralPlexer-generated\nbinding geometries could be locally optimized while preserving the\npredicted pose, with RMSD filtering used to exclude complexes that\nunderwent substantial pose rearrangement during refinement. The retained\ncompounds were subsequently filtered by predicted anticancer activity,\nfollowed by MD simulations and MM/GBSA calculations to evaluate dynamic\nstability and relative binding energetics. Thus, molecular docking\nand MM/GBSA were used as prioritization and enrichment tools, whereas\nthe final candidate selection was based on convergence across structural\nplausibility, pose stability, predicted anticancer relevance, MD/MM/GBSA\nbehavior, NRI-derived dynamic interaction signatures, and experimental\nvalidation.\nFigure  \n  shows a\nviolin plot of the docking score distribution for all 1294 compounds,\ntogether with the fractions of molecules exhibiting RMSD values above\nand below 2 Å. The 1294 docked compounds were then ranked according\nto their Glide scores, and the top-20 scored candidates were selected\nfor further analysis. These compounds were subsequently evaluated\nusing a binary QSAR model in MetaCore ( https://portal.genego.com/ ) to predict anticancer activity.  Among\nthe top-20 compounds, 10 of them with predicted anticancer activity\nscores greater than 0.5 and favorable Glide scores were retained for\nsubsequent studies. These selected candidates were subjected to 3\nindependent replicates of 200 ns all-atom MD simulations in Desmond,\nand their average binding free energies were estimated using the MM/GBSA\nmethod. For comparison, the same MD and MM/GBSA analyses were also\ncarried out for Venetoclax, the FDA-approved BCL-2 inhibitor used\nas a reference compound.\n(A) The violin plot of the distribution of all\n1294 molecules with\nthe score-in-place docking score and the percentage of molecules with\nRMSD over and under 2 Å; (B) 3D structure of the BCL-2 protein\nin complex with Venetoclax on the left and Relugolix on the right;\n(C) the interaction types of Venetoclax and Relugolix in the left\nand right, respectively, during the 200 ns MD simulations; (D) the\nheat map of the BCL-2 residues interacting with Venetoclax on the\nleft and Relugolix on the right.\nVenetoclax showed the most favorable average MM/GBSA\nbinding free\nenergy (−132.14 kcal/mol), exceeding those of the top candidate\ncompounds identified in this study. However, because Venetoclax is\nrelatively large (molecular weight 868 g/mol), MM/GBSA values were\nadditionally normalized on a per-non-hydrogen-atom basis to reduce\nsize-related bias and enable comparison in terms of ligand efficiency.\nLigand efficiency, which reflects binding energy relative to molecular\nsize, is a useful metric in drug discovery because compounds with\nfavorable efficiency often display more balanced drug-like properties.\nBased on this analysis, only Lathyrol may not show better ligand efficiency\nthan Venetoclax. The corresponding ligand efficiency value of Lathyrol\nis −1.87 kcal/mol, compared with −2.17 kcal/mol for\nVenetoclax. These results suggest that, despite weaker average MM/GBSA\nbinding energies than Venetoclax, these compounds may engage the BCL-2\nbinding site efficiently relative to their molecular size. All scores\nfor the 10 FDA-approved compounds and Venetoclax are summarized in  Table  \n , and their 2D structures\nare shown in  Figure  \n .\nNo docking score for Venetoclax.\n2D structures of selected FDA-approved compounds with Venetoclax\nas the reference molecule.\nThe NRI-based analysis was incorporated to provide\na trajectory-level\nmechanistic interpretation that is not captured by static docking\nscores or end point binding free-energy estimates. While docking and\nMM/GBSA can prioritize compounds according to pose compatibility and\nrelative energetic favorability, they do not directly describe how\nligand binding reshapes residue-level communication within BCL-2 during\nthe simulations. By extending NRI to include both protein residues\nand ligand heavy atoms, we were able to learn dynamic residue–residue\nand residue–ligand coupling patterns from the MD trajectories.\nThis allowed candidate compounds to be evaluated not only by their\npredicted binding strength, but also by the extent to which they reproduced\nthe dynamic interaction architecture observed for the reference inhibitor\nVenetoclax. Therefore, NRI served as an additional mechanistic filter\nto distinguish ligands that merely occupied the BCL-2 binding groove\nfrom those that induced a more reference-like dynamic interaction\nnetwork. The NRI model utilizes a variational graph autoencoder architecture\nto recreate the trajectories of MD simulations by learning latent\ncontact maps between protein residues and ligand heavy atoms. While\nprevious NRI models focused on understanding protein dynamics in apo\nand holo forms by tracking the carbon alpha atoms of amino acids,\nwe have extended the NRI model to track both carbon alpha atoms and\nligand heavy atoms. This enhancement allows us to quantify ligand–protein\ninteractions, protein stability, and interactions between ligand heavy\natoms. By training the NRI model to mimic MD trajectories, we obtain\na latent interaction matrix  z \n \n ij \n  (0 ≤  z \n \n ij \n  ≤ 1) that quantifies the contribution of atom  j  to the dynamics of atom  i . The protein–ligand\ninteraction energy is evaluated using the matrix  z \n \n ij \n  with details provided in the Methodology\nsection. We trained 10 NRI models for each identified candidate compounds\nand Venetoclax, using the same hyperparameters for the variational\ngraph autoencoder models.  Figure  \n A illustrates the interaction matrices  z \n \n ij \n  for all ligands with a better efficiency\nscore than Venetoclax, Relugolix as the lowest MM/GBSA score compound,\nand the reference compound. Etofylline’s heavy atoms significantly\naffect residue dynamics, according to the protein–ligand interaction\nmaps, which results in a higher normalized protein–ligand interaction\nenergy than any of the other compounds that were evaluated. Further\ndemonstrating its robust binding potential, Relugolix, which was chosen\nas the option with the lowest MM/GBSA binding energy, also has a greater\nprotein–ligand interaction energy than Venetoclax (except for\nresidue 81). The contact maps of BCL-2 in complex with Etofylline,\nRelugolix, and Venetoclax are shown in  Figure  \n B–D, respectively. Etofylline forms\na large number of interactions with the BCL-2 protein in spite of\nits comparatively small molecular size. However, the concentrated\nresidue–residue interactions that are primarily seen between\nresidues 70 and 100 suggest that it is unable to appreciably stabilize\nthe overall BCL-2 structure. Relugolix, on the other hand, has a contact\npattern that is quite similar to Venetoclax’s, indicating a\ncomparable binding mechanism. However, in some areas, Relugolix creates\nfewer intraprotein interactions in BCL-2 than Venetoclax, which would\nindicate a marginally diminished stabilizing impact after binding.\nResults\nof the NRI model are shown as follows: (A) the mean energy\ncontribution of the heavy atoms of all ligands with a better efficiency\nscore than Venetoclax, Relugolix as the lowest MM/GBSA score compound,\nand Venetoclax as reference with BCL-2 protein residues; (B) a contact\nmap of BCL-2 with Etofylline; (C) Relugolix, and (D) Venetoclax.\nAn ML-based prioritization\nbranch was also developed alongside the NeuralPlexer-guided workflow\nto provide an orthogonal framework for candidate selection. NeuralPlexer\noffers a structure-aware route for modeling ligand-specific BCL-2\ncomplexes, but its ranking of compounds is inherently linked to generated\nposes and subsequent physics-based evaluation. By contrast, Chemprop\nprovides ligand-based activity prediction directly from molecular\ngraphs and is therefore independent of explicit complex generation.\nThe integration of these complementary approaches was intended to\nminimize reliance on a single computational paradigm, improve the\nrobustness of hit prioritization, and reveal candidates supported\nby both structure-based and data-driven evidence. In this way, the\nML-guided arm served not only as an additional screening route, but\nalso as an internal cross-validation strategy prior to biochemical\nand cellular testing.\nAs described in the Methods section, we\nalso trained a Chemprop  model to predict\npIC 50  values using a data set of 1067 known BCL-2 inhibitors\nfrom ChEMBL ( https://www.ebi.ac.uk/chembl/ ) with experimentally determined activities. Model training was performed\nusing 5-fold cross-validation, with RMSE selected as the primary optimization\nmetric. Over 50 training epochs, the model showed progressive improvement\nin both RMSE and  R \n 2  across the validation\nsets. The best-performing fold achieved a validation RMSE of 0.636\nand an  R \n 2  of 0.889, indicating strong\nagreement between predicted and experimental pIC 50  values.\nThis optimal performance, reached at epoch 47, suggests that the model\ngeneralizes well across structurally diverse BCL-2 inhibitors. This\nvalidation performance was consistent with the results obtained on\nthe independent test set, where the final model achieved an RMSE of\n0.748 and an  R \n 2  of 0.858, supporting its\nrobustness and predictive reliability. The reproducibility of these\nresults across cross-validation folds further indicates that the Chemprop\narchitecture is stable for this regression task. Together, these metrics\nsuggest that the model is sufficiently accurate for downstream applications\nsuch as virtual screening of FDA-approved compounds and falls within\na reasonable range for ligand–target affinity prediction. Based\non the predicted pIC 50  values, compounds estimated to be\nmore potent than Venetoclax were selected for further evaluation.\nBecause pIC 50  values can be influenced by molecular size,\nligand efficiency was also calculated to enable size-normalized comparison\nof compound potency. Given the relatively high molecular weight of\nVenetoclax, all selected compounds showed better ligand efficiency\nvalues, thereby prioritizing smaller molecules with favorable predicted\npotency and allowing a more balanced comparison.\nGNINA,  which incorporates convolutional\nneural network (CNN)-based scoring, was then used to dock the selected\ncompounds into the BCL-2 binding site. The resulting poses were filtered\nusing CNN affinity scores, and only compounds with values greater\nthan 6.0 were retained, yielding a total of 10 candidates. Notably,\nwhen normalized for molecular size, all 10 compounds exhibited CNN-based\nligand efficiency values higher than that of Venetoclax, suggesting\nfavorable binding potential relative to their size. These top-ranked\ncompounds were subsequently evaluated by 3 independent replicates\nof 200 ns all-atom MD simulations, from which 2000 trajectory frames\nwere collected for each simulation (i.e., 6000 frames per protein–ligand\ncomplex). Binding free energies were then estimated using the MM/GBSA\nmethod implemented in the Schrödinger Prime module. To further\nassess their potential anticancer relevance, the 10 candidates were\nanalyzed using the MetaCore platform, which integrates cheminformatics\nand systems biology for functional prediction. Based on structural\nfeatures and predicted interaction patterns, 6 of the 10 compounds\nshowed anticancer activity scores above 0.5, the default MetaCore\nthreshold for probable anticancer potential. Among these candidates,\nBaricitinib, Trazodone, and Buspirone displayed strong CNN affinity\nscores (≥6.5), high predicted pIC 50  values (≥9.1),\nand anticancer activity scores of at least 0.8, consistent with favorable\nbinding, predicted potency, and potential therapeutic relevance. In\nparallel, Lanraplenib, Trazodone, and Cobimetinib yielded the most\nfavorable MM/GBSA binding free energies (−70.49, −62.29,\nand −58.34 kcal/mol, respectively), further supporting their\npotential as BCL-2 binders.  Table  \n  summarizes the predicted pIC 50  values,\nCNN docking scores, MetaCore activity probabilities, and MM/GBSA free\nenergies for all selected compounds, while their 2D structures are\nshown in  Figure  \n .\n2D structures of FDA-approved compounds selected through the ML-based\nprioritization workflow.\nTo determine whether this ML-guided strategy could\nidentify potent\nBCL-2 candidates complementary to those prioritized by the NeuralPlexer-based\nworkflow, we compared the outputs of the two pipelines. Both the NeuralPlexer-based\nand Chemprop-guided approaches identified some of the compounds with\npotential BCL-2 inhibitory activities. These selected candidates were\nthen prioritized for experimental follow-up and 12 of them subjected\nto in vitro enzyme inhibition and cell viability assays to assess\ntheir biological activity and therapeutic potential.\nA time-resolved fluorescence resonance energy\ntransfer (TR-FRET) assay was performed to assess whether the selected\ncandidate compounds could disrupt the interaction between BCL-2 and\nits binding ligand. This assay provided a biochemical readout of inhibitory\nactivity by measuring the extent to which each compound interfered\nwith ligand engagement. Initial screening of 12 compounds showed that\n4 molecules were able to compete with the known BCL-2 inhibitor, indicating\nmeasurable inhibitory activity. In most of the cases, inhibition of\nBCL-2–ligand binding increased with compound concentration,\nconsistent with a concentration-dependent effect. Among the tested\nmolecules, Vericiguat, Trazodone, Relugolix, and Chlorzoxazone produced\nthe strongest inhibition at 100 μM ( Figure  \n ). In contrast, the remaining 8 compounds\nshowed small or no detectable inhibition and did not display a meaningful\nconcentration-dependent response.\nInhibitory activity of the 12 selected\nhit compounds, with Venetoclax\nincluded as a positive control. Compounds were evaluated at 100 μM\nand 1 μM, and percentage inhibition was calculated using the\nequation described in the Materials and Methods section.\nBased on the initial screening results, four compounds,\nVericiguat,\nLathyrol, Relugolix, and Chlorzoxazone, were selected for more detailed\nevaluation across five concentration levels ( Figure  \n ). Although Trazodone showed notable activity\nin the first-pass assay, Lathyrol was prioritized for dose–response\nanalysis because its initial profile suggested the possibility of\ninhibitory activity at lower concentrations. For this reason, Lathyrol\nwas advanced to broader concentration-range testing in place of Trazodone.\nThis second-stage analysis enabled a more refined comparison of inhibitory\npotency among the selected hits. To examine whether these biochemically\nactive compounds also affected cancer cell growth, they were further\ntested in LN-18 glioma cells using a 24 h MTT cell viability assay.\nTreatment with the selected compounds led to reduced cell proliferation,\nproviding preliminary cellular support for their predicted BCL-2 inhibitory\nactivity ( Figure  \n ).\nTaken together, the TR-FRET and cell viability results strengthened\nthe prioritization of these compounds, particularly those that demonstrated\nboth direct inhibition of BCL-2–ligand binding and measurable\nantiproliferative effects in cancer cells.\nMaximum percentage inhibitory\nactivity of four hit compounds across\nthe tested concentration range (100 μM to 10 nM), together with\nthe corresponding IC 50  values determined after 3 h of incubation.\nMTT-based cell proliferation assay of compounds that showed\ninhibitory\nactivity. Each compound was tested over a concentration range from\n100 μM to 1 nM. Cell viability measurements were performed after\n24 h of treatment in LN-18 glioma cells.\nAmong the four compounds (Vericiguat, Lathyrol,\nRelugolix, and\nChlorzoxazone) evaluated in the TR-FRET assay, Relugolix showed the\nstrongest inhibitory activity, with an IC 50  value of 0.11\nμM, although the relatively low R 2  value (0.29) suggests\nsome uncertainty in the curve fit. Chlorzoxazone also demonstrated\nsubstantial BCL-2 inhibition, with an IC 50  of 0.65 μM,\nwhereas Vericiguat was less potent, with IC 50  values of\n17.71 μM. In terms of dose–response fitting, Vericiguat\nand Chlorzoxazone yielded the highest  R \n 2  values (0.87 and 0.64, respectively), while the fits for Lathyrol\nand Relugolix showed greater variability ( R \n 2  = 0.07 and 0.29, respectively). In the LN-18 glioma cell viability\nassay, Relugolix emerged as the most active compound, reducing cell\nviability with an IC 50  of 23.55 μM and a reasonable\ncurve fit ( R \n 2  = 0.66), in agreement with\nits biochemical activity in the TR-FRET assay. By contrast, Vericiguat\ndisplayed a much weaker antiproliferative effect, with an IC 50  of 261.1 μM, suggesting limited cytotoxic activity under the\ntested conditions. Chlorzoxazone and Lathyrol did not yield measurable\nIC 50  values within the tested concentration range, indicating\neither weak or inconsistent effects on cell viability.\nTaken\ntogether, these results suggest that Relugolix is a promising\ncompound, showing BCL-2 inhibitory activity in the TR-FRET assay and\na well-fitted antiproliferative response in cancer cells, although\nthe TR-FRET IC 50  estimate should be interpreted cautiously\ndue to its relatively low R 2  value.\n\nThis study demonstrates\nthe value of combining generative protein–ligand\ncomplex prediction, physics-based refinement, graph-based trajectory\nanalysis, and experimental testing in a unified BCL-2 drug-repurposing\nworkflow. Rather than relying on a single docking score or a purely\nligand-based screening strategy, we designed a multimodal pipeline\nin which NeuralPlexer was first used to generate ligand-specific BCL-2\ncomplex conformations, followed by docking refinement, MD simulations,\nMM/GBSA calculations, NRI, and  in vitro  validation.\nThis integrated design enabled prioritization of compounds on the\nbasis of structural plausibility, binding energetics, dynamic interaction\npatterns, and biological activity, thereby addressing several of the\nlimitations that commonly affect conventional repurposing studies.\nThe computational components of this workflow should be interpreted\nwithin the known strengths and limitations of each method. Docking\nand grid-based refinement provided an efficient means to evaluate\nlocal compatibility of NeuralPlexer-generated poses, but their scoring\nfunctions remain approximate and may not fully capture receptor flexibility,\nlong-time scale conformational adaptation, or solvent-mediated binding\neffects. MD simulations provided a dynamic description of ligand-bound\nBCL-2 complexes, but the 200 ns MD trajectories used here represent\nfinite sampling windows rather than exhaustive exploration of the\nconformational landscape. Similarly, MM/GBSA calculations were used\nto support relative prioritization rather than to estimate exact thermodynamic\nbinding affinities. For this reason, we did not base candidate selection\non any single computational metric. Instead, computational confidence\nwas assigned only when multiple orthogonal criteria converged, including\nstructural plausibility, local pose stability, ligand efficiency,\nMD/MM-GBSA behavior, NRI-derived similarity to the Venetoclax interaction\narchitecture, and subsequent biochemical and cellular validation.\nA key contribution of the present work is the application of NeuralPlexer\nas the entry point for large-scale, structure-guided screening against\nBCL-2. In contrast to conventional docking workflows, which often\nassume a largely rigid receptor and rank compounds primarily through\napproximate scoring functions, NeuralPlexer generates ligand-specific\ncomplex conformations in a manner that is more compatible with protein–ligand\nstructural cooperativity. In our study, this allowed the generation\nof bound-state models for a large FDA-approved drug set and yielded\n1294 structurally acceptable complexes for downstream analysis. Importantly,\nthe comparison with the co-crystallized Venetoclax–BCL-2 structure\nfurther supported the reliability of this approach, as NeuralPlexer\nreproduced the ligand pose with a binding-site RMSD of 0.35 Å\nwhile maintaining reasonable agreement at the protein level ( Figure  \n ). Taken together,\nthese findings support the use of NeuralPlexer-generated complexes\nas meaningful starting points for posthoc energetic and dynamical\nanalysis in repurposing studies.\nComparison of the cocrystallized BCL-2\nstructure of Venetoclax\n(PDB:  6O0K )\nwith NeuralPlexer predicted 3D structure. (A) The X-ray structure\nof BCL-2 is colored with cyan, while the NeuralPlexer predicted structure\nis represented with ice blue. (B) Venetoclax is colored as pink (X-ray\nconformer) and green (NeuralPlexer prediction).\nThe oral anticancer drug Doxifluridine had an average\nMM/GBSA score\nof −40.90 kcal/mol and a ligand efficiency score of −2.41\nkcal/mol. It also has antiangiogenic properties. Even if it is not\none of the best options, Doxifluridine is a worthwhile compound to\nlook into further in the context of BCL-2-driven tumors due to its\ncapacity to suppress angiogenesis and cause apoptosis in cancer cells.\nAfter Doxifluridine, the ligand efficiency of Fadrozole, an aromatase\ninhibitor, is −2.46 kcal/mol. It is approved to treat breast\ncancer and has an average MM/GBSA score of −41.88 kcal/mol.\nThis implies that Fadrozole might be able to maintain a relatively\nmodest molecular size while exhibiting a high affinity for binding\nto the BCL-2 protein binding site. Following Fadrozole, Chlorzoxazone\nhas a slightly better ligand efficiency score of −3.17 kcal/mol.\nThe possibility of the muscle relaxant chlorzoxazone as a BCL-2 inhibitor\nhas not been previously explored. With a ligand efficiency score of\n−2.88 kcal/mol, etofylline, a phosphodiesterase inhibitor used\nto treat cerebrovascular diseases, shows promise as a strong and selective\nBCL-2 inhibitor.\nThe initial NeuralPlexer-guided prioritization\nidentified several\ncompounds with favorable ligand–efficiency profiles, including\nEtofylline, Chlorzoxazone, Fadrozole, and Doxifluridine. These molecules\ndid not outperform Venetoclax in terms of average MM/GBSA binding\nfree energy, which remained most favorable for the reference inhibitor,\nbut they compared favorably after normalization for molecular size.\nThis distinction is important. Average MM/GBSA values tend to favor\nlarger ligands, whereas ligand efficiency can better reflect how effectively\na compound engages the target relative to its size. From a repurposing\nand medicinal chemistry perspective, smaller compounds with favorable\nsize-normalized binding characteristics may provide more tractable\nstarting points for optimization, particularly when the reference\nligand is relatively large, as in the case of Venetoclax. Accordingly,\nthe NeuralPlexer branch did not simply identify compounds with strong\nnominal binding scores, but also highlighted chemically compact scaffolds\nwith potentially useful binding efficiency. However, one of the most\nimportant findings of this study is that ligand efficiency alone was\ninsufficient to define the most promising candidate. The NRI analysis\nadded a mechanistic layer that changed how the top compounds were\ninterpreted. By extending NRI to include both protein residues and\nligand heavy atoms, we were able to move beyond static scoring and\nask how each ligand altered the dynamic interaction network of BCL-2\nduring MD simulations. Etofylline, for example, produced strong residue-level\ninfluence in the learned interaction maps and displayed a high normalized\nprotein–ligand interaction energy. Nevertheless, its contact\npattern appeared more spatially concentrated, particularly in the\nregion spanning residues 70–100, suggesting localized perturbation\nrather than broad stabilization of the bound protein state. Relugolix,\nin contrast, showed a contact pattern that more closely resembled\nthat of Venetoclax, even though some intraprotein interactions remained\nweaker. This observation is important because it suggests that the\nmost informative candidates are not necessarily those that maximize\none energy-derived metric, but those that best reproduce the interaction\narchitecture of a productive reference binder.\nThis mechanistic\ndistinction helps explain why Relugolix emerged\nas the most compelling hit from the NeuralPlexer/NRI arm of the workflow.\nRelugolix was not the strongest compound by ligand efficiency, yet\nit combined several favorable properties. A comparatively strong MM/GBSA\nprofile among the repurposed candidates, an NRI-derived interaction\npattern that closely paralleled Venetoclax, and consistent performance\nin subsequent biochemical and cellular assays. In other words, its\nprioritization was driven by multiparameter convergence rather than\nby any single ranking criterion. This is an important aspect of the\nworkflow, because one of the persistent problems in repurposing is\noverinterpretation of compounds selected only on the basis of docking,\nfree-energy estimates, or ligand-based prediction. Here, Relugolix\nadvanced because it was supported simultaneously by structure generation,\nenergetic refinement, dynamic network analysis, and experimental follow-up.\nClinically, relugolix is an orally active gonadotropin-releasing\nhormone (GnRH) receptor antagonist currently approved for the treatment\nof advanced prostate cancer and, in combination with estradiol and\nnorethindrone acetate, for the management of heavy menstrual bleeding\nassociated with uterine fibroids and moderate to severe pain associated\nwith endometriosis in premenopausal women.  As such, its established clinical utility is rooted in endocrine\nmodulation rather than in direct targeting of apoptotic machinery.\nThis is an important distinction, because the canonical clinical example\nof BCL-2 inhibition remains venetoclax, which is specifically approved\nas a BCL-2 inhibitor in hematologic malignancies. In this context,\nthe relevance of relugolix lies in the fact that, despite belonging\nto a therapeutically distinct class, it reproduced key Venetoclax-like\ninteraction features within the BCL-2 binding environment and demonstrated\nconvergent support across structure generation, energetic refinement,\nnetwork-level analysis, and downstream experimental evaluation. Taken\ntogether, these observations raise the possibility that relugolix\nmay possess off-target or secondary BCL-2-modulating potential, a\nhypothesis that warrants direct mechanistic validation through dedicated\ntarget-engagement and further apoptosis-focused studies.\nThe\nsecond major component of the study was the ML-guided prioritization\nbranch based on Chemprop and GNINA. This arm served an important role\nbecause it provided an orthogonal, ligand-centric perspective that\nwas independent of explicit NeuralPlexer complex generation. The Chemprop\nmodel achieved good predictive performance on known BCL-2 inhibitors,\nand the downstream GNINA/MM-GBSA/MetaCore pipeline identified several\nadditional compounds with plausible BCL-2 activity, including Lanraplenib,\nCobimetinib, Trazodone, Baricitinib, and Buspirone. While this branch\ndid not ultimately yield the experimentally strongest candidate, it\nstrengthened the overall study by reducing dependence on a single\ncomputational paradigm. In practical terms, the Chemprop-guided workflow\nfunctioned as an internal cross-check on the structure-based prioritization\nstrategy and demonstrated that the proposed screening framework is\nnot limited to one class of computational evidence. Notably, retrospective\napplication of the Chemprop model to the four experimentally validated\nNeuralPlexer-derived hits (Vericiguat, Lathyrol, Relugolix, and Chlorzoxazone)\nyielded predicted pIC 50  values of 6.69, 4.56, 6.50, and\n7.84, respectively. This discrepancy suggests that a purely ligand-based\nscreening strategy would likely have deprioritized these compounds,\nreinforcing the added value of the structure-based NeuralPlexer/NRI\narm in surfacing candidates that a conventional QSAR-type model might\noverlook.\nThe experimental results suggest that the NeuralPlexer-guided\nworkflow\nprovided the more effective route for identifying compounds that translated\ninto measurable biological activity. Among the compounds that showed\nthe strongest initial TR-FRET inhibition, most were selected from\nthe NeuralPlexer-based pipeline, and the final dose–response\nand cell-viability experiments further strengthened this trend, with\nRelugolix emerging as the potentially most convincing experimentally\npromising candidate and Chlorzoxazone and Vericiguat also showing\nbiochemical activity. By comparison, the ML-guided branch contributed\nat least one promising candidate, Trazodone, in the first-pass assay,\nindicating that ligand-based prediction can still enrich the screening\nset with biologically relevant compounds. However, the fact that the\nmost robustly supported compounds ultimately originated from the\nNeuralPlexer-guided arm suggests that explicit modeling of ligand-specific\nBCL-2 complexes, followed by MD/MM-GBSA and NRI-based dynamic analysis,\nmay provide higher practical precision for selecting experimentally\nsuccessful candidates in this target system.\nAt the same time,\nthe current results do not suggest that Relugolix\nis superior to Venetoclax. Venetoclax remained stronger in average\nMM/GBSA binding energetics and continues to outperform the repurposed\ncandidates under the tested conditions. The significance of Relugolix\nlies elsewhere: It represents an experimentally supported noncanonical\nscaffold that may engages BCL-2 and displays measurable cellular activity,\ndespite not having been developed as a BCL-2 inhibitor. For this reason,\nRelugolix is best viewed not as a replacement for Venetoclax, but\nas a promising starting point for medicinal chemistry optimization.\nIts experimentally supported partial activity, combined with a dynamic\ninteraction profile resembling that of Venetoclax, makes it a suitable\nscaffold for analog design aimed at strengthening target-specific\ncontacts, improving potency, and refining selectivity.\nThe overall\nfindings establish a useful conceptual framework for\nBCL-2-focused repurposing. The study shows that generative complex\nprediction can be combined with NRI-based dynamic analysis to prioritize\ncompounds not only by static fit or predicted affinity, but also by\ntheir ability to reproduce interaction patterns associated with a\nproductive reference inhibitor. The inclusion of an orthogonal Chemprop/GNINA\nbranch further improves the robustness of candidate selection, while\nexperimental follow-up provides an essential filter against computational\nfalse positives. In this context, the main advance of the work is\nnot simply the identification of Relugolix, but the demonstration\nthat structure-aware generation, dynamic graph inference, and ligand-based\nlearning can be integrated into a practical workflow for repurposing\nagainst BCL-2. This strategy should be transferable to other targets\nfor which ligand-induced conformational adaptation and residue-level\ncommunication are important determinants of productive binding.\nFrom a translational perspective, the most immediate next step\nis the optimization of Relugolix-derived analogs guided by the structural\nand dynamical information generated here. Compounds derived from the\nRelugolix scaffold could be designed to strengthen the interactions\nthat overlap with Venetoclax-like binding while compensating for the\nweaker intraprotein stabilization observed in the NRI and contact-map\nanalyses. More broadly, the present results support the idea that\nrepurposing campaigns should move beyond single-metric prioritization\nand toward integrated decision-making frameworks in which structure\ngeneration, dynamic interpretation, and experiment are explicitly\nlinked. Within that broader framework, Relugolix emerges as the most\nconvincing lead outcome of the current study and a reasonable starting\npoint for future BCL-2 inhibitor development.\n\nIn this study, we established\na multimodal computational-experimental\nworkflow for BCL-2-focused drug repurposing that integrates NeuralPlexer-based\nprotein–ligand complex generation, physics-based postprocessing,\nextended NRI, and an orthogonal ML-guided prioritization strategy.\nApplied to a library of 3094 FDA-approved compounds, this framework\nenabled large-scale generation of ligand-specific BCL-2 complex conformations,\nsystematic filtering of candidate binders, and dynamic comparison\nof prioritized compounds against the reference inhibitor Venetoclax.\nBy extending NRI to protein–ligand trajectories, we further\nmoved beyond static scoring and incorporated residue-level dynamic\ninteraction patterns into candidate selection, thereby providing an\nadditional mechanistic layer for prioritization.\nThe results\nhighlight two principal advances. First, the study\ndemonstrates that NeuralPlexer can serve as a practical entry point\nfor large-scale, structure-guided repurposing against BCL-2, generating\nchemically meaningful bound complexes suitable for downstream refinement\nand trajectory analysis. Second, the combination of NRI-based dynamic\ninference with conventional scoring metrics improved interpretation\nof the prioritized compounds by distinguishing candidates that merely\nscored favorably from those that more closely reproduced the interaction\narchitecture of a productive BCL-2 binder. In parallel, the Chemprop/GNINA\nbranch provided an orthogonal ligand-centric route for candidate prioritization,\nthereby increasing the robustness of the overall screening strategy\nand reducing dependence on a single computational paradigm. Among\nthe evaluated compounds, Relugolix emerged as the most compelling\nexperimentally supported hit. Its prioritization was not driven by\na single metric, but by convergence across multiple layers of evidence,\nincluding favorable binding energetics, a Venetoclax-like NRI interaction\nsignature, measurable BCL-2 inhibition in the TR-FRET assay, and antiproliferative\nactivity in LN-18 glioma cells. These findings support Relugolix as\na validated noncanonical BCL-2-interacting scaffold and a promising\nstarting point for future optimization.\nTaken together, this\nwork shows that generative complex prediction,\ngraph-based dynamic analysis, ligand-based ML, and experimental validation\ncan be combined into a coherent and effective framework for structure-guided\ndrug repurposing. Beyond the identification of Relugolix, the broader\nsignificance of the study lies in demonstrating a transferable strategy\nfor targets in which ligand-induced conformational adaptation and\nresidue-level communication are important determinants of productive\nbinding. We anticipate that this integrated workflow will be useful\nnot only for BCL-2-driven malignancies, but also for other therapeutic\ntargets where static scoring alone is insufficient for reliable candidate\nprioritization.\n\nIn this\nstudy, we generated high-quality conformations of the BCL-2 protein\nand the 3094 FDA-approved drugs from the DrugBank database using the\nprotein-ligand structure prediction tool, NeuralPlexer.  Using protein sequences and ligand molecules\nas input, the deep generative model NeuralPlexer can accurately predict\nthe 3D structures of protein–ligand complexes. The following\nare the main ways that NeuralPlexer  creates\nthe conformations: NeuralPlexer samples the 3D coordinates of the\nligand atoms iteratively by means of a diffusion process. This enables\nthe model to provide realistic and varied ligand conformations that\nmeet important biophysical requirements. Also, the model uses a hierarchical\nmultiscale architecture to improve the all-atom 3D coordinates step-by-step\nafter predicting residue-level contact maps. NeuralPlexer is able\nto capture the intricate structural cooperativity between the ligand\nand protein as a result. NeuralPlexer’s diffusion mechanism\nis engineered to take into account crucial biophysical limitations\nsuch as lengths, angles, and steric conflicts. This aids in the generation\nof physiologically plausible ligand conformations by the model. NeuralPlexer\ncan capture the conformational changes brought about by ligand binding\nsince it is trained to sample the protein in both its ligand-bound\nand ligand-free states. An ensemble of predicted ligand conformations\nis produced by NeuralPlexer, so they can be sorted and chosen according\nto confidence scores or other standards to determine the most promising\nbinding poses. The following particular settings that are applied\nwhen using NeuralPlexer: Produce 8 distinct forms for every protein–ligand\ncomplexes, divide the input into 4 segments to maximize memory utilization,\n40 diffusion stages are needed to produce the final conformations,\napply the sampling strategy of Langevin Simulated Annealing, find\nand manage connections between the ligand and the protein, then sort\nthe produced conformations according to the scores for confidence.\nNeuralPlexer produced high-quality conformations, which served as\nthe foundation for further computational investigation. This included\ncalculations of binding free energy, molecular docking, and MD simulations.\nThe success of this investigation was largely dependent on NeuralPlexer’s\ncapacity to precisely capture the intricate structural cooperativity\nbetween the BCL-2 protein and small molecule ligands.\nThe protein–ligand complex conformations\ngenerated by NeuralPlexer provide valuable structural insight. However,\nestimation of binding affinity is required to assess their therapeutic\npotential more rigorously. To enable this step, we developed a Python\nscript, which is available through our GitHub repository ( https://github.com/DurdagiLab/Neuralplexer_Ligand_Scoring ). The script first uses Maestro’s built-in tools to merge\nand prepare the protein–ligand complex structures generated\nby NeuralPlexer. A molecular docking grid is then defined by centering\nthe grid box on the ligand centroid, thereby specifying the binding\nsite for subsequent calculations. Glide is subsequently applied to\nestimate the Standard Precision (SP) docking score of each complex.  In addition, Glide’s “refine only”\nmode is used to locally optimize the ligand conformation starting\nfrom the NeuralPlexer-generated pose, without performing extensive\nresampling, and to report the corresponding optimized SP docking score\nfor each protein–ligand complex. By combining the robust scoring\nframework of Glide with the ligand-specific complex conformations\ngenerated by NeuralPlexer, this workflow enables efficient evaluation\nof the predicted binding affinities of candidate compounds. In this\ncontext, the “refine only” option provides a practical\nbalance between computational cost and scoring reliability by allowing\nlimited conformational adjustment while preserving the overall binding\ngeometry predicted by NeuralPlexer. Overall, the Python-based workflow\nstreamlines the screening process and supports high-throughput prioritization\nof large chemical libraries based on predicted binding affinity.\nAfter SP docking scores had been obtained, compounds\nwere grouped into individual  sdf  files according\nto their docking results, provided that the ligand conformation generated\nin Glide “refine only” mode did not deviate by more\nthan 2 Å from the original NeuralPlexer pose. These files were\nthen submitted to the MetaCore platform for screening with the “cancer\ntherapeutic activity prediction QSAR” module, which estimates\nthe likelihood that a compound possesses anticancer properties. In\nthis framework, therapeutic activity values (TAVs) are assigned by\ncomparing the query compounds with compounds known to have strong\nanticancer activity. \n , \n  TAV scores are normalized between\n0 and 1, and values above 0.5 are considered indicative of potential\nanticancer activity. The MetaCore cancer-QSAR model was developed\nusing a training set of 886 compounds and a defined descriptor set,\nachieving a sensitivity of 0.95, specificity of 0.92, accuracy of\n0.93, and a Matthews correlation coefficient (MCC) of 0.87. Validation\non an external test set of 167 compounds yielded a sensitivity of\n0.89, specificity of 0.83, accuracy of 0.86, and an MCC of 0.72. Based\non these criteria, a TAV threshold of 0.5 was applied for candidate\nselection, and among the compounds that passed this cutoff, the 10\nmolecules with the most favorable docking scores were retained for\nfurther analysis.\nProtein structures associated with compounds that\nsatisfied the binary QSAR cutoff of 0.5 were further refined using\nModRefiner,  a high-resolution protein\nstructure refinement method. ModRefiner can improve structural models\nat atomic resolution starting from Cα traces, backbone-only\nmodels, or full-atom structures. During refinement, both backbone\nand side-chain atoms remain flexible, and the conformational search\nis guided by a combination of physics-based and knowledge-based force\nfields. The method can also incorporate a reference structure to improve\nbackbone topology, side-chain orientation, and hydrogen-bonding geometry,\nthereby bringing the model closer to a physically realistic native-like\nstate. In addition, the stand-alone implementation supports  ab initio  full-atom relaxation without constraining the\nrefined model to the starting or reference structure. In the present\nstudy, this refinement step was included to correct potential local\nstructural inaccuracies, particularly proline-related geometric errors,\nin the NeuralPlexer-generated protein conformations.\nTo prepare the protein, the Protein Preparation Tool  was used. This included reassembling disulfide\nbonds, adding hydrogen atoms, forming zero-order bonds with metals,\nand allocating bond orders. With the help of PROPKA, the target protein’s\nresidues’ protonation states at physiological pH were assigned,\nand the OPLS3e force field  was used to\nreduce the side chain atoms.\nAll-atom MD simulations were subsequently performed, following\nthe preparation of the selected compounds and Venetoclax as the reference\nmolecule. The protein–ligand complex was positioned in a solvation\nbox with the TIP3P  water model in the\nsimulations, which were run using the Desmond.  A buffer zone measuring 10 Å in an orthorhombic box\nencircled the compound. Protonation states were assigned at pH 7.4\nusing Epik based on predicted p K \n a  values.\nThe system was neutralized with Na +  counterions and brought\nto physiological ionic strength (0.15 M NaCl). Ligand force-field\nparameters were assigned using the OPLS3e  force field, which provides atom-type coverage for a broad range\nof drug-like heterocycles without requiring custom parametrization.\nEach system was subjected to the Desmond default relaxation protocol\nprior to production simulation. This protocol consists of five sequential\nstages: (i) energy minimization of solvent and ions with solute restrained\n(2000 steps); (ii) 12 ps NVT MD at 10 K with solute restrained; (iii)\n12 ps NPT MD at 10 K with solute restrained; (iv) 12 ps NPT MD at\n310 K with solute restrained; and (v) 24 ps NPT MD at 310 K with all\nrestraints released, allowing the full system to equilibrate prior\nto production. Production MD was run for 200 ns per compound under\nthe  NPT  ensemble at 310 K and 1.01325 bar, controlled\nby the Nosé-Hoover thermostat  and\nMartyna–Tobias–Klein barostat,  respectively. Equations of motion were integrated using the RESPA\nintegrator with 2 fs (bonded), 2 fs (short-range nonbonded), and 6\nfs (long-range nonbonded) time steps. Short-range electrostatic and\nvan der Waals interactions were evaluated with a 9 Å cutoff;\nlong-range electrostatics were treated by the particle mesh Ewald\n(PME) method  under periodic boundary\nconditions. Each compound was simulated as 3 independent replicate\nsimulations. A total of 2000 frames were saved at equal intervals\n(one frame per 100 ps) from each 200 ns production trajectory for\nsubsequent analysis.\nBinding free energies were estimated using the MM/GBSA approach\nimplemented in the Prime  module of Maestro\n(Schrödinger). For each compound, 2000 evenly spaced snapshots\nwere extracted from the full 200 ns production trajectory (one frame\nper 100 ps), covering the entire simulation to ensure representative\nconformational sampling. The stability of binding-energy profiles\nacross the trajectory was used as an indicator of convergence. Three\nindependent MD replicates were performed. The VSGB 2.0  implicit solvation model was applied, with an\nexternal dielectric constant fixed at 80 and an internal dielectric\nconstant that varies between 1.0 and 4.0 under the OPLS3e force field,  consistent with standard practice for protein–ligand\nMM/GBSA calculations. Mean binding free energies and standard deviations\nare reported for each compound.\nThe MM/GBSA values were interpreted\nas comparative binding-energy estimates for candidate prioritization\nrather than as rigorous absolute binding free energies. Because MM/GBSA\nrelies on an implicit solvent approximation and does not fully account\nfor exhaustive conformational sampling, entropic contributions, or\nall solvent-mediated effects, the calculated values were used primarily\nto rank compounds within the same computational protocol and force-field\nenvironment. To reduce size-dependent bias, particularly in comparisons\nwith the large reference inhibitor Venetoclax, binding-energy values\nwere additionally evaluated in terms of ligand efficiency normalized\nby the number of non-hydrogen atoms. Accordingly, MM/GBSA was not\nused as an isolated determinant of activity, but as one component\nof a broader decision framework that also incorporated docking-pose\nstability, MD-derived interaction behavior, NRI-based dynamic coupling\npatterns, ML-guided prioritization, and experimental assay outcomes.\nIt is crucial to comprehend\nthe intricate relationships found in nature, particularly as they\nrelate to physical dynamical systems. Complex interacting atoms are\nproduced by MD simulations, and conventional examination of pairwise\nrelationships between residues requires the assumption of linear correlation.\nThe NRI model has an unsupervised variational graph-based autoencoder\narchitecture and have been proposed for the investigation of complex\ninteractions in dynamical systems, particularly in protein dynamics. \n , \n  By comprehending the nonlinear linkages as an interaction graph\nor latent representation  z \n \n ij \n , which represents the strength of the interaction between\nresidues  i  and  j , the NRI model\nseeks to recover dynamic trajectories. To evaluate protein–ligand\nenergy scores, we use the formula \n E z = 1 n l i g a n d ∑ i = 1 n r e s i d u e ∑ j = n r e s i d u e + 1 n r e s i d u e + n l i g a n d z i j δ i j \n where  n \n residue  and  n \n ligand  represent protein residue\nnumber and ligand heavy atom number,  z \n \n ij \n  is the learned interaction matrix and δ \n ij \n  = 1 for ∥ \n x \n \n \n i \n \n 0  –  \n x \n \n \n j \n \n 0 ∥\n≤ 12 Å and δ \n ij \n  = 0 for\n∥ \n x \n \n \n i \n \n 0  –  \n x \n \n \n j \n \n 0 ∥ > 12 Å. We followed\nthe\nliterature to select 12 Å as the threshold distance while evaluating\nthe energy score  E \n z . \n , \n  Higher energy scores mean that the ligand stabilizes the protein\nmore strongly. In other words,  E \n z  measures\nhow much ligand’s heavy atoms affect the residue dynamics on\naverage.\nIn this study, we extended the original NRI framework\nto reconstruct not only residue trajectories but also ligand trajectories,\nthereby enabling simultaneous learning of residue–residue,\nresidue-ligand, and ligand–ligand interactions within protein–ligand\ncomplexes. Using this approach, we analyzed 2000-frame MD trajectories\nfor a total of 11 BCL-2-ligand systems. The mathematical formulation\nand original implementation of the NRI model have been described previously. \n , \n  Here, we adapted the published source code to support protein–ligand\ntrajectory analysis. To facilitate automation and reproducibility,\nthe revised NRI workflow for protein–ligand complexes has been\nmade available through our GitHub repository ( https://github.com/DurdagiLab/Automate-Neural-relational-inference-NRI- ).\nThis implementation automatically organizes the trajectory\nfiles,\ninitializes model training, analyzes the resulting interaction maps,\nand exports the outputs as plots and CSV files.\nThe parallel ligand-based arm, combining a Chemprop\nmessage-passing neural network for pIC 50  prediction with\nGNINA docking, was included as an orthogonal validation strategy.\nChemprop operates directly on molecular graphs without relying on\npredefined fingerprints, allowing it to capture structure–activity\nrelationships from training data in a manner that is complementary\nto the physics-based scoring used in the primary structure-based pipeline.\nConvergence between the two independent arms, structure-based and\nligand-based, provides higher confidence that identified hits are\nnot artifacts of any single scoring scheme.\nA data set of 1067\ncompounds with experimentally determined pIC 50  values was\nobtained from the ChEMBL database. ( https://www.ebi.ac.uk/chembl/ ) Each molecule was represented by its canonical SMILES string, and\nactivity data reported as IC 50  values in molar units were\nconverted to pIC 50  using the relation pIC 50  =\n−log 10 (IC 50  [ M ]). To\nensure reproducibility, the data set was split into training (80%),\nvalidation (10%), and test (10%) sets using a fixed random seed (4523).\nModel development was carried out with Chemprop, a graph-based deep\nlearning framework that employs message-passing neural networks (MPNNs)\nto learn directly from molecular graph representations.  The model was configured for regression to predict\npIC 50  values, with key hyperparameters including a message-passing\ndepth of 3, two feed-forward hidden layers of size 300, and a batch\nsize of 50. Training was performed for 50 epochs using 5-fold cross-validation,\nwith mean squared error (MSE) as the optimization loss function and\na learning-rate schedule ranging from 1 × 10 –4  to 1 × 10 –3 . Model performance was monitored\nprimarily using root-mean-square error (RMSE), while  R \n 2  was also reported as an additional measure of predictive\naccuracy. The final trained model was then applied to a data set of\nFDA-approved drugs provided in CSV format with canonical SMILES strings,\nand Chemprop’s command–line interface was used to generate\npredicted pIC 50  values for downstream analysis. To support\nreproducibility, the SMILES-based data set splits and model checkpoints\nwere retained.\nMolecular docking was performed\nusing GNINA  an AutoDock Vina–based\nmethod that incorporates convolutional neural network (CNN) scoring\nto improve binding affinity prediction. The prepared crystal structure\nof human BCL-2 (PDB ID:  6O0K ) was used as the receptor, and ligand structures were\ngenerated from SMILES strings using LigPrep to obtain 3D conformers\nand protonation states appropriate for physiological pH. The binding\npocket was defined by a grid box centered on the known BCL-2 binding\nsite and sized to encompass the key interacting residues. Molecular\ndocking simulations were carried out with GNINA v1.0 using the default\nCNN scoring model and an exhaustiveness value of 8, generating multiple\nposes for each ligand. For each compound, the pose with the highest\nCNN affinity score was selected for further analysis, while both CNN-based\naffinity scores and conventional Vina scores were recorded. Compounds\nwith CNN affinity scores greater than 6.0 were retained as potential\nhigh-affinity binders. In addition, CNN-based ligand efficiency values\nwere calculated by normalizing the predicted binding score by the\nnumber of nonhydrogen atoms in each molecule. The resulting top-ranked\ncandidates were then prioritized for MM/GBSA binding free-energy calculations\nand MD simulations.\nThe BCL-2 TR-FRET Assay Kit (50222-1, BPS, USA) was used to measure\nthe inhibition of BCL-2 (B-cell lymphoma 2) binding to its ligand\nin the presence of BCL-2 inhibitory molecules in a homogeneous 96-well\nformat. The assay protocol for TR-FRET analysis was performed based\non the suggestions of the manufacturer. Briefly, a sample containing\nantihis terbium-labeled donor, dye-labeled streptavidin acceptor,\nBCL-2 protein, peptide ligand, and each inhibitor were incubated at\nroom temperature for 3 h. All samples and controls were examined in\nduplicate. Venetoclax was used as a positive control drug. After 3\nh of incubation at room temperature, the plates were read on a microplate\nreader (Varioskan Lux, Thermo Fisher, USA) with 340 ± 20 nm laser\nexcitation, a first emission filter at 620 ± 10 nm, and a second\nemission filter at 665 ± 10 nm. The percentage inhibitory activity\nof tested molecules was calculated by \n % a c t i v i t y = F R E T s − F R E T n e g F R E T p − F R E T n e g × 100 % \n where FRETs, FRETneg, and FRETp are sample\nFRET, negative control FRET, and positive control FRET, respectively.\nLN-18 (#CRL-2610, ATCC, USA) glioma cell line was used for cell\nculture experiments. Cells were seeded with high glucose Dulbecco’s\nModified Eagle Medium (DMEM) medium (Capricorne) supplemented with\n10% fetal bovine serum (FBS) (Gibco) and 1× penicillin/streptomycin\n(Multicell). 24 h prior to molecule treatment, 10,000 cells were seeded\ninto each well of 96-well cell culture plates. Values of half-maximal\ninhibitory concentration (IC 50 ) were determined by 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium\nbromide (MTT) cell proliferation assays. Different concentrations\nof molecules ranging between 100 μM and 1 nM were tested. Absorbance\nwas measured at 570 nm with microplate reader (Varioskan Lux, Thermo\nFisher, USA), and IC 50  values were calculated by dose–response\ninhibition curves and nonlinear regression analysis on GraphPad Prism\n8 software. For cell proliferation assays, we performed 24 h experiments\nperformed in triplicate.","source_license":"CC-BY-4.0","license_restricted":false}