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 RET^G810A mutant. Recent studies have identified another TKI nintedanib as an effective molecule to inhibit vandetanib-resistant RET^G810A. However, there is no clear evidence of why nintedanib and vandetanib displayed different inhibitory activities towards RET^G810A. Here, we exploited molecular dynamic (MD) simulations to compare the interactions of the RET^G810A mutant with nintedanib and vandetanib. A higher structural flexibility of the activation loop was observed in the nintedanib-bound RET^G810A, 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 RET^G810A 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 RET^G810A 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 RET^G810A, 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.

Keywords:

• RET kinase
• molecular dynamics simulations
• drug resistance
• oncogenic mutations

1. 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–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–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–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.

2. Methods

2.1. 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.

2.2. 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 Cl^- counterions were added to neutralise system charges and then to the concentration of 0.154 mol·L⁻¹ to simulate the physiological environment. Second, for each system, two rounds of energy minimizations were carried out for structural optimizations [29–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 μs 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.

2.3. Dynamic cross-correlation matrix (DCCM) analysis

The dynamic cross-correlation matrix (DCCM) of all protein Cα atoms were calculated to reflect the inter-residue correlations [36–40]. The cross-correlation coefficient C_ij was calculated by: (1) $$C_{i,j}=\frac{〈(r_i-〈r_i〉〉·〈(r_j-〈r_j〉〉}{\sqrt{{〈(r_i-〈r_i〉〉}^2·{〈(r_j-〈r_j〉〉}^2}}$$(1) where r_i and r_j represents the positions of the ith and jth Cα atoms.

2.4. 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].

2.5. Determination of representative structures

The representative structures from certain groups of snapshots in the MD trajectories were extracted using the cluster analyses method [42–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 Cα atoms to eliminate overall rotations and transitions.

2.6. 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–48]. The Gibbs free energy change $\Delta G_{\mathit{binding}}$ of the binding process is defined by Eq. (2)(2) $$\Delta G_{\mathit{binding}}=G_{\mathit{complex}}-G_{\mathit{ligand}}-G_{\mathit{receptor}}$$(2) : (2) $$\Delta G_{\mathit{binding}}=G_{\mathit{complex}}-G_{\mathit{ligand}}-G_{\mathit{receptor}}$$(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_{\mathit{solv}}$) and the entropic components ($-TS$): (3) $$\Delta G_{\mathit{binding}}={(E}_{MM,\mathrm{}\mathit{complex}}-E_{MM,\mathrm{}\mathit{ligand}}-E_{MM,\mathrm{}\mathit{receptor}})+{(G}_{\mathit{solv},\mathrm{}\mathit{complex}}-G_{\mathit{solv},\mathrm{}\mathit{ligand}}-G_{\mathit{solv},\mathrm{}\mathit{receptor}})-({TS}_{\mathit{complex}}-{TS}_{\mathit{ligand}}-{TS}_{\mathit{receptor}})$$(3)

By defining molecular mechanical energy change ${\Delta E}_{MM},$ solvation energy change $\Delta G_{\mathit{solv}},$ and entropy change $\Delta S,$ Eq. (3)(3) $$\Delta G_{\mathit{binding}}={(E}_{MM,\mathrm{}\mathit{complex}}-E_{MM,\mathrm{}\mathit{ligand}}-E_{MM,\mathrm{}\mathit{receptor}})+{(G}_{\mathit{solv},\mathrm{}\mathit{complex}}-G_{\mathit{solv},\mathrm{}\mathit{ligand}}-G_{\mathit{solv},\mathrm{}\mathit{receptor}})-({TS}_{\mathit{complex}}-{TS}_{\mathit{ligand}}-{TS}_{\mathit{receptor}})$$(3) can be simplified to: (4) $$\Delta G_{\mathit{binding}}={\Delta E}_{MM}+\Delta G_{\mathit{solv}}-T\Delta S$$(4)

Then, ${\Delta E}_{MM}$ is ulteriorly decomposed into the van der Waals component ${\Delta E}_{\mathit{vdW}},$ the electrostatic component ${\Delta E}_{\mathit{ele}}$ and the internal component ${\Delta E}_{\mathit{int}}.$ Because ligands and RET proteins adopted same configurations in both unbound and bound states within a single snapshot, ${\Delta E}_{\mathit{int}}$ is always equal to zero in our calculations. (5) $$\Delta E_{MM}={\Delta E}_{\mathit{vdW}}+\Delta E_{\mathit{ele}}+\Delta E_{\mathit{int}}$$(5)

For the calculation of ${\Delta G}_{\mathit{solv}}$ in Eq. (4)(4) $$\Delta G_{\mathit{binding}}={\Delta E}_{MM}+\Delta G_{\mathit{solv}}-T\Delta S$$(4) , we divided it into the polar part ${\Delta E}_{PB}$ and the nonpolar part ${\Delta E}_{\mathit{nonpolar}}$ and applied Poisson-Boltzmann continuum solvent model: (6) $$\Delta G_{\mathit{solv}}={\Delta E}_{PB}+\Delta E_{\mathit{nonpolar}}$$(6)

The nonpolar component ${\Delta E}_{\mathit{nonpolar}}$ was calculated based on the solvent-accessible surface-area (SASA) with the following format, where the surface tension parameter γ was set to 0.00542 kcal·mol⁻¹·Å⁻², and the solvation parameter b was set to 0.92 kcal/mol. (7) $$\Delta E_{\mathit{nonpolar}}=\gamma \mathit{SASA}+b$$(7)

The conformation entropy component ($-T\Delta S$) 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. Results

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 vandetanib-resistant 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.

Figure 1. Conformational dynamics of RET^G810A in complex with nintedanib and vandetanib. (A) The alignment of the initial structures of RET^G810A−nintedanib and RET^G810A−vandetanib based on PDB ID 6NEC (chain A) and 2IVU. (B) RMSDs of all Cα atoms of RET^G810A bound to nintedanib (orange) and vandetanib (green) along simulation time. (C) RMSFs of RET^G810A bound to nintedanib (orange) and vandetanib (green) during equilibrium simulations. Transparencies represent the standard deviations. Dynamic cross-correlation matrix of (D) RET^G810A−nintedanib and (E) RET^G810A−vandetanib. Red colour denotes correlation and blue colour denotes anti-correlation. The density of colours implies the correlation extent.

The root-mean-square deviations (RMSDs) of RET^G810A C_α 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_α 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 αD and αE 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–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 αC-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 αC-helix to A loop, the interactions between A loop and αC-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.

Figure 2. Interactions between RET A loop and αC-helix. (A) The alignment of the representative structures of RET^G810A−nintedanib (orange) and RET^G810A−vandetanib (green). Hydrogen-bonds in two systems are shown by cyan and yellow dashed lines, respectively. (B) Distance distributions between Glu768 CG atom and Arg912 CZ atom in two systems. (C) Distance distribution between Asp771 CG atom and Arg912 CZ atom in two systems. H-bond frac.% represents the fractions of frames in which the corresponding hydrogen-bond is detected.

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–51]. Herein, we focussed on the two loop regions with high flexibilities recognised by RMSF analysis. In the nintedanib complex, the loop between αD and αE 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.

3.2. 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 $\Delta G_{\mathit{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.

Table 1. Binding free energy analysis between nintedanib/vandetanib and RET^G810A (all units in kcal/mol)*.

	RET^G810A−nintedanib	RET^G810A−vandetanib
$\Delta E_{\mathit{vdw}}a$	−57.3585 (2.8745)	−53.7228 (3.0229)
$\Delta E_{\mathit{ele}}b$	−27.6610 (6.6645)	−18.9090 (4.2395)
$\Delta E_{PB}c$	43.2488 (6.1599)	42.6741 (4.6422)
$\Delta E_{\mathit{nonpolar}}d$	−5.1445 (0.1979)	−4.4240 (0.1820)
$\Delta E_{MM}e$	−85.0194 (7.4984)	−72.6318 (5.7238)
$\Delta G_{\mathit{solv}}f$	38.1044 (6.0392)	38.2501 (4.5701)
$\Delta {\mathit{G}}_{\mathit{binding}}$	−46.9151 (3.7371)	−34.3817 (3.6543)

*a. The van der Waals component of the total molecular mechanical energy change. b. the electrostatic component of the total molecular mechanical energy change. c. The electrostatic component of solvation energy change determined by Poisson Boltzmann (PB) equation. d. The nonpolar component of solvation energy change. e. Total molecular mechanical energy change. f. Total solvation energy change. All numbers in the parentheses represent standard deviations.

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 solute-solvent interactions upon the binding process. Each item was further divided into the polar (i.e. electrostatic) and nonpolar components, namely $\Delta E_{\mathit{ele}}$ and $\Delta E_{\mathit{vdw}}$ for the molecular mechanic interactions, and $\Delta E_{PB}$ and $\Delta E_{\mathit{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.

3.3. 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).

Figure 3. Structural and energetic comparisons on the interactions of nintedanib/vandetanib with RET^G810A. LigPlot diagram of the (A) RET^G810A−nintedanib complex and (B) RET^G810A−vandetanib complex. Hydrogen bonds between ligands and protein residues are represented by green dashed lines and bond length are labelled in green (Å). Hydrophobic contacts are represented by an arc with spokes radiating towards the ligand atoms they contact. The contacted atoms in the ligands are shown with spokes radiating back. (C) Per-residue decomposition of binding free energy for the RET^G810A−nintedanib (orange) and RET^G810A−vandetanib (green) systems. (D) Front and (E) back view of residues with large binding free energy contributions (>0.5 kcal/mol) in the representative structures. Residues that favour the binding of nintedanib are labelled in red while residues favour the binding of vandetanib are labelled in blue.

Table 2. Hydrogen-bonds between RET^G810A and the inhibitors*.

System	Acceptor^a	Donor^b	Fraction^c	AvgDist^d	AvgAng^e
RET^G810A−Nintedanib	NINT O10	Ala807	53.15%	2.84	159.8
RET^G810A−Nintedanib	NINT O26	Lys788	42.00%	2.84	158.0
RET^G810A−Nintedanib	Glu805	NINT HN9	21.71%	2.88	165.1
RET^G810A−Vandetanib	VAND N1	Ala807	26.78%	2.92	157.6

*a. Hydrogen-acceptor in the hydrogen-bond. b. Hydrogen-donor in the hydrogen-bond. c. The fraction of frames in the trajectories in which hydrogen-bond is detected. d. Average distance (Å) of heavy atoms involved in the hydrogen-bonds. e. Average bond angle (°) of the hydrogen-bonds.

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 αC-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.

Figure 4. The representative structures of the hydrogen-bond networks involving Lys758 in (A) RET^G810A−nintedanib and (B) RET^G810A−vandetanib complex. Distance distributions between (C) Lys758 NZ atom and Glu775 CD atom, (D) Lys758 NZ atom and Asp892 CG atom in two systems. (E) Distance distribution between Lys758 NZ atom and nintedanib O26 atom in the RET^G810A−nintedanib complex. H-bond frac.% represents the fractions of frames in which the corresponding hydrogen-bond is detected.

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 vandetanib-resistant G810A mutant.

3.4. Both nintedanib- and vandetanib-bound RET^G810A maintained the active ‘DFG-in, αC-helix-in, a loop-open’ conformation

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–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 αC-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, αC-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 αC-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 αC-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.

Figure 5. Detailed conformational insights into several RET functional regions. (A) 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.

The ‘αC-helix-in’ conformation is characterised by the formation of a polar interaction between RET Glu775 of αC-helix and the Lys758 at its N-terminus β-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 ‘αC-helix-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, αC-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 α-helices in A loop. Here, using the dictionary of secondary structure of proteins (DSSP) algorithm [61], we validated that no α-helices was formed in the A loop in both systems during all simulations (Table S1 & S2). This confirmed that despite the disparities in the A loop flexibility, it persisted the ‘open’ state conformation in both systems.

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.

4. 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 hydrogen-bond 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, αC helix-in, A loop-open’ conformation, where two residues from the αC 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.

Acknowledgements

The authors thank the Supercomputer Center at Shanghai Supercenter for providing computational resources.

Disclosure statement

The authors have no conflicts of interest.