Binding site hotspot map of PI3Kα and mTOR in the presence of selective and dual ATP-competitive inhibitors

Abstract The PI3K/Akt/mTOR signaling pathway plays a pivotal role in cellular metabolism, growth and survival. PI3Kα hyperactivation impairs downstream signaling, including mTOR regulation, and are linked to poor prognosis and refractory cancer treatment. To support multi-target drug discovery, we took advantage from existing PI3Kα and mTOR crystallographic structures to map similarities and differences in their ATP-binding pockets in the presence of selective or dual inhibitors. Molecular dynamics and MM/PBSA calculations were employed to study the binding profile and identify the relative contribution of binding site residues. Our analysis showed that while varying parameters of solute and solvent dielectric constant interfered in the absolute binding free energy, it had no effect in the relative per residue contribution. In all complexes, the most important interactions were observed within 3–3.5 Å from inhibitors, responding for ∼75–100% of the total calculated interaction energy. While closest residues are essential for the strength of the binding of all ligands, more distant residues seem to have a larger impact on the binding of the dual inhibitor, as observed for PI3Kα residues Phe934, Lys802 and Asp805 and, mTOR residues Leu2192, Phe2358, Leu2354, Lys2187 and Tyr2225. A detailed description of individual residue contribution in the presence of selective or dual inhibitors is provided as an effort to improve the understanding of molecular mechanisms controlling multi-target inhibition. This work provides key information to support further studies seeking the rational design of potent PI3K/mTOR dual inhibitors for cancer treatment. Communicated by Ramaswamy H. Sarma

