Deciphering the resistance mechanism of RET kinase mutant against vandetanib and nintedanib using molecular dynamics simulations

Abstract The RET protein is a transmembrane receptor tyrosine kinase (RTK) whose oncogenic mutations or fusions are closely related to human cancers such as thyroid and non-small cell lung cancer. Vandetanib as a clinical-approved protein-tyrosine kinase inhibitor (TKI) exhibits anti-cancer efficacy by blocking the RET ATP-binding site, but drug resistance was observed for the RETG810A mutant. Recent studies have identified another TKI nintedanib as an effective molecule to inhibit vandetanib-resistant RETG810A. However, there is no clear evidence of why nintedanib and vandetanib displayed different inhibitory activities towards RETG810A. Here, we exploited molecular dynamic (MD) simulations to compare the interactions of the RETG810A mutant with nintedanib and vandetanib. A higher structural flexibility of the activation loop was observed in the nintedanib-bound RETG810A, which may result in discrepant autophosphorylation activity in the nintedanib- and vandetanib-bound RET kinase, causing differentiated pharmacological effects of the two compounds. Molecular mechanics/Poisson-Bolzmann surface area method suggested that nintedanib had a higher affinity towards RETG810A over vandetanib, accounting for its better inhibitory effect as an ATP-competitive compound. These results depicted the underlying mechanism for the different inhibitory efficacy of nintedanib and vandetanib on RETG810A from both conformational and energetic aspects. Furthermore, we also found that both compounds maintained the ‘DFG-in, αC-helix-in, and activation loop-open’ conformation of RETG810A, which is the characteristic of the active state. Together, our results provide comprehensive mechanistic insights into nintedanib’s capability in inhibiting vandetanib-resistant RET mutant and enlighten future structural-based optimisation of RET TKIs to overcome drug resistance.


Introduction
The RET protein, encoded by the gene RET (abbreviation for the 'REarranged during Transfection'), is a transmembrane receptor belonging to the receptor tyrosine kinases (RTKs) family [1,2]. As an important promoter for the development of human kidneys and enteric nervous system [3][4][5], RET receives the extracellular signals such as the glial cell line-derived neurotropic factor (GDNF) with co-receptors [6], activates its kinase activity through homodimerization and autophosphorylation [7], and ultimately facilitates the down-stream signalling cascades including MAPK, PI3K, and PKA pathways to activate cell proliferation [8][9][10][11]. Aberrant activation of RET resulting from either oncogenic mutations of RET gene or chromosomal rearrangements can lead to the occurrence of solid tumours such as medullar thyroid cancer (MTC) and non-small cell lung cancer (NSCLC) [12][13][14].
RET structure consists of an extracellular domain (residues 1-635) containing four cadherin-like domains and a cysteine-rich domain, a transmembrane domain, a juxtamembrane region, and the kinase domain (residues 724-1016) [15]. The intracellular kinase domain, which shows high structural similarity with other RTKs, sets up a suitable target for the inhibition of RET activities using small-molecule inhibitors in the treatment of RET-associated diseases [16,17]. Several multi-kinase TKIs that possess broad-spectrum RTK inhibitory effects have shown activities against RET, including cabozantinib, lenvatinib, sorafenib, vandetanib, ponatinib, sunitinib, and alectinib [15,18,19]. Among them, vandetanib has received the approval from US Food and Drug Administration for the treatment of MTC in patients with oncogenic RET mutations.
However, the emergence of drug-resistant mutations represents an upcoming challenge, as a vast amount of cases in clinic and in laboratory have been reported [20]. For example, the RET G810A mutant that has an alteration of glycine residue to alanine at the ATP-binding pocket has recently identified as a vandetanib-resistant target [21]. Since mutations on RET kinase domain may result in the drug resistance to one or several specific RET TKIs, current solutions to overcome the challenge are by screening other potentially active compound among the clinically relevant RET TKIs. As such, another TKI nintedanib, which was previously designed for targeting angiokinases (e.g. VEGFR), was found to be effective in inhibiting RET G810A [21,22].
Until now, the molecular basis of how a single mutation in the RET kinase domain can result in drug resistance and why various TKIs have different inhibitory effects on the same drug-resistant mutant remains poorly understood [23]. Recent studies have resolved the co-crystal structures of wild-type RET kinase domain in complex with vandetanib [17] or nintedanib [20]. Nonetheless, there are no direct structural evidence of how nintedanib and vandetanib differently interact with the RET G810A mutant. In addition, whether the two compounds may cause distinct allosteric effects on the RET conformational dynamics is also unknown.
Here, using computational methods we established the RET G810A Ànintedanib and RET G810A Àvandetanib models and exploited molecular dynamic (MD) simulations to investigate discrepancies in the interaction modes between the two compounds towards the RET G810A mutant. We compared the binding affinities of the two compounds towards the ATP-binding pocket and the global conformational changes of RET G810A kinase domain, thus revealing the underlying mechanisms for the differentiated inhibitory effects of nintedanib and vandetanib. This study may promote future optimisation of RET TKIs.

