Repurposing of FDA-approved drugs as autophagy inhibitors in tumor cells

Abstract Autophagy and apoptosis are the two crucial processes of programmed cell death found in all eukaryotic cells; however, the elevated physiological stress in the tumor microenvironment leads to uncontrolled up-regulation in the process of autophagy. Available literatures suggest that inhibiting up-regulated autophagy in the cancerous cells may lead to the apoptosis and thereby culminate to tumor clearance. Several studies have been performed to design autophagy-inhibitors using either Beclin-1 or Bcl-2 as a target in isolation. However, to overcome the constraints of the availability of small and potent autophagy inhibitors, we have attempted extensive computational approach of repurposing the FDA-approved drugs from the ZINC database in order to inhibit the interaction between the Beclin1 and Bcl-2. Out of 1565 FDA-approved drugs used in our computational work, we sorted the drugs Ponatinib, Simeprevir, and Nilotinib through various methods viz. molecular docking, Lipinski’s filter, MD simulation and MM/PBSA, and we found these aforementioned drugs to show good binding energy and favorable interaction with the BH3 domain of Beclin1. We anticipate from our computational results that these drugs may become potent candidates to inhibit autophagy and exhibit the apoptosis in the tumor microenvironment and combat the current limitation of potent autophagy inhibitors; however, to substantiate our in-silico results, further experimental validations of these drug molecules are currently in progress. Communicated by Ramaswamy H. Sarma


Introduction
Autophagy and apoptosis are the two crucial processes of programmed cell death found in all eukaryotic cells. Autophagy serves as an important conserved cellular process of cell death, both, constitutively as a recycling pathway for aged proteins as well as an up-regulated stress response mechanism. In cancer, the physiological stress goes beyond the threshold value which leads to uncontrolled up-regulated autophagy and the tumor cells use this elevated process of autophagy as a savior for meeting its exaggerated growth requirements. The mechanism behind this beyond-threshold physiological stress-dependent up-regulated autophagy involves an interaction of a key autophagic protein Beclin1 through its BH3 domain with an anti-apoptotic protein Bcl-2.
During physiological stress such as hypoxia and starvation, the cells adapt to a class II type of programmed cell death, referred as, autophagy and this adaptive autophagy behaves as a tumor suppressive process. As the physiological stress gets prolonged and upon the chemotherapeutic interventions, the tumor cells start using the process of autophagy as a mean of tumor promoting pathway, wherein, the induced autophagy generates the metabolic and biosynthetic catabolites in effort to meet the nutritional requirement imposed by the extensively growing tumor cells. Moreover, a plethora of studies have been reported where the process of adaptive autophagy exhibited in the tumor cells transforms into apoptosis as the stress goes beyond the threshold value (Sun et al., 2010).
Human Beclin1 is an essential and central component of class III PI3K complex formed during the initiation of autophagy. It is composed of 450 amino acids and containing three domains, out of which, one is BH3 domain (amino acid residues 108-127) which has been found to interact with the Bcell lymphoma 2 (Bcl-2), an anti-apoptotic protein and therefore gives an evidence of molecular interaction between autophagy and apoptosis (Liang et al., 1998;He and Levine, 2010). These two proteins have well-established roles in the development of cancers such as monoallelic deletion of human Beclin1 gene has been found in 40% of human breast cancers reported at isolated instances, and overexpression of the antiapoptotic protein Bcl-2 has been shown as a result of chemotherapeutic resistance in different forms of cancers Rashad et al., 2019). Therefore, the interaction of both Beclin1 and Bcl-2 may provide a crucial target for combatting chemotherapeutic resistance and promoting tumor clearance.
In stressed conditions, where the cell tries to suppress the tumor microenvironment, the Beclin1-Bcl-2 interaction may reduce the formation of Beclin1/Vps34 complex thereby reducing the enzymatic action of Beclin1/class III PI3K complex and ultimately inhibiting the Beclin1-dependent autophagy (Liang et al., 1998). As there are other proteins bound to Beclin1 in the Beclin1/Vps34 complex which altogether contribute in promoting autophagy, hence, Beclin1-Bcl-2 interaction-dependent inhibition of autophagy finally culminates to a controlled autophagy in that particular tumor-suppressive microenvironment (Pattingre et al., 2005;Pattingre & Levine, 2006). As the stress prolongs, the interaction between Beclin1 and Bcl-2 gets weakened and the inhibitory action of Bcl-2 on Beclin1-induced autophagy gets abolished, which in turn leads to excessive autophagy (Maiuri et al., 2007). However, a previous study has demonstrated that BH3 mimetics also interact with the Beclin1 BH3 domain and inhibit the Beclin1-Bcl-2 interaction, thereby increasing autophagic activity and finally lead to inhibit the tumor growth by suppressing the tumor microenvironment (Akar et al., 2008;Lian et al., 2010;Oltersdorf et al., 2005;Wang et al., 2017). A potential regulatory mechanism behind the inhibited tumor growth may involve a process wherein, upon the dissociation of Bcl-2 from Beclin1 takes place, the dissociated Bcl-2 gets phosphorylated by c-Jun N-terminal kinase 1 (JNK1) and as the concentration of phosphorylated Bcl-2 rises, it does not bind to the BH3 domain of the pro-apoptotic proteins which ultimately triggers to successive apoptosis and culminates to tumor clearance and inhibit the tumor growth .
Due to the limited availability of potent small molecule autophagy inhibitors, the further research in inhibiting tumor growth and tumor clearance is constrained. In this study, we used in silico approach for integrating the molecular docking, Molecular Dynamics simulation and MM/PBSA method to screen the FDA-approved drugs for their repurposing ability to inhibit the cellular autophagy in tumor cells. Moreover, the dual key attributes of the anti-apoptotic protein Bcl-2 in inhibiting both autophagy and apoptosis, make this protein as a potential chemotherapeutic target. Therefore, our current study may act as a new insightful approach and help in further investigating our suggested pharmacological molecules for their ability to behave as an autophagy-inhibitor(s), which may further be used as a tumor-sensitizing treatment to combat the chemotherapeutic resistance and subsequently to inhibit the tumor growth by targeting this Beclin1-Bcl-2 complex.