Despite pharmacological efforts to disrupt the downstream cell survival signaling through single-target therapy using PI3K or mTOR inhibitors, the complex biochemical network regulation in cancer pathology has made such task extremely challenging (Hillmann & Fabbro, 2019;Shukuya et al., 2019). The use of pan-class I PI3K inhibitors have elicited high toxicity, and compounds selectively targeting PI3K isoforms have also failed due to the demand for high dosage during treatment (Cheson et al., 2019;Patsouris et al., 2019). In addition, the inhibition of PI3K/AKT is not enough to completely inactivate the mTOR downstream pathway as mTORC1 can be independently activate by the signalling pathways TSC1/TSC2/Rheb, LKB1/AMPK and Rag GTPases (Chi, 2012;Dowling et al., 2010). In this scenario, similarities between the ATP-binding site of PI3Ks and mTOR (shared homology $32%) represents a unique opportunity to explore multi-target pharmacological approaches and some studies have focused on the multikinase inhibitory therapy to tackle prostate cancer Luszczak et al., 2020), lung cancer (Wu et al., 2019), lymphoma (Bresin et al., 2020;Tarantelli et al., 2014Tarantelli et al., , 2020, gastric cancer , ovarian cancer , breast cancer (Kennedy et al., 2020), neuroblastoma (Mohlin et al., 2019;O'Neil et al., 2014), haematological cancer (Reidy et al., 2014) and rare genetic diseases (Hillmann & Fabbro, 2019), among other malignancies.
In this work, to shed light in the binding mechanism controlling dual-PI3K/mTOR inhibition and investigate the role of binding-site residues interacting with selective or dual inhibitors, we took advantage of crystallographic coordinates of PI103 in complex with PI3Ka (PDB ID: 4L23) (Zhao et al., 2014) or mTOR (PDB ID: 4JT6) (Yang et al., 2013), alpelisib complexed with PI3Ka (PDB ID: 4JPS) (Furet et al., 2013) and torin-2 complexed with mTOR (PDB ID: 4JSX) (Yang et al., 2013). Alpelisib is a high selective PI3Ka inhibitor approved for the treatment of metastatic breast cancer in 2019 (Narayan et al., 2021), torin-2 is a potent mTOR selective inhibitor  and PI103 is a potent dualinhibitor targeting class I PI3Ks and mTOR proteins (Park et al., 2008). Molecular dynamics simulations were employed for structural relaxation and sampling while the molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) approach was employed to calculate the individual amino acid binding interactions, providing a description of the dependence of the total binding energy with the binding site radius. Our results characterize, under an energy-based point of view, key residues responding for the selective inhibition of PI3Ka and mTOR and highlight those participating in the mechanism of dual inhibition. Results from this work have potential to guide the design of new therapeutics to be used in the target-therapy against cancer and other diseases.

Molecular dynamics
To account for protein flexibility, we performed molecular dynamics simulations, for each system, using GROMACS 2019 (Abraham et al., 2015) package with AmberSB99 force field Figure 1. Schematic representation of PI3K/Akt/mTOR signalling pathway. mTOR is present in two distinct complexes, mTORC1 and mTORC2, which are also subject to activation by alternative pathways. (Hornak et al., 2006). Simulation consisted in $190,000 atoms for PI3Ka-inhibitor systems, and $292,000 atoms for mTORinhibitor systems. In all systems tip3p water model (Jorgensen et al., 1983) was employed to describe water molecules and net charge was neutralized by adding Na þ and Cl À ions at 0.15 M concentration. Total energy minimization was accomplished by combining steepest-descent and conjugate gradient algorithms, in sequence. At total, each system was equilibrated using integration steps of 1 fs through 2 ns following gradually decreasing ligand restraint force constants along 4 steps of 250 ps each. Long-range electrostatics was modelled with the Particle Mesh Ewald (PME) method (Essmann et al., 1995), temperature coupling was done with Nose-Hoover thermostat at 310.15 K and the Parrinello-Rahman barostat (Parrinello & Rahman, 1981) with a reference pressure of 1 bar. LINCS algorithm was used to constrain covalent bonds to their equilibrium length (Hess et al., 1997). The integration steps of all simulations during the production step were set to 2 fs. At the end, a trajectory of 500 ns was obtained for each system and the distance root-mean-square deviation (rmsd) and the root-mean-square-fluctuation (rmsf) were calculated to access structural stability of the complex.

Relative binding free energy calculations
Molecular Mechanics Poisson-Boltzmann Surface Area (MM/ PBSA) calculations were performed through the g_mmpbsa (Kumari et al., 2014) code. Relative binding energies were calculated from 100 snapshots representing the last 10 ns of each 100-ns trajectory, using the MM/PBSA approach (Eq. (1)). Additional comparison between calculated values using 100 and 1000 snapshots were performed by sampling the last 10 ns of 500-ns trajectories.
where DG bind (the total binding free energy) is the free energy difference between the bound complex (formed by the association between the protein and the ligand), and the isolated protein and ligand in solvent. The free energy (G) for each individual entity is estimated as in Eq. (2): where E MM represents the sum of the Molecular Mechanics internal energy of the molecule, G solvation is the free energy of solvation, T is the temperature in units of Kelvin and S represents the conformational entropy of the system. The Molecular Mechanics internal energy is given by the bounded and the nonbounded terms: where, for each individual molecular entity, E bounded represents bonded interactions consisting of bond-stretch, anglebend, torsion and improper-dihedral interactions, while E nonbounded represents the sum of both, the electrostatics and the van der Waals interactions calculated through Coulomb and Lennard-Jones potential functions, respectively. G solvation represents the amount of energy spent to transfer a solute from vacuum into solvent and is calculated using an implicit solvent model: where G polar is the polar solvation energy of a molecule, estimated by solving the Poisson-Boltzmann (PB) equation and G nonpolar is the term for the nonpolar solvation energy which is approximated by a solvent-accessible surface area (SASA) (Eisenhaber et al., 1995) term, based on the assumption it has a linear dependency on the G nonpolar term: where c is a coefficient related to surface tension dependent of the solvent employed, A represents SASA and b is a fitting parameter.
In this work, we used a single-trajectory protocol where only the protein-ligand complex is simulated and the conformations of protein and ligand in the bound and unbound states are assumed to be identical, and therefore their E bounded and E bounded terms cancels out. In addition, to assess the dependence of DG bind on the solute dielectric constant e solute value, the solute dielectric constant was set to 2 or 8, in combination with three distinct values (20, 40 and 80) for the solvent dielectric constants, to account for the effect of solvent in the binding cleft. Other parameters in g_mmpbsa were left as default (Kumari et al., 2014).

Interaction entropy (IE)
Interaction entropy (IE) has emerged as a fast and accurate method to estimate the entropic contribution directly from MD simulations, at negligible computational cost (Duan et al., 2016;. The entropic contribution can be calculated by the IE method using the following term: where DE int pl represents the fluctuation of protein-ligand interaction energy E int pl around the average interaction energy E int pl D E , as depicted in Eq. (7): The term E int pl D E is obtained as shown below, While ⅇ bDE int pl is calculated as follows: 2.6. Total interaction energy and relative amino acid contributions One of the great advantages in the MM/PBSA method is its ability to perform reliable comparison among distinct ligands interacting with the similar target. In this work, python scripts 'MmPbSaStat.py' and 'MmPbSaDecomp.py' (Kumari et al., 2014) were employed to perform per-residue energy decomposition based on MM/PBSA calculated energy components E MM , G polar , and G nonpolar were performed in order to identify key amino acid contributions to the binding of distinct ligands in the binding site of PI3Ka and mTOR. Total interaction energy is represented as DG MMPBSA , or DG MMPBSA-IE when the entropic contribution, calculated through the IE method, is considered. Following, to observe the fluctuation of the total interaction energy, the individual interaction energy of every amino acid residue within 10 Å from the ligand was summed up and plotted against its radial distance to the ligand. Distances were calculated using the last snapshot of 100 ns-MD trajectories. To observe the relative contribution of each amino acid residue to the total binding energy, hot spots which normally respond for large changes in this curve are plotted in horizontal bar charts representing (i) the interaction energy (in kJ/mol) of ligands with individual amino acids residues, from which it is possible to visually assess its attractive or repulsive nature in the given system; and, (ii) the most important residues contributing to the stability of the complex in the column of residues at the left side. In addition, to better compare the action of conserved and non-conserved residues across PI3Ka and mTOR binding sites, calculated energies were depicted through a heatmap plot representation, where similarities and differences are easily assessed.

Results and discussions
Initially, to better analyze similarities and differences between PI3Ka and mTOR ATP-binding site residues, crystallographic structures 4L23 and 4JT6 were superimposed to identify the spatial arrangement of conserved and non-conserved residues located within 6 Å from inhibitors. At total, it was identified 10 conserved and 07 non-conserved residues. A visual representation of the superimposed ATP-binding site is shown in Figure 2, and a list of conserved and non-conserved PI3Ka-mTOR pair of residues, numbered from 1 to 17, is shown in Table 1, supplementary material. Following, each protein-inhibitor complex was protonated at physiological pH and submitted to molecular dynamics simulations and relative binding energy calculations through MM/PBSA approach as described below.

Molecular dynamics simulation and structural stability
Due to the crystallographic nature of the selected protein-ligand complexes, it is of paramount importance to allow them to reach conformational equilibrium under explicit solvent at physiological pH prior to further energy analysis. Such approach is necessary because during crystal packing, sometimes spurious contacts are formed and both ligand and side chains are trapped in local/global minima which do not represent the equilibrium conformation at physiological pH. Therefore, to ensure a correct analysis of relative amino acid contribution, each selected protein-ligand complex must be adequately protonated and solvated, and any molecular clash must be removed, improving, and stabilizing their internal atomic contact network. In this study, molecular dynamics (MD) using the AmberSB99 force field and TIP3P water model were employed to sample protein-ligand conformations used during MM/PBSA calculations. As observed in Figure 4(A) and 4(B), which represents the rmsd fluctuation for the entire structure, all complexes were stable at 100 ns. Among them, the complex mTOR-PI103 took more time to stabilize, reaching plateau at about 82 ns of simulation. For the sake of comparison, MM/PBSA were calculated using snapshots from 10 ns of MD simulations extracted from the time 90 ns and 490 ns (yellow shadow). In addition, it was verified the root-mean square fluctuation (rmsf) of ATPbinding pocket residues in close contact with ligands along the first 100 ns and compared with fluctuation along the whole simulation (500 ns). As observed in Figure 4(C)-(F) only small displacements of binding site residues were observed, suggesting that most of the rearrangement observed in Figure 4(A) and 4(B) are related to domains external to the ATP-binding pocket.

Calculated interaction energy
In the past decades, the MM/PBSA approach has been successfully used to calculate the interaction energy between a wide range of small ligands and biological molecules (Poli et al., 2020). Different from single point energy calculations, which is commonly applied to crystallographic structures, MM/PBSA requires a set of snapshots from molecular dynamics simulations to accurately access and average the effect of small fluctuations in intermolecular interactions. Nevertheless, while the accuracy of the calculated absolute binding free energy in MM/ PBSA approach is disputable and very system-dependent, the analysis of per residue interaction energy among similar biological targets bound to distinct ligands has been of great value to rank individual amino acid contributions, the effect of mutations and to guide rational drug design. In this work, to ensure the consistence of MM/PBSA results during our analysis, we sampled conformations for MM/PBSA calculations from two distinct points along the simulation time: the first, immediately after the stabilization of all structures (from 90 to 100 ns), and the second at the end of the simulation (from 490 to 500 ns). Both ranges are depicted in yellow shadow in Figure 4(A) and 4(B). As observed in  (2021), who observed an increase in the entropy with number of energies considered. Following, we tested distinct  Table 1, supplementary material. combinations of protein and solvent dielectric constant parameters under SAV or SASA nonpolar models during the free energy of solvation calculations. As shown in Table 4, supplementary material, changes in the dielectric constant had only small impact in the total free binding energy and showed consistence under distinct nonpolar models. On the other hand, the use of nonpolar model SAV produced results 2-times the magnitude of SASA model. Such overestimation in SAV energy were expected, as previously observed by Kumari and colleagues (2014), who compared calculated MM/PBSA binding free energy with a set of experimental data. Therefore, in this work we reported MM/PBSA and IE (MM/PBSA-IE) energies  calculated from 90 to 100 ns, using 100 snapshots, and applying SASA model to describe the per residues interactions in each complex. In addition, considering that the ATP-binding site in these systems consist in inner cavities, with low charge and limited exposition to the external medium (Li et al., 2013), the parameters of protein and solvent dielectric constant were set as 2 and 20, respectively. By summing up the interaction energy of each residue within 10.0 Å from the ligand, it was possible to analyse the fluctuation in the total interaction energy as a function of the radial distance. As shown in Figures 5 and 6 (for PI3Ka and mTOR complexes, respectively), despite the changes observed in the calculated total binding energy under distinct values of the solvent dielectric constant (20, 40 or 80), for all complexes, the fluctuations on the binding energy along the radius distance kept the same behaviour, indicating that the residue relevance do not change upon distinct parameters.
MM/PBSA calculations with the solute (solvent) dielectric constant set to 2 (20) showed that PI103 interacts with PI3Ka through À44.10 kJ/mol, with residues located withing 3.5 Å responding for about 75% (À32.92 kJ/mol). The calculated total energy for PI3Ka-alpelisib complex was À52.13 KJ/mol, with residues within 3 Å responding for an interaction of À53.46 kJ/mol. After these distances, in both structures, it was also observed oscillations, with values returning to the same level after 7.5 Å.
In the same fashion, it was calculated total interaction energy as a function of the radius distance for mTOR complexes, shows that mTOR interacts with PI103 with a strength of À55.87 kJ/mol, with residues located withing 3.5 Å responding for about 97% (À54.32 kJ/mol). The calculated total energy for mTOR-torin2 complex was À38.12 kJ/mol, with residues within 3.0 Å responding for an interaction of À42.35 kJ/mol. After these distances, in both structures, it was also observed oscillations, with values returning to the same level after 7.0 Å.
It worth to note that, in both PI3Ka complexes, the binding of inhibitors was stabilized mainly by hydrophobic residues, while amino acids with charged side chains are responding for repulsive effects. Nevertheless, despite remarkable behavioral similarities (attractive or repulsive effect) among most of the residues in PI3Ka complexes, significant differences were observed when Ile932 (DE ¼ 4.84 kJ/mol, Lys802 (DE ¼ 6.63 kJ/ mol), Asp810 (DE ¼ 3.26 kJ/mol), Val851 (DE ¼ 2.75 kJ/mol), Arg852 (DE ¼ 2.73 kJ/mol) and Tyr836 (DE ¼ 2.19 kJ/mol) and interacted with one or another inhibitor. Indeed, as shown in Figure 8, the non-polar amino acid Ile932 is close to the pyridofuro-pyrimidine group of PI103 ( Figure 3B, region ii), while in the PI3Ka-alpelisib complex it is close to the methyl-trifluoromethylpropanyl-pyridinyl group ( Figure 3A, region i). In the PI3Ka-PI103 complex, residue Val851 makes a hydrogen bond with the to the oxygen atom (O2) in the PI103 morpholine ring ( Figure 3B, region iii), while in the PI3Ka-alpelisib complex, Val851 forms a nitrogen bond with atom (N1) in the thiazole ring ( Figure 3A, region ii) and are close to the nitrogen atom (N2) in the pyrrolidine-dicarboxamide group ( Figure 3A, region iii). Interestingly, such interaction is conserved in mTOR-PI103 complex, with the oxygen atom (O2) in the PI103 morpholine forming a hydrogen bond with the mTOR residue Val2240. In addition, as shown in Figure 10, Va2240 also forms a hydrogen bond to the nitrogen atom (N2) in the benzo-naphthyridine ring of torin-2 ( Figure 3C, region ii).
Interestingly, it is believed that selective inhibition of PI3K relies on a small number of non-conserved residues located in the periphery of the ATP-binding site, being four of them ( 852 RNSH 855 , in PI3Ka) located near the hinge (called region 1), and few others located at the P-loop (called region 2) (Frazzetto et al., 2008). Our results show that non-conserved residues in the ringe region hold weak interactions but distinct behaviours when interacting with PI103 or alpelisib. For instance, while the attractive strength of residue His855 was larger for alpelisib (almost 8x) than for PI103, residues Arg852 and Ser854 showed attractive interaction with PI103, but repelled alpelisib. The relevance of His855 for PI3Ka selectivity was also shown previously by Zheng et al. (2011), using in vitro mutagenesis to evaluate the binding of the inhibitor PIK-75 (Zheng et al., 2011). In addition, residue Gln859, which is known as an important player in PI3K isoform selectivity (Furet et al., 2013;Zheng et al., 2012), showed attractive interaction for alpelisib, making two hydrogen bonds with the oxygen (O1) and nitrogen (N3) atoms, but acted repulsively in the presence of the dual-inhibitor PI103 (see Figures 7 and 8).
Trp2239 in mTOR binding pocket shields the effect of Lys2171. Indeed, the spatial orientation of Trp780 and Trp2239 seems to be pivotal for the binding of single and dual inhibitors as shown in Figures 8 and 10. In addition, our analysis shows the strong repulsive effect of negative charged residues Asp933 and Asp2357 and highlight the strongest attractive interaction of hydrophobic (PI3Ka-mTOR) residues Ile932-Ile2356, followed by Ile848-Ile2237, Met922-Met2345 and Val851-Val2240. Among non-conserved hydrophobic residues, the strongest attractive interactions were observed with Ile800-Leu2185 and Val850-Trp2239.
Interestingly, when interacting with the dual inhibitor PI103, the conserved and positively charged (PI3Ka-mTOR) residues Asp810-Asp2195 and Lys802-Lys2187 showed opposite behaviour. For Asp810 and Asp2195, the distances to the phenol group of PI103 (region i, Figure 2B) impacts the magnitude of the calculated interaction energy, with Asp810 (in PI3Ka) strongly repelling PI103, whilst Asp2195 (in mTOR) holds a near neutral to weak attractive interaction (DE ¼ 10.36 kJ/mol). Similarly, the conserved PI3Ka residue Lys802 showed stronger repulsion to PI103 than did Lys2187, in mTOR complex (DE ¼ 7.46 kJ/mol). For non-conserved residues Asp805-Glu2190, the spatial arrangement of PI103 into the binding pocket results in distinct interaction pattern (DE ¼ 10.88 kJ/mol), while for Trp780-Lys2171 (DE ¼ 2.55 kJ/mol), the presence of residue Trp2239 in mTOR binding pocket shields the interaction between Lys2171 and PI103. In addition, there is a clear opportunity for improving the total binding energy of new inhibitors by increasing their interaction with PI3Ka residues Asp933 and Leu814 and, mTOR residues Asp2357 and Met2199. To better visualize the attractive or repulsive behaviour of binding site residues, the calculated per residue interaction energies were plotted onto the cartoon representation of spatial coordinates for each complex, as shown in Supplementary Figure 1, supplementary material.

Conclusions
In this work, MM/PBSA approach was employed to calculate the relative individual amino acid contribution to the binding of selective or dual inhibitor into PI3Ka and mTOR ATP-binding site. Calculations identified hydrophobic PI3Ka residues Ile932, Ile848, Ile800, Val850, Val851 and Met922 as the most important for the binding of PI103 and alpelisib. In addition, calculations also identify hydrophobic mTOR residues Figure 11. Heatmap showing similarities and differences between PI3Ka and mTOR residues due to the ligand effect. Left column shows PI3Ka residues, whilst mTOR residues are shown in the right. As shown in the scale, attractive energies are shown in red, while neutral and repulsive interactions are shown in green and blue, respectively. Trp2239, Ile2356, Leu2185, Ile2237 and Met2345 as the most important for the binding of torin-2 and PI103 inhibitors. Interestingly, PI3Ka residues Asp810 and Lys802 showed opposite behaviour when compared with their conserved mTOR counterparts Asp2195 and Lys2187. While Asp810 repelled alpelisib and PI103, Asp2195 showed neutral behaviour with torin-2 and PI103. On the other hand, Lys802 repelled alpelisib and PI103, while Lys2187 showed neutral behaviour with PI103 and repelled torin2. In addition, while Asp805 showed neutral behaviour in the presence of alpelisib but acted attractively in the presence of PI103, its counterpart mTOR residue Glu2190 repelled both, torin2 and PI103. These results suggest that improving/balancing molecular contacts with these residues could be beneficial to tune selectivity in drug design studies. When our attention was turned to the binding of the dual-inhibitor PI103, we identify key interactions in PI3Ka residues Phe934, Lys802, Asp805 and mTOR residues Leu2192, Phe2358, Leu2354, Lys2187 and Tyr2225 as their behaviour differs to that observed with respective selective inhibitors. Among interactions that could be improved during dual inhibition are those with PI3Ka residues Asp810, Leu807, Leu814 and Lys802 and, with mTOR residues Asp2195, Met2199, Lys2187, Glu2190 and Lys2171. Therefore, the nature of each of those interaction should be carefully analyzed during the design of new multi-target inhibitors as they represent a clear opportunity for enhancing the binding strength in both proteins. Altogether, our results depict a clear panel representing the interaction network of selective and dual ATP-binding site inhibitors and highlight the key similarities and differences among their binding mechanism. Therefore, this work provides an improved understanding in the complex behaviour of PI3Ka and mTOR residues controlling the binding of selective and dual ATP-binding site inhibitors and pave the way for the design of new multi-target single agents with better affinity and reduced side effects.

Disclosure statement
The authors declare no conflicts of interest.  Authors' contribution FFNA performed the calculations and analysis. FJFS and FLSO prepared figures and tables. GZ supervised the work and performed analysis. GZ and JVC wrote the initial manuscript. All authors discussed the results, modified the manuscript, and approved the final version.