System preparation
Based on the crystal structures, two systems were built for MD simulations, namely RET G810A Ànintedanib and RET G810A Àvandetanib, The crystal structures of wild-type RET kinase bound to nintedanib (PDB ID: 6NEC) [20] and vandetanib (PDB ID: 2IVU) [17] were used to define the initial coordinates of RET G810A mutant in complex with nintedanib or vandetanib. The RET G810 residue in both complexes were mutated to alanine using Discovery Studio 2019. The non-terminal missing residues were repaired using Modeller [24] implanted in the UCSF Chimera software.

All-atom MD simulations
Initial parameter files for minimisation and simulation were prepared first using Amber ff14SB force field [25] and general Amber force field (GAFF) [26]. Hydrogens of the RET proteins were added using the LEaP program in Amber18 package. The pharmaceutical compounds nintedanib and vandetanib were first protonated to electroneutral, and then their force field parameters were built by antechamber [26,27]. Both systems were solvated to truncated octahedron transferable intermolecular potential three point (TIP3P) water box [28]. Na þ and Clcounterions were added to neutralise system charges and then to the concentration of 0.154 molÁL À1 to simulate the physiological environment. Second, for each system, two rounds of energy minimizations were carried out for structural optimizations [29][30][31][32][33][34]. All atoms of RET protein and the inhibitor were restrained in the first round of minimisation while were without any constraints in the second round. Then, both systems were heated from 0 K to 300 K within 300 ps in a canonical ensemble (NVT), followed by 700 ps equilibration simulations. Finally, 300 ns production runs were performed under isothermal isobaric condition using the CUDA version of PMEMD. For each system, the production run was repeated in three replicates by setting random initial velocities, adding up a total of $1.8 ls simulation time for the two systems. During both simulations, the particle mesh Ewald method [35] was employed for long-range electrostatic interactions, and the cut-off distance for short-range electrostatic interactions as well as van der Waals interactions was set to 10 Å. Covalent bonds involving hydrogens were restricted using SHAKE methods.

Dynamic cross-correlation matrix (DCCM) analysis
The dynamic cross-correlation matrix (DCCM) of all protein Ca atoms were calculated to reflect the inter-residue correlations [36][37][38][39][40]. The cross-correlation coefficient C ij was calculated by: where r i and r j represents the positions of the ith and jth Ca atoms.

Hydrogen-bond analysis
Inter-and intra-molecule hydrogen-bonds were determined using simple geometric criteria. A hydrogen bond is determined as a hydrogen atom being in the between of two heavy atoms among N, O, and F, with the distance of heavy atoms less than 3.5 Å and the bond angle greater than 135 . A comprehensive scan of hydrogen-bonds throughout MD trajectories was realised by using the cpptraj software [41].

Determination of representative structures
The representative structures from certain groups of snapshots in the MD trajectories were extracted using the cluster analyses method [42][43][44]. Here, the hierarchical agglomerative (bottom-up) algorithm was adopted, with a distance cut-off value of 2.0 Å. Before cluster analysis, the snapshots were superimposed using all Ca atoms to eliminate overall rotations and transitions.

Molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) calculations
To evaluate the binding affinities between nintedanib/vandetanib and RET G810A , we used MM/PBSA plugin from MMPBSA.py in Amber18 package to calculated the binding free energies [45][46][47][48]. The Gibbs free energy change DG binding of the binding process is defined by Eq. (2): where the Gibbs free energy of each system is the algebraic sum of three compartments, including the molecular mechanical energy (E MM ), solvation energy (G solv ) and the entropic components (ÀTS): By defining molecular mechanical energy change DE MM , solvation energy change DG solv , and entropy change DS, Eq. (3) can be simplified to: Then, DE MM is ulteriorly decomposed into the van der Waals component DE vdW , the electrostatic component DE ele and the internal component DE int : Because ligands and RET proteins adopted same configurations in both unbound and bound states within a single snapshot, DE int is always equal to zero in our calculations.
For the calculation of DG solv in Eq. (4), we divided it into the polar part DE PB and the nonpolar part DE nonpolar and applied Poisson-Boltzmann continuum solvent model: The nonpolar component DE nonpolar was calculated based on the solvent-accessible surface-area (SASA) with the following format, where the surface tension parameter c was set to 0.00542 kcalÁmol À1 ÁÅ À2 , and the solvation parameter b was set to 0.92 kcal/mol.
The conformation entropy component (ÀTDS) was omitted here for the simplification of our calculations since no significant conformational shifts on the pharmaceutical compounds or RET backbones were observed during our MD simulations.
For each system, trajectory snapshots from the last 180 ns of the three independent replicates were selected and averaged for energy calculations. The dielectric constant was set to 1.0 and 80.0 for the solute and solvent, respectively. Furthermore, a per-residue scheme was applied to determine the energetical contributions of each residue of the binding process.