Retrieval and preparation of protein
For the purpose of current study, the experimentally determined structure of BH3 domain of Beclin1 complexed with the anti-apoptotic protein Bcl-2 (PDB ID: 5VAX) with very high resolution of 2.00 Å was selected and retrieved from RCSB Protein Data Bank and except for the sequence of amino acid residues 109-129 of chain E of Beclin1_BH3 domain, all sequences present in the structure 5VAX were removed. Then after, the remaining structure (Beclin1_BH3 domain) was energy minimized with the help of the GROMOS96 implementation of Swiss-PdbViewer (Guex & Peitsch, 1997). In this study, we purposefully removed the amino acid residue phosphothreonine (TPO) at 108 th position so as to improve the dissociation of Beclin1_BH3 domain and Bcl-2 (Maejima et al., 2013). After energy minimization, the minimized structure of the protein was checked for the polar hydrogen based on the atom types, avoidance of any steric clashes and hydrogen bond formation, and then assigned for the partial atomic charges, atomic solvation parameters and fragmental volumes and then further optimized using the established protocol of AutoDockTools4.2 (ADT 4.2). The default configuration for protonation of receptor (net charge 0) was employed during the process of molecular docking. Finally, the pdb format of the protein was converted into the pdbqt format using the same tool (Morris et al., 2009).

Retrieval and preparation of ligands
For this study, the 3-D SDF structures of the 1565 FDAapproved drugs were downloaded from the ZINC database, a public repository of commercially-available compounds for virtual screening (Irwin et al., 2012). The SDF structures of the ligands were subjected to energy minimization using 1 st order common algorithms; i.e. Steepest Descent (SD) and Conjugated Gradients (CONJ), followed by addition of hydrogen and assignment of Gasteiger charges using Assisted Model Building with Energy Refinement (AMBER) force field (Wang et al., 2006) using UCSF Chimera, an extensible molecular modelling system. After completing the preparation steps, the prepared ligands were converted into the mol2 format (Ef et al., 2004).

Docking used in drug screening study
During the process of blind docking, the protein molecule was considered as a rigid molecule while the ligand was as a flexible entity. Defining the grid parameter file (GPF) is one of the crucial stages among all the steps to be followed in the docking procedure. The aforementioned parameter file was generated through the ADT 4.2 tool. This tool follows the Lamarckian Genetic Algorithm (LGA) which uses the AMBER force field in order to run the docking between the receptor and the ligand (Morris et al., 2009). To prepare the GPF, a 3-D grid was defined with À12.46, 7.208 and À12.938 as X-, Y-and Z-coordinates respectively so as to evaluate the favourable binding energy between the protein and the ligand. Additionally, in the GPF, the number of grid points was set as 48, 74, and 84 in X-, Y-, and Z-coordinates respectively and the spacing value was set as 0.375 Å. Further, in silico drug screening among the selected 1565 FDA-approved drugs were performed using AutoDock Vina with default configuration parameters (Trott & Olson, 2009) with exhaustiveness value of 8 to evaluate potentially effective drug candidates based on the favourable binding energy of the individual drug.

