| Review | Open Access |
|---|
Molecular Dynamics Investigation of Bcl-xL Interactions with Potential Cancer Inhibiting Compounds |
|---|
The Bcl-xL protein belongs to the Bcl-2 family of proteins. Bcl-xL is reportedly actively involved in cancer, along with various other factors. Since it has been revealed that this protein plays an active role in cancer, it can be a potential drug target to inhibit cancer activity. Its overexpression can lead to the evasion of apoptosis, allowing cancer cells to persist and proliferate uncontrollably. Hence, targeting Bcl-xL is a promising strategy in cancer therapy, while aiming to induce apoptosis in cancer cells and inhibiting their survival mechanism. In the current study, reported Bcl-xL inhibitors were evaluated for their stability, binding specificity, and interaction dynamics by performing molecular dynamics (MD) simulations. The compounds 24 and 54 revealed stable hydrogen bonding and hydrophobic interactions with active residues in the hydrophobic binding pocket. The binding sites were elucidated for the Bcl-xL protein. Exploring the molecular structural properties involved in binding provides mechanistic insights. This understanding aids in unraveling the protein-ligand interactions crucial for inhibiting anti-apoptotic proteins like Bcl-xL, offering potential for developing targeted inhibitors with anticancer properties.
B-cell lymphoma extra-large (Bcl-xL) is a protein that belongs to the Bcl-2 family of proteins. It is actively involved in the survival of various types of cancer and hematopoietic disorders [1, 2]. This family of proteins plays a crucial role in regulating apoptosis through the mitochondrial apoptotic pathway [3]. This particular protein encompasses both anti-apoptotic and pro-apoptotic members, sharing homology domains known as BH domains. Bcl-xL forms associations to create heterodimers, either with multi-domain pro-apoptotic members including BAX and BAK or with single-BH3 domain members such as BAD, BIM, and BID. The binding of these pro-apoptotic members results in the prevention of mitochondrial outer membrane permeability (MOMP), thereby impeding the release of cytochrome c and hindering the activation of caspases and counteracting apoptotic effects [4]. Apoptosis is intricately regulated by a balance between anti- and pro-apoptotic members. Any disturbance in their levels contributes to apoptotic dysfunction, associated with a cancer phenotype. Overexpression of Bcl-xL disrupts this balance in cancer cells that allows these cells to elude programmed cell death. This leads to uncontrolled proliferation, accumulation of mutations, and resistance to chemotherapeutic agents. Notably, Bcl-xL is observed to be overexpressed in various cancers including lung cancer, non-Hodgkin’s lymphoma, adenocarcinoma of reproductive tissues, pancreatic cancer, leukemia, and multiple myeloma [5–7]. The phenotypic expression of Bcl-xL in cancer cells is associated with multidrug resistance, undermining the efficacy of chemotherapeutic agents [8]. Bcl-xL also supports tumor progression by regulating other signaling pathways, such as PI3K/Akt and NF-κB pathways, and by binding to Beclin-1 to control autophagy mediated cell death [9].
The human Bcl-xL protein comprises four BH domains (BH 1-4) and includes a transmembrane region facilitating its span across the mitochondrial membrane. Structurally, all four domains exhibit an alpha-helical configuration with amphipathic characteristics [10]. In its folded conformation, domains 1-3 of this protein create a hydrophobic groove, serving as a binding site for BH3-only peptides [11]. These peptides bind to a hydrophobic pocket, adopting an alpha-helical conformation within the complex. Bcl-xL has shown its activity in three main defense pathways of cells, namely apoptosis, senescence, and autophagy. All these three pathways are interconnected and, depending upon the condition, they can be either beneficial or detrimental. The above protein has been identified to localize in the mitochondria, where it plays a crucial role in promoting cell survival. Interestingly, it promotes the ATP synthesis in mitochondria and regulates the Ca2+ flux from the endoplasmic reticulum. In this way, it contributes to energy balance which is essential for healthy aging [12].
Bcl-xL is actively involved in cancer biology, which makes it a potential therapeutic target against cancer. Studies have revealed that Bcl-xL, or more precisely Bcl-2 family of proteins, have multiple functions in addition to apoptosis, including the regulation of other cellular functions. Previously, it was believed that they have the ability to reduce the oncogenic properties and oppose apoptosis. It has been suggested in various studies that this protein has adopted new functions other than apoptosis [13–15]. It has a role in both in vitro and in vivo invasion of colorectal carcinoma, glioma [15], and breast carcinoma [13, 16]. In addition, several studies in pancreatic cancer showed an accelerated tumor level formation and growth due to it [17]. The versatile roles of this protein in cell proliferation makes it an attractive target for anticancer drug design. Given the flood of information regarding its role in disease endurance and movement, this research work is an effort to identify a distinguishing inhibitor against Bcl-xL to restrain malignant growth cell endurance. It is worth mentioning here that research on this protein and cancer treatment methodologies remains constantly evolving [18].
Computational approaches are utilized to identify potential strategies for the inhibition of Bcl-xL activity [19–26]. Molecular docking is a predictive technique that determines the optimal alignment of a ligand within a protein’s active site, forming a stable complex. MD simulation, on the other hand, is a computational method employed to study the physical movements of atoms and molecules. In MD simulations, atoms and molecules interact over a set duration, providing insight into the dynamic evolution of the system. These simulations are used to estimate structural, dynamic, and thermodynamic properties of the systems under observation [27].
In a previous study [11], various potential inhibitors of Bcl-xL were identified through virtual screening and molecular docking. However, the dynamic nature of intermolecular interactions and conformational flexibility cannot be comprehended completely through molecular docking alone. Hence, the current study is focused on an exclusive molecular dynamics study of previously reported potentially active compounds. The target is to gain mechanistic insights into ligand binding through the analysis of binding dynamics, bringing stability and causing key residue interactions. This would help to prioritize potential compounds for further optimization and experimental validation. Understanding these details can expand rational drug design efforts against Bcl-xL, promoting the development of potent compounds that help to restore apoptotic pathways in cancerous cells.
There is a variety of research currently being performed on some Bcl-xL inhibitors, including compound 54, based on MD analyses and binding affinity prediction. However, the current study extends the findings by performing comparative MD simulations across various reported inhibitors and provides novel mechanistic insights into the interaction dynamics allowing the prioritization of promising candidates. By highlighting key interacting residues and the differences in binding stability, this study contributes additional knowledge that may aid structure-based drug design of more effective Bcl-xL-targeted inhibitors.
The protein structure Bcl-xL was retrieved from the RCSB PDB website (http://www.rcsb.org) using a PDB ID 1R2D.
2.1. Preparation of Bcl-xL StructureThe protein structure was preprocessed using AutoDockTools (ADT). All missing hydrogen atoms were added, inappropriate bond orders were corrected, and molecular groups were aligned to achieve proper geometry and structural optimization, prior to further computational analysis [28]. During energy minimization, steric hindrances and strained bonds or angles were relaxed to obtain a stable molecular conformation. This allowed flexibility in heavy atoms up to 0.3 Å and up to 0.5 Å in others.
2.2. Ligand PreparationLigands were selected from the already reported research. Bcl-xL was reported to have the best binding affinity with ligands in the selected active site. A total of five ligands were selected from a previous study [11]. The compounds with numbers 24, 46, 50, 54, and 100 showed the best binding affinity with Bcl-xL protein. Hence, they were selected for the current study. ChemDraw was used to draw the structure of all these five ligands which were prepared before docking with the help of AutoDockTools [29].
2.3. Molecular Docking by AutoDock VinaMolecular docking was performed with AutoDock Vina [30]. The grid box covered the binding pocket using AutoGrid with a box of 40*40*40 Å along xyz directions to enclose the ligand. The box spacing was 0.3 Å and the grid center was designated in three dimensions (x, y, z) 23.048 Å, 36.224 Å, and 51.357 Å. The configuration file containing the information about protein and the ligand was used to execute AutoDock Vina. The pose with the lowest binding affinity value was selected for further analysis, which shows the most favorable free energy of the binding state. LigPlot and BIOVIA Discovery Studio Visualizer was used to visualize the 2D and 3D structures of ligands and target protein Bcl-xL, respectively [31, 32].
Figure 1. Schematic Outline for Molecular Docking
2.4. Molecular Dynamics SimulationsThe molecular dynamics (MD) simulation of Bcl-xL with the selected active compound was carried out using AMBER (Assisted Model Building with Energy Refinement) software (version 2021). The ligand was preprocessed with the help of antechamber. Different force fields were used for both ligands and the protein. The General Amber Force Field (GAFF) was used for ligand and ff14SB was used for protein. The tLeap module was used to record the topology of both the ligands and the protein. Sodium was added as counterion to the system to bring it into electrostatic neutrality. The TIP3P water model was used to solvate the system for geometric simplicity [33]. Energy minimization was performed for 2000 steps by executing 1500 steps of the steepest descent, followed by conjugate gradient steps. The cut off for non-bonded interactions was set to 8.0Å. The system was set for heating from 0 to 300 K for 25000 steps.
Equilibration was performed for a series of intervals at a constant temperature of 300K and a pressure of 1 atm. The production run was performed for 55ns with an isothermal-isobaric ensemble (T=300 K; P=1 atm). The periodic boundary conditions (PBC) were used to employ long-range electrostatic effects with particle-mesh-Ewald (PME) with a cut off non-bonded set of 8.0 Å. The SHAKE algorithm was used for hydrogen bonds. Production simulations for the ligands 54 and 24 were performed for a duration of 55 ns each. Root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), B-factor, and the radius of gyration of all the trajectories were analyzed with CPPTRAJ used in AMBER21.
Figure 2. Schematic Outline for MD Simulations
2.5. ADME InvestigationADME analysis was performed using SwissADME online software [34]. The five ligands were submitted to the software and a series of ADME factors were considered and reported.
Molecular docking studies were conducted to assess the binding affinity of ligands with protein active site residue (which in this case was Phe97) involved in the complex formation. Molecular docking was conducted to confirm that ligands fit properly into the protein pocket and form a protein-ligand complex. Autodock Vina gave optimal interaction and determined the lowest possible binding affinity which gave an ideal conformation from the 20 generated docked conformations for each of the five ligands. Subsequently, all selected docking compounds were subjected to docking into the Bcl-xL protein structure, with the obtained results presented in tables 1, 2, and 3, respectively.
Both hydrogen bonding and hydrophobic bonding play a very important role in the biological system [35]. The selected ligands interacted with the active site residue Phe97 and resulted in hydrophobic interaction, as the active site residue lies in the hydrophobic pocket of the ligand. Docking analysis indicated that Bcl-xL inhibitors occupied the hydrophobic groove within the BH3 domain, forming interactions with the key active-site residue Phe97. Among all five selected compounds, compound 54 displayed the lowest binding affinity or the lowest docking energy as (-11.2 kcal/mol) with Bcl-xL protein, whereas compound 50 showed the highest docking score (-2.5kcal/mol). The docking pose of all five compounds is shown in the following figures (3 to 7).
MD simulations aimed to explore the dynamic interaction between the target protein Bcl-xL and the inhibitor ligands 24 and 54. Keeping in view the time and computational resources required for extended MD simulations, only selected compounds were advanced at this stage. Top scoring compounds were prioritized for further analysis after docking. However, compound 46 was omitted because it is a structural analogue of compound 50 and its inclusion would not have provided additional mechanistic insight. The system’s stability was assessed through the analysis of the root mean square deviation (RMSD) of specific Phe97 atoms. The RMSD values of compound 24 did not exceed 2.5 Å. It did not exceed also in the last 400 ps of the simulation.
The average RMSD value for compound 24 with Bcl-xL was 1.8 Å. On the other hand, the average RMSD value for compound 54 with Bcl-xL was 1.3 Å. These results indicate the reliability of the data obtained in the final nanoseconds of the simulation, qualifying it for subsequent analysis. The time-dependent calculation of RMSD values for protein-ligand complexes offered comprehensive insights into the structural and conformational alterations of the protein [36].
Additionally, the root mean square fluctuation (RMSF) plot illustrates local variations within the protein chain [37]. The peaks of the plot point out the regions that fluctuate the most when undergoing MD simulations. Typically, the fluctuation phenomenon works in a way that the protein’s N terminal exhibit greater fluctuations when compared to the overall stability of the protein structure [38, 39]. In ligand24 Bcl-xL complex, RMSF shows that the binding region displays lower fluctuations in comparison to the freely moving side chains, where ligands do not bind [40] as shown in Figure 3b. The system displays high stability due to a strong binding between the ligand and the protein which is only possible when there are hydrophobic interactions or high hydrogen bonds. In this case, ligand54 binds to Phe97 which lies in a hydrophobic pocket of the protein.
The RMSF plot is used to characterize the local changes that occur along the chain of the protein [41]. The peaks of the plot indicate the regions that fluctuate the most during the MD simulation run. Observations typically reveal that the terminals of the protein exhibit greater fluctuation as compared to the remaining protein structure [38, 39]. Analyzing the ligand54 Bcl-xL complexes, RMSF indicates that the binding region experiences fewer fluctuations in contrast to the unbound side chains where ligands are not engaged, as shown in Figure 3b.
Figure 3. RMSD and RMSF Plots of Ligand 24-Bcl-xL Protein Complex
Figure 4. RMSD and RMSF Plots of Ligand 54-Bcl-xL Protein Complex
Figure 5. a. Docking Pose of Compound 24 b. Compound 24 2-Dimensional Structure. c. Compound 24 2-Dimensional Structure
The three-dimensional docking pose of compound 24 within the Bcl-xL binding pocket is shown below in Figure 5a. The protein is displayed in ribbon representation and the ligand is shown in bond form to highlight its orientation and key contacts. The corresponding two-dimensional interaction diagram is shown in Figure 5b, illustrating the specific residues involved in ligand binding. Notably, Phe97—located within the hydrophobic pocket—was identified as an important active-site residue contributing to the interaction profile [42]. The 2-dimensional pose docked compound 24 with the most active residues of the active site is shown below in Figure 5c. It highlights the hydrogen bond and its length. The selected Phe97 is also shown, which lies in the hydrophobic pocket of the protein along with other active site residues.
6. a. Docking Pose of Compound 54. b. Compound 54 2-Dimensional Structure. c. Compound 54 2-Dimensional Structure
Figure 6a shows the binding of compound 54 in the active site, indicating the binding residues of the protein. The protein is in the form of a ribbon and the ligand is in bond style, as shown in below Figure 6a. The docked compound 54 with the active residues of protein highlighting significant interactions is shown in Figure 6b. It is important to highlight that Phe97—located within the hydrophobic pocket—was identified as an important active site residue contributing to the interaction profile. The 2-dimensional pose docked compound 54 with the most active residues of the active site is shown below in Figure 6c. It highlights the hydrogen bond and its length. The selected Phe97 is also shown which lies in the hydrophobic pocket of the protein with other active site residues.
Table 1. Interaction Details of Active Compounds
|
Interactive Details of Active Compounds |
Interacting Atoms of the Ligand |
Interacting Atom of the Protein |
Distance (A˚) |
|---|---|---|---|
|
Compound 24 |
O |
ARG100:NH2 |
2.80 |
|
H |
ARG132:NE |
3.37 |
|
|
O |
GLN111:NE2 |
3.36 |
|
|
O |
ARG100:NH1 |
3.78 |
|
|
Compound 54 |
H |
LEU130:O |
3.40 |
|
H |
ARG132:NH1 |
3.29 |
|
|
H |
ARG132:NH2 |
3.92 |
|
|
H |
ARG132:NH1 |
3.46 |
|
|
H |
ARG132:NH2 |
3.85 |
|
|
O |
ARG103:NH2 |
3.11 |
|
|
O |
ARG132:NE |
3.94 |
|
|
O |
ARG132:NH1 |
3.00 |
|
|
O |
ARG132:NH2 |
3.45 |
|
|
O |
ARG132:NH1 |
3.30 |
|
|
O |
ARG132:NH2 |
3.43 |
In drug design, in silico prediction of druglikeness was performed through the calculation of ADME properties [43]. For this purpose, SwissADME was employed in this study. The physicochemical properties, lipophilicity, solubility, pharmacokinetics, and druglikeness were evaluated. Out of the five compounds presented in this study, compound 54 is the only one that fulfilled both Lipinski’s rule of five and Veber’s rule, making it the most druglike and favorable compound. The detail description of ADME properties is given below in Table 3.
Table 2. Binding Affinity and the Structure of Selected Compounds
|
Binding Affinity (kcal/mol) |
Compound 50 |
Compound 54 |
Binding Affinity (kcal/mol) |
||
|---|---|---|---|---|---|
|
-2.5 |
-11.2 |
||||
|
Binding Affinity (kcal/mol) |
Compound 46 |
Compound 100 |
Binding Affinity (kcal/mol) |
||
|
-10.7 |
-9.2 |
||||
|
Binding Affinity (kcal/mol) |
Compound 24 |
|
|||
|
-10.5 |
|
||||
Table 3. ADME Properties
|
|
24 |
46 |
50 |
54 |
100 |
|---|---|---|---|---|---|
|
Physicochemical Properties |
|||||
|
Formula |
C42H48O5 |
C34H14Br2O2 |
C37H20N4O4S |
C34H16O4 |
C28H21FN4O5S2 |
|
Molecular weight |
632.83 g/mol |
614.28 g/mol |
616.64 g/mol |
488.49 g/mol |
576.62 g/mol |
|
Num. heavy atoms |
47 |
38 |
46 |
38 |
40 |
|
Num. arom. heavy atoms |
24 |
32 |
40 |
32 |
24 |
|
Fraction Csp3 |
0.38 |
0.00 |
0.00 |
0.00 |
0.07 |
|
Num. rotatable bonds |
14 |
0 |
2 |
0 |
8 |
|
Num. H-bond acceptors |
5 |
2 |
6 |
4 |
7 |
|
Num. H-bond donors |
2 |
0 |
0 |
2 |
2 |
|
Molar Refractivity |
192.24 |
160.79 |
184.83 |
149.43 |
154.19 |
|
TPSA |
83.83 Ų |
34.14 Ų |
122.09 Ų |
74.60 Ų |
167.13 Ų |
|
Lipophilicity |
|||||
|
Consensus Log Po/w |
8.42 |
8.02 |
5.14 |
6.10 |
4.82 |
|
Water Solubility |
|||||
|
Log S (ESOL) |
-9.98 |
-10.20 |
-9.20 |
-8.10 |
-7.50 |
|
Solubility |
6.67e-08 mg/ml; 1.05e-10 mol/l |
3.88e-08 mg/ml ; 6.31e-11 mol/l |
3.93e-07 mg/ml ; 6.37e-10 mol/l |
3.85e-06 mg/ml; 7.88e-09 mol/l |
1.82e-05 mg/ml; 3.15e-08 mol/l |
|
Class |
Poorly soluble |
Insoluble |
Poorly soluble |
Poorly soluble |
Poorly soluble |
|
Pharmacokinetics |
|||||
|
GI absorption |
Low |
Low |
Low |
Low |
Low |
|
BBB permeant |
No |
No |
No |
No |
No |
|
P-gp substrate |
Yes |
Yes |
Yes |
No |
No |
|
CYP1A2 inhibitor |
Yes |
No |
No |
No |
No |
|
CYP2C19 inhibitor |
Yes |
No |
No |
No |
Yes |
|
CYP2C9 inhibitor |
No |
No |
No |
No |
Yes |
|
CYP2D6 inhibitor |
No |
No |
No |
No |
No |
|
CYP3A4 inhibitor |
No |
No |
No |
No |
No |
|
Druglikeness |
|||||
|
Lipinski |
2 violations: MW>500, MLOGP>4.15 |
2 violations: MW>500, MLOGP>4.15 |
2 violations: MW>500, MLOGP>4.15 |
0 violation |
1 violation: MW>500 |
|
Veber |
1 violation: Rotors>10 |
Yes |
Yes |
Yes |
1 violation: TPSA>140 |
|
Bioavailability Score |
0.56 |
0.17 |
0.17 |
0.55 |
0.55 |
Although previous studies have reported the binding affinity and MD simulations of compound 54 with Bcl-xL [11], the current study provides comparative insights across multiple reported inhibitors, allowing a more comprehensive evaluation of their relative stability and binding behavior. The MD simulation revealed differences in interaction patterns, hydrogen bonding, and active site dynamics that were not highlighted in earlier work. By systematically analyzing multiple compounds, this study prioritizes candidates with the maximum potential for further experimental validation, hence contributing to rational drug design.
This study extends the findings of previous studies by performing MD simulations of different Bcl-xL inhibitors. Although literature provides substantial evidence that Bcl-xL actively contributes to cancer cell survival, the discovery of novel inhibitors remains a pressing need. In this study, the best-docked protein-inhibitor complex was subjected to MD simulations to analyze the system for further stability. The set of MD simulations was performed to check whether the selected ligand is stable enough to be used as an inhibitor, when placed in a dynamic environment for a particular time in the presence of the target protein. The RMSD and RMSF graphs were generated to analyze the stability of protein when it is bound to the inhibitors. The simulation was run for 55 ns and its maximum value recorded for both inhibitors was 2.5 Å for inhibitor24 and 2.0 Å for inhibitor54, which shows that the prepared system is stable. The RMSF plot also predicted that the system showed less fluctuation in sites where inhibitors bind, as compared to the free chain where inhibitors do not bind due to hydrophobic interactions with the Phe97 active residue on inhibitors and protein. The results verify that the prepared system shows stability when subjected to MD simulation and the selected ligand can be used as a potent compound against Bcl-xL to inhibit cancer activity. This study provides valuable insights into stability, protein-ligand interactions, and key residue contributions for a clear understanding of anti-apoptotic protein inhibition. Such insights hold potential for further experimental validation and the development of targeted Bcl-xL inhibitors, offering promise as anticancer agents. Further, this study offers a mechanistic basis to guide the rational design of more effective Bcl-xL-targeted anticancer agents.
Asma Abro: conceptualization, writing-original draft, formal analysis. Sumra Wajid Abbasi: data curation, visualization, writing-review and & editing. Muazma Nabi: formal analysis, writing-review and & editing
The authors of the manuscript have no financial or non-financial conflict of interest in the
subject matter or materials discussed in this manuscript
The data associated with this study will be provided by the corresponding author upon request.
This research received no external funding.
During the preparation of this work, the authors used ChatGPT in order to improve the flow and language. After using this tool/service, the authors reviewed and edited the content as needed and takes full responsibility (legal, moral, etc.) for the content of the manuscript.