3.
1. An overview on the overall conformational dynamics of RET G810A Ànintedanib and RET G810A Àvandetanib To investigate the underlying mechanism how nintedanib can inhibit the vandetanibresistant RET G810A mutant, the structural basis between RET G810A Ànintedanib and RET G810A Àvandetanib was compared. To the best of our knowledge, there is no direct structural evidence of how nintedanib or vandetanib interact with the RET G810A mutant, and the crystal structures of RET kinase domain shared great similarities between RET WT Ànintedanib and RET WT Àvandetanib ( Figure 1A). Thus, we constructed the nintedanib-and vandetanib-bound RET G810A complexes using computational model based on X-ray crystal structures [17,20] and conducted MD simulations to compare the difference of the ligand-protein interactions.
The root-mean-square deviations (RMSDs) of RET G810A C a atoms relative to the corresponding initial crystal structures were calculated to reflect the overall conformational changes of RET during 300 ns simulations. After $100 ns simulations the RMSD values of the RET G810A Ànintedanib system and RET G810A Àvandetanib system were stabilised at 2.83 ± 0.46 Å and 2.78 ± 0.54 Å, respectively ( Figure 1B). This suggested that both systems reached equilibrium where no major conformational shifts occurred on the backbone of RET G810A . The reference structures of the two systems showed a low relative RMSD value of 0.36 Å, indicating that no substantial overall conformational distinctions were observed among the two systems during simulations.
Next, by calculating the root-mean-square fluctuations (RMSFs) of each RET residue using the coordinates of the C a atoms, the mobilities of corresponding residues during the equilibrium-state simulations were compared. As shown in Figure 1C, in most regions of RET G810A of the two systems, the backbone atoms showed low flexibilities with the RMSF values under 1.0 Å. The most significant distinctions between the two systems lie in the loop region that links the aD and aE helices (residues 821-846) and the activation loop (A loop) (residues 892-921) of RET G810A . Since the former loop region was missing in the initial crystal structures due to its hyperdynamic nature and may have little impact on the kinase activities [20,[46][47][48], we focussed on the analyses on the distinctions of the A loop region between the two systems. The positively charged sidechain of Arg912 in the A loop is involved in polar interactions with the aC-helix acidic residues of Glu768 and Asp771 (Figure 2A). In the nintedanib complex, Arg912 formed hydrogen bonds with either Glu768 or Asp771 for 75.15% of the simulation trajectories, while in the vandetanib complex the hydrogen bonds only took up 52.34% snapshots ( Figure 2B, C). Also, during the simulations of the RET G810A Ànintedanib complex, the distance between Arg912 and Glu768/Asp771 distributed a broader range than that of RET G810A Àvandetanib ( Figure 2B, C), which accounted for the observation that the A loop in the nintedanib system exhibited a higher flexibility from the RMSF analysis. Considering that a previous model suggested the auto-inhibited state of RET by the tether of the glycine-rich loop (GRL) and the aChelix to A loop, the interactions between A loop and aC-helix observed here may be associated with the higher inhibitory efficacy of nintedanib over vandetanib. Moreover, the disparities of the A loop RMSF may vary the activities of RET's auto-phosphorylation process, in which tyrosine residues in the A loop serve as the substrate for another RET molecule [16], further regulating its catalytic activity and substrate presenting ability.
Furthermore, the dynamic cross-correlation matrix (DCCM) was analyzed to reveal the correlations between residues during simulation. As shown in Figure 1D, E, regions in red denote correlation between the corresponding residues and regions in blue denotes anti-correlation. The regions with low density of colours indicate a low extent of correlation [49][50][51]. Herein, we focussed on the two loop regions with high flexibilities recognised by RMSF analysis. In the nintedanib complex, the loop between aD and aE helices was found to be anti-correlated with several other regions of RET, while the same region in the RET G810A Àvandetanib was poorly correlated. This explained the RMSF observation that the flexibility of corresponding residues was much lower in the nintedanib complex.

Nintedanib exhibits higher affinity towards RET G810A than vandetanib
Since nintedanib and vandetanib are both ATP-competitive inhibitors of RET, their binding affinities towards RET G810A are a key determinant for their pharmacological efficacies. Here, we exploited MM/PBSA calculations to compare the DG binding of two compounds towards RET G810A . For each system, a total of 180 frames from three independent replicates after the RMSD values reaching equilibrium were employed to calculate the average Gibb's free energies. As shown in Table 1, the binding free energies between nintedanib or vandetanib and the RET G810A are À46.92 ± 3.74 kcal/mol and À34.38 ± 3.65 kcal/mol, respectively, exhibiting significant differences. The $12 kcal/mol gap in the binding free   energy indicated a much higher affinity of nintedanib towards RET G810A , which suggested a greater ability for it to outcompete ATP and to inhibit a larger proportion cellular RET proteins than vandetanib.
The binding free energy was calculated by the algebraic summation of system energy changes on the inter-molecular mechanic interactions in the gas phase and the solutesolvent interactions upon the binding process. Each item was further divided into the polar (i.e. electrostatic) and nonpolar components, namely DE ele and DE vdw for the molecular mechanic interactions, and DE PB and DE nonpolar for the solvation energies, respectively. A detailed insight into the energetic components suggested that the predominant discrepancy between the two systems is the electrostatic attraction among the ligand and the RET G810A polar residues. This may mainly attribute to the more polarised methylcarboxy group and the indolin-2-one ring at the one end of the nintedanib molecule, which stretches deep into the ATP-binding pocket of RET G810A .

Combined structural and energetical analyses on nintedanib/ vandetanib À RET G810A interactions
To better understand the mechanisms for the higher affinity of nintedanib, we then gave a more detailed insights into the interaction modes between the inhibitors and RET G810A from both structural and energetical aspects. First, the representative structures of two systems during simulations were extracted using cluster analysis and the LigPlot program [52] was used to more explicitly display the interactions between the ligand's functional groups and the residues in RET G810A . As shown in Figure 3A,B, there were three hydrogen-bonds recognised in the representative structure of RET G810A Ànintedanib complex, while there was only one observed in the vandetanib complex. In addition, through hydrogen-bond analysis we revealed that the polar interactions between the ligand and RET G810A had a longer lifetime and shorter average hydrogen-bond length during simulations in the nintedanib-bound system than in the vandetanib-bound system ( Table 2). This suggested more stable and strengthened polar interactions between the inhibitor and RET in the RET G810A Ànintedanib complex, in consistence with the results from the above MM/PBSA analysis. Moreover, as more RET residues are recognised to participate in the hydrophobic contacts with the ligand, more binding interactions between nintedanib and RET G810A were observed ( Figure 3A, B).
Then, using the MM/PBSA decomposition method, we analyzed the residue-wise contributions on the binding affinity between the ligands and RET G810A ( Figure 3C). Residues that had significant different contributions (> 0.5 kcal/mol) between the two systems were shown as sticks in Figure 3D, E. As expected, the main contributors of the inhibitor-RET binding energy derived from the hinge strand (residues 806-811), which is also important for the binding of ATP [53,54]. This is because two residues from the hinge strand, Glu805 and Ala807, bolster the binding of nintedanib through hydrogen bonds with the N9 and O10 atoms of nintedanib, respectively. In contrast, the interaction was weaker between the hinge region and the vandetanib, since only one hydrogen-bond (between the N1 atom of vandetanib and RET Ala807) was detected ( Figure 3A, B). Moreover, Tyr806 of the hinge strand also had a significant higher binding contribution in the nintedanib system, which may be the result of a favoured hydrophobic interactions with the central phenyl ring of the nitedanib backbone.
In addition to the hinge region at the solvent front of the ATP-binding pocket, there was also a solvent-inaccessible residue Lys758 inside the deep ATP pocket that formed hydrogen bonds with nintedanib. Interestingly, the sidechain orientation of Lys758 is stabilised by the formation of hydrogen-bonds with Glu775 in the aC-helix and Asp892 in the DFG motif ( Figure 4). Lys758, Gly775, and Asp892 all favoured the binding of nintedanib over vandetanib, further confirming that the tether of these residues and the methylcarboxy group of nintedanib largely strengthened the docking of the inhibitor to the ATP-binding pocket. In contrast, no polar interactions comprising RET Lys758 were observed in the vandetanib complex throughout the simulation trajectories.   Additionally, we also found that at the mutation site the Ala810 blocked the polar interactions of the vandetanib towards the enzyme of the solvent by hydrophobic burying of the vandetanib O and O1 atoms ( Figure S1), while this effect was reduced to a lesser extent in the case of nintedanib. This finding can explain the drug resistance arose in the RET G810A mutant towards vandetanib. Together, these results strengthened the interpretation that nintedanib had a higher affinity towards RET G810A over vandetanib, which probably was the main reason that nintedanib persisted its efficacy to the vandetanibresistant G810A mutant. Apart from direct binding to the ATP pocket to inhibit ATP loading, orthosteric TKIs may also induce allosteric effects on specific regions of protein kinases, which further perturb the active-state conformations and diminish the binding of down-stream substrates or allosteric regulators [55][56][57][58][59][60]. To explore whether the different inhibitory effects on RET G810A of nintedanib and vandetanib are also a result of the allosteric regulation, a more detailed insight into the conformational dynamics of several specific regions of RET that is associated with RET activities was analyzed. The configurations of the DFG motif (residues 892-894), the aC-helix (residues 766-781), and A loop have been reported to be highly correlated to the activation process of protein kinases. In the initial crystal structures of RET À nintedanib and RET À vandetanib, the wild-type RET protein adopted a 'DFG-in, aC-helix-in, and A loop open' conformation ( Figure S2A-C), which suggested an active state of RET suitable for substrate binding and phosphorylation catalysis. We compared the conformational changes in these three regions in the nintedanib-and vandetanib-bound RET G810A mutant during the simulations.
According to the above RMSF analysis, the DFG motif and aC-helix were rigid in both systems during simulations. We observed that the DFG motif maintained in an 'in' conformation in both nintedanib-or vandetanib-bound RET G810A , with the Asp892 inserting into the ATP-binding pocket and the benzene ring of Phe893 pointing away from the pocket to form hydrophobic contacts with aC-helix ( Figure 5A). This was confirmed by the stabilisation of the heavy-atom RMSD of the DFG motif at low levels during simulations, using the initial crystal structures as references ( Figure 5B). Thus, no 'DFG-in' to 'DFG-out' conformational transitions was detected.
The 'aC-helix-in' conformation is characterised by the formation of a polar interaction between RET Glu775 of aC-helix and the Lys758 at its N-terminus b-sheet. As shown in Figure 4, the distance between CD atom of Glu775 and NZ atom of Lys758 was compactly distributed at a low range in both systems. Also, through hydrogen-bond analysis we observed polar interactions between Glu775 and Lys758 in 95.63% and 91.29% of the trajectory frames in the nintedanib-and vandetanib-bound systems, respectively. These results demonstrated that neither nintedanib nor vandetanib was able to disrupt the 'aChelix-in' conformation of RET G810A . In addition to enhance the binding of nintedanib towards RET G810A , the above-mentioned hydrogen-bond network featuring Lys758, Glu775, Asp892, and the methylcarboxy group of nintedanib also stabilised the 'DFG-in, aC-helix-in' conformation in RET G810A Ànintedanib (Figure 4).
A key determinant for the substrate binding ability of RET depends on whether the A loop is in the 'open' or 'closed' conformation. The 'open' conformation features an extended A loop, which is compatible with substrate docking, while in the 'closed' conformation A loop is compressed. For example, in the catalytically inactive state of the RET-analogous kinase EGFR, the A loop adopted a 'closed' conformation with the formation of two short a-helices in A loop. Here, using the dictionary of secondary structure of proteins (DSSP) algorithm [61], we validated that no a-helices was formed in the A loop in both systems during all simulations ( Additionally, the arrangement of the regulatory spine (R spine), consisting of four residues (Leu779, Leu790, His872 and Phe893), represents the activation status of the RET [62]. The RMSDs of the R spine also remained at a low value during simulations in both systems, implying the stabilisation of the active-state conformation ( Figure 5C, D). Moreover, in contrast to the previous crystal structure observations that the GRL adopted different conformations in the nintedanib-and vandetanib-bound RET [20], the GRL regions of RET G810A was both in the open conformation in our simulations ( Figure 5E), suggesting the active conformation.

Conclusions
Using MD simulations, we revealed the underlying mechanisms of the discrepant inhibitory effects of nintedanib and vandetanib on the kinase activities of RET G810A mutant. As ATP-competitive inhibitors, nintedanib exhibited a much higher affinity towards RET G810A ATP-binding pocket than vandetanib, which is due to the stronger hydrogenbond polar interactions between nintedanib and RET. Through H-bond analysis and energetic decomposition calculations, we found that nintedanib formed three stable hydrogen-bonds with RET G810A , involving RET Glu805 and Ala807 from the hinge region and Lys758 from the N-lobe. In comparison, there was only one hydrogen-bond observed between RET G810A and vandetanib, and was less durable in MD simulations. In the RET G810A Ànintedanib complex, the Lys758 was stabilised by the 'DFG-in, aC helix-in, A loop-open' conformation, where two residues from the aC helix and DFG motif (Glu775 and Asp892, respectively) tethered the sidechain of Lys758 and made it better accommodated to the binding of nintedanib. Moreover, the A loop in the RET G810A Ànintedanib complex displayed higher flexibility, which may also be associated with RET's auto-regulation and substrate presenting abilities. Collectively, these results untangled the molecular basis for the experimental findings that nintedanib was efficient to inhibit the vandetanib-resistant RET G810A mutant, inspiring future structure-based optimizations of TKIs targeting RET to overcome drug resistance. The alignment of the representative structures of RET G810A Ànintedanib (orange) and RET G810A Àvandetanib (green) displaying the DFG motif. The conformation of the DFG motif in the initial crystal structure of RET WT Ànintedanib is shown in white sticks. (B) Heavy-atom RMSD of the DFG motif relative to the initial crystal structures. Transparencies represent the standard deviations. (C) An alignment of the representative structures of RET G810A Ànintedanib and RET G810A Àvandetanib displaying R spine residues (Leu790, Leu779, Phe893, and His872). The R spine residues in the active assembly in the RET WT Ànintedanib crystal structure is also shown in white sticks. (D) Heavy-atom RMSD of the R spine residues relative to the initial crystal structures. Transparencies represent the standard deviations. (E) An alignment of the representative structures of RET G810A Ànintedanib and RET G810A Àvandetanib highlighting the GRL region. The GRL in RET G810A Ànintedanib and RET G810A Àvandetanib were both in the 'open' conformation, in comparison to the 'closed' conformation in the RET WT Ànintedanib crystal structure.