Assessment of druglikeness of the ligands
Among the 1565 FDA-approved drugs, the top 25 drugs were checked for their role in the cellular autophagy by going through various documented literature reports. By this, we could sort four molecules among the 25 drugs and then we assessed these four drugs (For 3-D structures of these four drugs, refer Supplementary Table S1) for their ability to follow Lipinski's rule of five (Lipinski, 2004) using the online SwissADME analysis server (Daina et al., 2017). The SMILES string of all the four compounds was retrieved from the input available on the PubChem database. This Lipinski's rule of five defines a criterion according to which a newly discovered orally effective drug should comply to four out of five suggested criteria; which are: molecular mass, lipophilicity of the molecule, number of hydrogen donors and acceptors and molar refractivity. Furthermore, this online SwissADME analysis server was also used to predict the pharmacokinetic properties of a drug such as absorption, distribution, metabolism, and excretion (ADME).

Molecular dynamics (MD) simulation study
After sorting the drugs on the basis of already available study reports of involvement of the drugs in cellular autophagy and by subjecting them to Lipinski's filter analysis, the chosen 4 drugs were subjected to perform their molecular dynamics simulation studies. This is a powerful computational tool to investigate the interactions between the biomolecules using the approach of atomic resolution based on the assumptions of Newtonian physics and thus giving an insight on dynamics of the biomolecule at distinct time scale (Durrant & McCammon, 2011;Borkotoky et al., 2016). The calculations for classical MD simulation were performed for the prepared apoprotein and the selected docked proteinligand complexes through the GROningen MAchine for Chemical Simulation (GROMACS 18.8) (Lemkul, 2019) MD simulation package using GROMOS96 as the force field (Oostenbrink et al., 2004). The PRODRG server was employed as a parameter generator the and topology files builder for the ligands (Sch€ uttelkopf & Van Aalten, 2004). Then after, the biomolecular complex was centred in a dodecahedron periodic box with a distance of 1.0 nm from the edge of the box as satisfying the periodic boundary conditions. The defined area was then subjected to solvation using the simple point charge (SPC), a three-point water model so that the biomolecular complex within the box may get appropriately filled with water molecules. Sodium ions or chloride ions were added in the system in effort to make the overall charge of protein and ligand complex to a neutral condition. After solvation and making the system electroneutral, any steric clashes present within the system was removed by relaxing the structure by subjecting the system to energy minimization. So, the energy minimization of the system was performed with a convergence threshold of 1000 kJ mol À1 nm À1 for 50,000 steps using the steepest descent method. Then after, the solvent and ions around the protein were equilibrated using the methods of modified Berendsen thermostat coupling (Bussi et al., 2007) and Parrinello-Rahman pressure coupling (Marton ak et al., 2003) at temperature 300 K with the coupling constant set at 0.1 ps and pressure 1.0 bar with the coupling constant set at 2.0 ps under an isothermal-isochoric ensemble (NVT ensemble or constant Number of particles, Volume and Temperature) and isothermal-isobaric ensemble (NPT ensemble or constant Number of particles, Pressure and Temperature) respectively. The particle mesh Ewald (PME) algorithm (Kawata & Nagashima, 2001) was used in the system to calculate the long-range electrostatic interactions, however, the cut-off distance for short-range electrostatic (rcoulomb) and Van der Waals (rvdW) interactions was set at 1.0 nm. The bond parameters were restricted using the LINear Constraint Solver (LINCS) algorithm (Hess et al., 1997). The time step (dt) was set to 0.002 ps. Each of the equilibration phases was carried out for 1000 ps. Hereafter, the production of MD simulation was allowed to run for 100 ns with the time step of 0.002 ps for both the apoprotein and the protein-ligand complexes. The generated trajectories were preserved for every 10 ps for further analysis.

MM/PBSA for free energy calculations
After performing the MD simulation of the 4 drugs, the same set of drugs were subjected to calculation of the contribution energy by each residue in individual peptide-drug complex. A highly efficient computational method to compute the free energy interaction between the biomolecules is molecular mechanics/Poisson-Boltzmann surface area (MM/ PBSA) method which is widely carried out in combination with MD simulation. This method is based on a Poisson-Boltzmann equation (PBE) model which summarises the association of electrostatic potential with the solute's and the solvent's dielectric properties, along with the solution's ionic strength and the reach of the ions to the interior of the solution and the allocation of the solute's atomic partial charges (Baker et al., 2001). Hence, the PBE model is a type of implicit solvent method which is required to estimate the interaction between the polar solute and the solvent by using the solvent as a simple dielectric continuum model (Wagoner & Baker, 2006). The estimated free energy calculated by a three-step calculation of the MM/PBSA method which basically comprises to three energy components i.e. change in the potential energy, desolvation of species and configurational entropy of the complex in the gas phase. In the current study, we used an open-source tool 'g_mmpbsa' which employs the simulation trajectories to estimate the energy due to the binding interaction of the protein-ligand complexes. This tool uses a theory to estimate the enthalpic components of the MM/PBSA interaction which is explained as follows; Where, G complex is the total free energy of the biomolecular complex and G protein and G ligand are the free energies of the protein and ligand in the solvent, respectively.
Where, E MM is the average potential energy based on molecular mechanics in a vacuum. E bonded and E nonbonded are the bonded and nonbonded interactions respectively. E vdW and E elec are the nonbonded interactions and are represented for Van der Waals and Electrostatic interactions respectively.
Where, G x is the free energy for each individual specie (protein/ligand/protein-ligand complex). -TDS term represents change in the conformational entropy upon the drug binding to the receptor. This entropy term is not taken into consideration due to its tedious and time-consuming calculation. G solvation denotes the solvation free energy. The components G polar and G nonpolar stand for the polar and nonpolar contributions to the free energy of the solvation. In order to calculate the G polar term, the solute dielectric constant (e in ) was set at 2, the value of solvent dielectric constant was set at 80, and rest of the parameters were kept at its default settings in order to calculate the solvation free energy using pbsa method. The 'g_mmpbsa' tool is also used to decompose the total binding energy (BE) on a per residue (x) basis which can be calculated using following equation; are the energies of i th atom from the residue 'x' in complexed and free forms respectively. 'n' is the total number of atoms in a residue.

Analyses and visualization
The docking results were analysed using the AutoDock Vina and the interactions between the protein-ligand complexes were visualized through LigPlot þ v.2.1, which is a graphical user interface for the LIGPLOT program (Laskowski & Swindells, 2011) and the Protein-Ligand Interaction Profiler (PLIP) (Salentin et al., 2015). For general visualization purposes, PyMOL2.3.3 was used. After the MD run production, the analyses of the structural properties; such as, root-meansquare deviation (RMSD), root mean square fluctuation (RMSF) and radius of gyration (Rg) of the protein and its bound ligand complexes were performed with the help of the trajectory files generated through the 'g_rms', 'g_rmsf', and 'g_gyrate' utilities respectively, the in-built functions of GROMACS 18.8. The number of hydrogen bonds formed between the protein and the ligand was assessed, which is based on bond length smaller than 3.6 Å distance between the protein and the ligand and the angle between the donor-hydrogen and acceptor larger than 90 through the 'g_hbond' utility of GROMACS 18.8. The data for threedimensional backbone of the RMSD, the RMSF, the Rg, and number of H-bonds were plotted using the QtGrace version of GRaphing Advanced Computation and Exploration (GRACE) program.

Molecular docking
In the current study, 1565 FDA-approved drugs were subjected to molecular docking with the Beclin1_BH3 domain as a receptor using the AutoDock Vina based VEGA ZZ server as a docking platform. After docking, the top 25 compounds among the 1565 FDA-approved drugs were sorted on the basis of the VEGA ZZ-generated highest negative binding energy with their best binding model along with their 3-D representative structures. Then after, further sorting was done on the basis of available documented literature reports of involvement of these 25 drugs into the process of cellular autophagy. Based on this, we could sort the 4 complexes among the top 25 drugs. Although, we got more than 4 drugs in the top 25 drugs list to be involved in autophagy, but most of the drugs among the top 25 drugs list were not able to be submitted successfully in the PRODRG server during the ligand preparation step in Molecular Dynamics simulation. To solve this problem, we downloaded the 3-D SDF format of the top 25 drugs from the PubChem database (Wang et al., 2009) and among those, only 4 drugs were able to be successfully submitted in the PRODRG server during the ligand preparation step in Molecular Dynamics simulation. Therefore, we are showing only those 4 complexes among the top 25 drugs list which got a smooth submission in ligand preparation step during MD simulation and also are involved in the process of autophagy. Further, these 4 drugs were assessed for their ability to pass through Lipinski's filter and then after the same set of drugs were subjected to MD simulation studies, followed by calculation of binding energy through the MM/PBSA method. In our current study, after following the whole aforementioned protocol of sorting the 1565 drugs, we are focusing on giving the docking result of only the chosen 4 drugs (Figure 1). The binding energies of all the docked complexes were calculated through VEGA ZZ server and all the ligands showed a good binding energy. However, among the chosen 4 hits, Digoxin was found to be interacting with the peptide with the highest negative binding energy of À6.9 kcal/mol, followed by Nilotinib (-6.7 kcal/mol), Ponatinib (-6.5 kcal/ mol), and Simeprevir (-6.5 kcal/mol). Furthermore, out of an interesting observation, it was found that, mostly, the C-terminal residues of the BH3 domain of the Beclin1 were involved in the hydrogen bond formation with all of the aforementioned ligands. In case of the complex Digoxin (PubChem CID: 2724385, C 41 H 64 O 14 ), the residue Ser 127 of the Beclin1_BH3 domain formed the hydrogen bond with the ligand. However, other residues such as Asp 121 , Asp 124 , Gly 128 , and Lys 117 of the peptide were also found to form hydrogen bonds with the ligand Digoxin (Figure 2A). Similarly, the residues Asp 121 and Asp 124 of the peptide were involved in the formation of hydrogen bond with the ligand Nilotinib (PubChem CID: 644241, C 28 H 22 F 3 N 7 O) ( Figure 2B). Moreover, in case of complex Ponatinib (PubChem CID: 24826799, C 29 H 27 F 3 N 6 O), none of the residues were involved in the hydrogen bonds formation with the peptide ( Figure 2C), whereas, the complex Simeprevir (PubChem CID: 24873435, C 38 H 47 N 5 O 7 S 2 ) showed the hydrogen bond formation with the residues Lys117 and Asp 121 of the peptide ( Figure 2D). Furthermore, other than hydrogen bonds, there were also hydrophobic as well as a few of non-covalent (salt bridges and pi-Stacking) interactions observed between the peptide and the ligand (Table 1). Noticeably, the residue Phe123 of the peptide was found to be more consistent in making hydrophobic interaction with almost each of the four ligands. The drug Digoxin was found to be interacting with the residues Lys117 and Phe123 of the peptide through hydrophobic interaction. Similarly, Nilotinib formed hydrophobic interaction with the residues Leu116 and Asp124 of the peptide. The ligand Ponatinib was also involved in the hydrophobic interaction with the residues Leu116, Lys117, and Phe123 of the peptide and the drug Simeprevir showed maximum number of hydrophobic interactions with the residues Leu116, Lys117, Thr119, Phe123, and Asp124 of the BH3 domain of the protein Beclin1. A very few non-covalent interactions were also observed such as salt bridge formation between the drug Digoxin and Lys117 of the peptide and the pi-Stacking interaction between the drug Nilotinib and Phe123 of the peptide. However, among all the 4 ligands, only Ponatinib was found to form halogen bond through its Fluorine atoms with the Met109 of the BH3 domain of the Beclin1 protein.
Hydrogen bond along with hydrophobic interaction and other non-covalent interactions (such as salt bridge and pi-Stacking interaction) play important role in determining the favourable interactions between the protein and the ligand. However, the presence of hydrogen bonds is not necessarily required to build a favourable association between the biomolecular complexes, as according to a reported literature, removal of hydrogen-forming group leads to a strong interaction between a protein and a ligand (Zhao & Huang, 2011).
In a recent study, scientists have reported that the deathassociated protein kinase 2 (DAPK 2) enzyme binds to Beclin1 and phosphorylates the 119 th positioned Threonine residue of Beclin1 which leads to dissociation of the Beclin1-Bcl-2 binding and ultimately leads to activation of autophagy (Shiloh et al., 2018). In case of Simeprevir, a hydrophobic interaction was observed with the Thr119 of the BH3 domain of Beclin1 protein which suggests that owing to this hydrophobic interaction, this ligand Simeprevir may act as a potent candidate to inhibit the interaction between Beclin1 and Bcl-2, which will lead to inhibit autophagy which may further culminate to apoptosis of the tumor cells.
The overall molecular docking results reveal an observation of contribution of some common amino acid residues such as Lys117, Asp121, and Asp124 in hydrogen bonding, hydrophobic interactions and other non-covalent interactions in almost all protein-ligand complexes; signifying that the aforementioned amino acids may account for a critical importance in building a favourable interaction with all the 4 hit ligands. Then after, these drugs were subjected for further analyses.

Assessment of druglikeness of ligands
The molecular docking process brought us to sort the 25 compounds among the 1565 drugs based on their highest negative binding energies. Further, among those 25 compounds, we selected 4 drugs which were involved in the process of cellular autophagy, based on the literature survey. Then after, all the 4 ligands were tested for their ability to follow Lipinski's rule of five and subjected to prediction of their ADME properties using an online analysis server SwissADME. The molecule which follows the four out of five parameters of Lipinski's rule, i.e. molecular weight less than or equal to 500 g/mol, MLogP less than or equal to 4.15, number of H-bond acceptors (N or O) less than or equal to 10 and number of H-bond donors (NH or OH) less than or equal to 5, is considered as an orally consumable drug and is potential enough to be used in the biological system. However, there are many drugs reported so far, which have immense pharmacological potentials but are less druggable and not likely to pass through the Lipinski's filter, have been approved by FDA as the potential drugs to increase the chemical discovery space and subsequently to be used for the clinical purposes (Degoey et al., 2018). We found the three out of four drugs: Nilotinib, Ponatinib, and Simeprevir, following the parameters for Lipinski's rule of five. All the above-mentioned drugs represented values indicating their ability to be used as a potential drug in the clinical systems. However, the drug Digoxin did not show its ability to follow the Lipinski's rule of five, as it showed violations with respect to molecular weight, number of H-bond donors and number of H-bond acceptors. Moreover, all the four ligands showed the value of MLogP less than 4.15, which indicates that the molecules are more likely to be in the hydrophilic environment and are favourable for their drulikeness. Overall, the drugs used in this study satisfy the filtering criteria based on Lipinski's rule of five and follow the ADME properties for being a potent orally active drug. The detailed comparable analysis of all other parameters of Lipinski's filter and ADME properties for all the four drugs are given in Table 2 and  Table 3, respectively.

MD simulation analysis
After molecular docking and assessment of the drugs to be involved in the process of cellular autophagy and their druglikeness, the best 4 hits were subjected to MD simulation to check the binding stability for a certain period of time. After the simulation, the obtained trajectories for RMSD, RMSF, Rg, and intermolecular hydrogen bonds were analysed to have an understanding about the stability of the protein-ligand complexes for 100 ns course of duration. The generated trajectory was used to calculate the RMSD after choosing the 'backbone' for least-squares fitting with aid of the 'g_rms' module of the GROMACS. The RMSD reveals the dynamism in the ligand's position with respect to the protein which indicates the stability of the poses of the docked complexes during the course of entire simulation period. From the RMSD graph ( Figure 3A), it is evident that the drug Ponatinib is most stable among all the 4 hit drugs with maximum deviation at 0.7 nm and showed a stable pattern after 50 ns. Followed by Ponatinib, the drug Simeprevir showed RMSD at $0.75 nm, whereas, the drugs Nilotinib and Digoxin attained the maximum RMSD value of $0.8 nm and $0.85 nm respectively and all the drugs exhibited a steady pattern after 50 ns. The drugs Ponatinib and Simeprevir showed lower RMSD value than that of the peptide alone ($0.85 nm) from the time scale of 80 ns ( Figure 3B), however, the RMSD values of drugs Nilotinib and Digoxin exhibited a 'converged'  Simeprevir interacts with the BH3 domain of Beclin1 through various covalent as well as non-covalent interactions. The bond between the ligand atoms is shown with maroon lines, whereas, the comparatively larger spheres represent the atoms of the amino acid of the peptide (BH3 domain of Beclin1). The dotted green line with the shown distance depicts the hydrogen bond between the atoms of the ligand and the peptide. The blue radiating arcs represent the amino acid residue of the peptide involved in hydrophobic interaction(s) (for details of various interactions in the docked complexes, refer to Table 1). pattern with that of the peptide alone. Hence, it is evident enough to suggest that all the drugs show a stable interaction with the peptide, however, the drug Ponatinib stands the best among all.
To further analyse the dynamic pattern of residues of the peptide with respect to its interaction with the four hit drugs, the RMSF for time scale of 100 ns was calculated with the help of the 'g_rmsf' module of GROMACS. The RMSF values for all the residues showed a range of fluctuation from $0.1-0.9 nm in the 100 ns period of simulation. From the RMSF graph ( Figure 3C), it can be seen that the N-terminal (residue 110-114) of the peptide complexed with the drug Ponatinib showed lower RMSF than that of the peptide alone. Moreover, the complex Simeprevir exhibited lower RMSF from the residue 119-126 as compared to that of the peptide alone, whereas, the complex Nilotinib showed lower RMSF than that of the RMSF of the peptide from the residue 110-123. Comparatively, slightly minimal fluctuations were also observed from residues 111-116 and from 123-124 which is an indicative of the fact that amino acids of these regions may contribute to hydrophobic interaction during the process of molecular docking ( Table 1). All of the residues of the complex Digoxin showed a higher RMSF than that of the peptide alone, suggesting that the hydrogen bonds formed during the molecular docking were not stable enough to persist for entire course of simulation. To verify this, a comparative analysis of formation of number of hydrogen bonds between the peptide and the drugs during the entire course of simulation for 100 ns was performed.
Hydrogen bonds are considered to play a crucial role in imparting the overall stability to the structure of the protein.
We assessed the intermolecular hydrogen bond formation for all the four complexes during the entire course of simulation. From the Figure 4, it is evident that in case of the drug Ponatinib complexed with the peptide (Figure 4A), the average number of hydrogen bond formed was four, whereas, the number reached up to a maximum of five. Moreover, in case of Simeprevir ( Figure 4B), on an average, two hydrogen bonds could be seen to be formed, however, exhibited a maximum of three hydrogen bonds. Another drug Nilotinib ( Figure 4C) formed two hydrogen bonds with the peptide, whereas, the drug Digoxin ( Figure 4D) formed an average of two hydrogen bonds while attained a maximum of three during the 100 ns period of simulation. Interestingly, it was observed that, in the peptide, Arg114, Arg115, Lys117, and Gln129 were the key residues which were imparting the Comparative analysis of all the four compounds with respect to parameters of predicting the ADME properties. TPSA represents topological polar surface area, LogS indicates the water solubility, PAINS is the acronym used for Pans assay interference structures and BS is the bioavailability score. Comparative representation of the Rg analysis for all the drug complexes as well as for the peptide alone. The analyses for all the essential parameters were done using the GROMACS simulation package. maximum contribution to the hydrogen bond formation with the drugs, suggesting that these residues have an important role in maintaining the stability of the protein-ligand interaction. Hence, by comparing both the RMSF and hydrogen bond graphs, we can suggest that the drugs Ponatinib, Simeprevir, and Nilotinib show lower RMSF values than that of the peptide alone which is an indication of exhibiting high degree of conformational stability owing to higher number of hydrogen bond formation, unlike in the case of the drug Digoxin.
Furthermore, we also analysed the level of compaction through calculating the radius of gyration for the peptide as well as its bound form with all of the four complexes with the help of using the 'g_gyrate' module of GROMACS. By definition, Rg accounts for the distance of all the atoms of a molecule of mass 'm' from the centre of mass of the molecule. From the graph, it can be visualized that initially up to 50 ns, the value for Rg varied for the peptide as well as for its bound form with all of the four complexes. However, as the simulation approached beyond 50 ns, the value for Rg was found to be similar and ranged from $0.75-0.85 nm for the peptide as well as for the docked complexes. Therefore, by analysing the graph ( Figure 3D), it can be suggested that the level of compaction is similar in the peptide as well as its drug-bound forms after the 50 ns course of simulation. Further, to analyse the variations obtained after the simulation in the overall binding energy of the peptide-drug complexes, the MM-PBSA method was employed.

Analyses of results of MM/PBSA method after MD simulation
In effort to compute the details of energetic investigations of all the protein-ligand complexes post-simulation, a calculation for the average binding energy was performed using the 'g_mmpbsa' tool. Free energy, being a state function, depends on the initial energy in the solution prior to binding of the ligand to the peptide and final energy in the solution post ligand binding to the peptide. Moreover, the residues of the peptide contributing to the average binding energy were also analysed using the same module. After performing the MM/PBSA method, a comparative analyses of the energy terms E MM , G polar , and G nonpolar of all the peptide-drug complexes were performed with the help of 201 snapshots stored from the production trajectories at time scale of every 0.1 ns from 80-100 ns course of simulation (Table 4). From the list of the peptide-bound drugs given in the Table 4, it is evident that the drug Nilotinib upon binding to the peptide, gave maximum binding energy of À15.468 kcal/mol, followed by the drug Simeprevir (-13.230 kcal/mol), Ponatinib (-11.001 kcal/mol), and Digoxin (-3.457 kcal/mol). However, it has been reported that the MM/PBSA method can also be used to re-score the docking score of the docked complexes so that there can be a clearer differentiation between an active or an inactive drug molecule (Thompson et al., 2008). Other than binding energy component of all the peptidedrug complexes, Van der Waals energy, electrostatic energy, polar solvation, and SASA energy are also mentioned in the Table 4 for the detailed investigation into the different energetic estimators of the protein-ligand complexes.
Furthermore, we analysed the average contribution of peptide's residues to the binding energy to have a clearer understanding about the binding interaction of the peptide and the drug. For this analysis, we executed the binding energy decomposition command for all the four peptidedrug complexes using 'g_mmpbsa' tool. As it can be seen from the Figure 5, the major energy contributor residues of the complex Ponatinib ( Figure 5A) were Glu110, Asp121, Asp124, and Gln129. The complex Simeprevir ( Figure 5B) was found with residues Val118, Leu122, Ile125, and Gln129 as the average participator residues to the binding energy. Moreover, the residues Arg115, Val118, Leu122, and Phe123 of the complex Nilotinib ( Figure 5C) were found as the maximum energy contributor residues, whereas, the complex Digoxin ( Figure 5D) showed the residues Arg114, Arg115, and Val118 as major contributors to binding energy. This observation suggests that mainly the residues Arg114, Arg115, and Gln129 show maximum contribution to the binding energy in case of all the four peptide-drug complexes. This observation is in agreement with the residues observed in imparting contribution to the post-simulation hydrogen bond formation and thereby, in the stability of the docked complexes.
The current study allows us to acquire a deeper analysis by investigating the various interactions between the drugs viz. Ponatinib, Nilotinib, Simeprevir, and Digoxin, and the BH3 domain of the protein Beclin1. This interaction study suggests Note. The table gives a comparative analysis for all the energetic components formed between the peptide-drug complexes after simulation. The different forms of energy were calculated by MM/PBSA method for all the drugs bound with the peptide. Figure 5. Contribution of individual residues to the binding energy. The plots depict the contributor residues of each of the four peptide-drug complexes which impart the maximum to the contribution energy. The contribution of the residues to the overall binding energy was calculated through MM/PBSA method. The residues contributing maximally are labelled in each plot.
the importance of a few key residues of the peptide which are consistent throughout the process of molecular docking, MD simulation, and MM/PBSA which lead to strengthening of the favourable binding interaction with the drugs used in this study.

Conclusion
We have reported structural aspects of various proteins related to diseases previously (Dubey et al., 2005;Singh et al., 2008). The objective of the present study focuses upon the structure based drug design for inhibition of Beclin1-Bcl-2 interaction thereby inhibiting the process of cellular autophagy, which will lead to activation of apoptosis and culminate to tumor clearance. A large number of FDAapproved drugs were selected for the current study and finally, four of them were chosen with the highest binding affinity which interact with the BH3 domain of Beclin1; the actual site of interaction of Bcl-2. Furthermore, considering only the essential parameter viz. RMSD, RMSF, Rg, and formation of intermolecular hydrogen bond generated after the MD simulation, the drug Ponatinib among all the four drugs showed better performance in the peptide-complex ranking. However, further analyses in terms of post-simulation binding energy, the drug Nilotinib exhibited the highest negative binding energy and the drug Ponatinib showed at third position in the peptide-complex ranking. Hence, overall, based on different investigating parameters, we can strongly suggest that the top three drugs viz. Ponatinib, Simeprevir, and Nilotinib can be the potent regimen candidates, however, the drug Ponatinib can be used as a remarkable drug for further in vitro and in vivo research in the domain of autophagy-apoptosis crosstalk. Moreover, the mechanism of inhibition of autophagy and activation of apoptosis by using these drugs is still elusive, however, together with our computational study on targeting the protein-protein interaction of Beclin1 and Bcl-2 for inhibiting the autophagy using the FDA-approved drugs and thereby activating apoptosis, may become an important consideration for performing further in vitro and in vivo investigations and ultimately provide a novel strategy which may give an insight into the protein dynamics and play a crucial role in the drug discovery.