L1198F Mutation Resensitizes Crizotinib to ALK by Altering the Conformation of Inhibitor and ATP Binding Sites

The efficacy of anaplastic lymphoma kinase (ALK) positive non-small-cell lung cancer (NSCLC) treatment with small molecule inhibitors is greatly challenged by acquired resistance. A recent study reported the newest generation inhibitor resistant mutation L1198F led to the resensitization to crizotinib, which is the first Food and Drug Administration (FDA) approved drug for the treatment of ALK-positive NSCLC. It is of great importance to understand how this extremely rare event occurred for the purpose of overcoming the acquired resistance of such inhibitors. In this study, we exploited molecular dynamics (MD) simulation to dissect the molecular mechanisms. Our MD results revealed that L1198F mutation of ALK resulted in the conformational change at the inhibitor site and altered the binding affinity of ALK to crizotinib and lorlatinib. L1198F mutation also affected the autoactivation of ALK as supported by the identification of His1124 and Tyr1278 as critical amino acids involved in ATP binding and phosphorylation. Our findings are valuable for designing more specific and potent inhibitors for the treatment of ALK-positive NSCLC and other types of cancer.


Introduction
Lung cancer is one of the most life-threatening malignancies worldwide and is the leading cause of cancer-related death for both men and women [1]. Anaplastic lymphoma kinase (ALK) is a member of the receptor tyrosine kinases (RTKs), which belong to the insulin receptor kinase superfamily [2]. Chromosomal rearrangements in the ALK gene lead to the deregulation of ALK kinase activity, which in turn alters the downstream signaling pathways in cancer biology [3]. Abnormal expression of fused ALK genes has been implicated in the pathogenesis of several types of cancer, including non-small-cell lung cancer (NSCLC), anaplastic large-cell lymphoma, glioblastoma,

Root-Mean-Square Deviation Analysis of the Protein Backbones in Crizotinib/Lorlatinib Associated ALKs
We performed molecular dynamics simulation of the ALK-inhibitor complexes for 30 ns with GROMACS software. We first analyzed the root-mean-square deviation (RMSD) of the protein backbones in crizotinib or lorlatinib associated wild type, C1156Y, L1198F, and C1156Y-L1198F mutants. As shown in Figure 1A, the RMSD of ALK-crizotinib complexes quickly reached a steady state after 5 ns of simulation. The fluctuation of the wild type ALK was slightly higher than the other mutants. The C1156Y-L1198F mutant experienced a leap of RMSD up to 0.2 nm from around 20 ns to 25 ns. Up to the end of the 30 ns simulation, the RMSD of C1156Y, L1198F, and C1156Y-L1198F mutants were steady around 0.15 nm, while the value of the wild type protein was moderately higher than the others mutants. The RMSD of ALK mutants complexed with lorlatinib were fairly stable throughout the whole course of simulation ( Figure 1B). There was no significant difference among the protein backbones analyzed. According these results, crizotinib or lorlatinib association did not significantly affect the RMSD protein backbones.

Crizontinib and Lorlatinib Binds to ALK (C1156Y-L1198F) with Different Affinity
In order to elucidate the mechanism that rendered the resensitization of ALK (C1156Y-L1198F) to crizotinib, we decided to focus on comparing ALK-crizotinib and ALK-lorlatinib complexes. We will refer to ALK (C1156Y-L1198F) double mutants simply as ALK in the future discussion. Structurally, crizotinib and lorlatinib bond to the same pocket of ALK protein (Figure 2A,B). The major difference is that lorlatinib is supposed to have a higher selectivity, which is achieved by the targeting of L1198 presented in only 25% of the kinases. Binding energy analysis of ALK-inhibitor complex by AUTODOCK software revealed that crizotinib binds to ALK with a much lower energy (25.9 kJ/mol) compared with lorlatinib (123.7 kJ/mol) (Table S1). This data indicated that crizotinib bound to ALK protein tighter than lolatinib. The prediction was in line with the inhibition constant (Ki) and half-maximal inhibitory concentration (IC50) values presented by Shaw et al. [22]. Further, energy decomposition by Hawkins generalized Born surface area (Hawkins GB/SA), molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) binding energy analysis also supported this conclusion (Table S1). As shown in Figure 3, the predicted electrostatic energy of ALK-crizotinib complex using the above-mentioned methods were 4.8-, 6.1-, and 8.7-folds lower than the ALK-lorlatinib complex, suggesting a higher affinity of crizotinib to the ALK (C1156Y-L1198F) mutant protein than lorlatinib.

Crizontinib and Lorlatinib Binds to ALK (C1156Y-L1198F) with Different Affinity
In order to elucidate the mechanism that rendered the resensitization of ALK (C1156Y-L1198F) to crizotinib, we decided to focus on comparing ALK-crizotinib and ALK-lorlatinib complexes. We will refer to ALK (C1156Y-L1198F) double mutants simply as ALK in the future discussion. Structurally, crizotinib and lorlatinib bond to the same pocket of ALK protein (Figure 2A,B). The major difference is that lorlatinib is supposed to have a higher selectivity, which is achieved by the targeting of L1198 presented in only 25% of the kinases.

Crizontinib and Lorlatinib Binds to ALK (C1156Y-L1198F) with Different Affinity
In order to elucidate the mechanism that rendered the resensitization of ALK (C1156Y-L1198F) to crizotinib, we decided to focus on comparing ALK-crizotinib and ALK-lorlatinib complexes. We will refer to ALK (C1156Y-L1198F) double mutants simply as ALK in the future discussion. Structurally, crizotinib and lorlatinib bond to the same pocket of ALK protein (Figure 2A,B). The major difference is that lorlatinib is supposed to have a higher selectivity, which is achieved by the targeting of L1198 presented in only 25% of the kinases. Binding energy analysis of ALK-inhibitor complex by AUTODOCK software revealed that crizotinib binds to ALK with a much lower energy (25.9 kJ/mol) compared with lorlatinib (123.7 kJ/mol) (Table S1). This data indicated that crizotinib bound to ALK protein tighter than lolatinib. The prediction was in line with the inhibition constant (Ki) and half-maximal inhibitory concentration (IC50) values presented by Shaw et al. [22]. Further, energy decomposition by Hawkins generalized Born surface area (Hawkins GB/SA), molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) binding energy analysis also supported this conclusion (Table S1). As shown in Figure 3, the predicted electrostatic energy of ALK-crizotinib complex using the above-mentioned methods were 4.8-, 6.1-, and 8.7-folds lower than the ALK-lorlatinib complex, suggesting a higher affinity of crizotinib to the ALK (C1156Y-L1198F) mutant protein than lorlatinib. Binding energy analysis of ALK-inhibitor complex by AUTODOCK software revealed that crizotinib binds to ALK with a much lower energy (25.9 kJ/mol) compared with lorlatinib (123.7 kJ/mol) (Table S1). This data indicated that crizotinib bound to ALK protein tighter than lolatinib. The prediction was in line with the inhibition constant (Ki) and half-maximal inhibitory concentration (IC 50 ) values presented by Shaw et al. [22]. Further, energy decomposition by Hawkins generalized Born surface area (Hawkins GB/SA), molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) binding energy analysis also supported this conclusion (Table S1). As shown in Figure 3, the predicted electrostatic energy of ALK-crizotinib complex using the above-mentioned methods were 4.8-, 6.1-, and 8.7-folds lower than the ALK-lorlatinib complex, suggesting a higher affinity of crizotinib to the ALK (C1156Y-L1198F) mutant protein than lorlatinib.

Root-Mean-Square Deviation Comparison of ALK-Crizotinib and ALK-Lorlatinib Complexes
ALK protein in the ALK-crizotinib and ALK-lorlatinib complexes were comparably stable during the whole course of 30 ns simulation, except for that ALK-crizotinib complex fluctuated during 20-25 ns ( Figure 4A). Since this fluctuation only lasted for a relatively short period of time, it might result from the adjustment of the binding pocket of ALK to the crizotinib in the process of protein-ligand association. This was supported by the RMSD analysis of the two ligands. We found that the RMSD of crizotinib was fluctuating violently throughout the course of simulation ( Figure 4B), indicating that the drug binding pocket of ALK was reorganized considerably upon the association of crizotinib. Such reorganization was confined to the binding pocket of crizotinib, because the ALK protein as a whole remained relatively steady. The RMSD of lorlatinib was extremely stable from the start to the end ( Figure 4B). Generally, a narrower fluctuation of the RMSD means the system is more stable. The RMSD analysis results strongly indicated that the binding affinity change resulted from L1198F mutation was not enough to explain the observed drug sensitivity shift. The deregulation of the other key events, such as ATP association and the substrate phosphorylation, might also contribute to the change of drug sensitivity.

Identification of the Key Amino Acid Residues Affecting ALK Activity
In order to determine the key amino acid residues that are affecting the ALK activity upon inhibitor binding, we investigated the electrostatic energy change trends of ALK-crizotinib and ALK-lorlatinib complexes for three periods of time in MD simulation (5-10 ns, 20-25 ns, 25-30 ns) (Table S2). Only the amino acid residues that kept the same trend of energy change, either constantly

Root-Mean-Square Deviation Comparison of ALK-Crizotinib and ALK-Lorlatinib Complexes
ALK protein in the ALK-crizotinib and ALK-lorlatinib complexes were comparably stable during the whole course of 30 ns simulation, except for that ALK-crizotinib complex fluctuated during 20-25 ns ( Figure 4A). Since this fluctuation only lasted for a relatively short period of time, it might result from the adjustment of the binding pocket of ALK to the crizotinib in the process of protein-ligand association. This was supported by the RMSD analysis of the two ligands. We found that the RMSD of crizotinib was fluctuating violently throughout the course of simulation ( Figure 4B), indicating that the drug binding pocket of ALK was reorganized considerably upon the association of crizotinib. Such reorganization was confined to the binding pocket of crizotinib, because the ALK protein as a whole remained relatively steady. The RMSD of lorlatinib was extremely stable from the start to the end ( Figure 4B). Generally, a narrower fluctuation of the RMSD means the system is more stable. The RMSD analysis results strongly indicated that the binding affinity change resulted from L1198F mutation was not enough to explain the observed drug sensitivity shift. The deregulation of the other key events, such as ATP association and the substrate phosphorylation, might also contribute to the change of drug sensitivity.

Root-Mean-Square Deviation Comparison of ALK-Crizotinib and ALK-Lorlatinib Complexes
ALK protein in the ALK-crizotinib and ALK-lorlatinib complexes were comparably stable during the whole course of 30 ns simulation, except for that ALK-crizotinib complex fluctuated during 20-25 ns ( Figure 4A). Since this fluctuation only lasted for a relatively short period of time, it might result from the adjustment of the binding pocket of ALK to the crizotinib in the process of protein-ligand association. This was supported by the RMSD analysis of the two ligands. We found that the RMSD of crizotinib was fluctuating violently throughout the course of simulation ( Figure 4B), indicating that the drug binding pocket of ALK was reorganized considerably upon the association of crizotinib. Such reorganization was confined to the binding pocket of crizotinib, because the ALK protein as a whole remained relatively steady. The RMSD of lorlatinib was extremely stable from the start to the end ( Figure 4B). Generally, a narrower fluctuation of the RMSD means the system is more stable. The RMSD analysis results strongly indicated that the binding affinity change resulted from L1198F mutation was not enough to explain the observed drug sensitivity shift. The deregulation of the other key events, such as ATP association and the substrate phosphorylation, might also contribute to the change of drug sensitivity.

Identification of the Key Amino Acid Residues Affecting ALK Activity
In order to determine the key amino acid residues that are affecting the ALK activity upon inhibitor binding, we investigated the electrostatic energy change trends of ALK-crizotinib and ALK-lorlatinib complexes for three periods of time in MD simulation (5-10 ns, 20-25 ns, 25-30 ns) (Table S2). Only the amino acid residues that kept the same trend of energy change, either constantly

Identification of the Key Amino Acid Residues Affecting ALK Activity
In order to determine the key amino acid residues that are affecting the ALK activity upon inhibitor binding, we investigated the electrostatic energy change trends of ALK-crizotinib and ALK-lorlatinib complexes for three periods of time in MD simulation (5-10 ns, 20-25 ns, 25-30 ns) (Table S2). Only the amino acid residues that kept the same trend of energy change, either constantly increasing or decreasing, were likely to affect the activity of ALK (i.e., contributing to the drug sensitivity change). According to this criterion, we identified 174 amino acid residues (Table S2), among which His1124, Lys1150, Met1199, Asp1203, and Glu1210 were the residues involved in the domains of inhibitor and ATP binding ( Figure 5A). We calculated the folds of energy change for the five key residues identified (Table S3). The average folds of the relatively energy change for His1124, Lys1150, Met1199, Asp1203, Glu1210 were 5.60, 1.14, 1.14, 7.86, and 8.55, respectively ( Figure 5B). increasing or decreasing, were likely to affect the activity of ALK (i.e., contributing to the drug sensitivity change). According to this criterion, we identified 174 amino acid residues (Table S2), among which His1124, Lys1150, Met1199, Asp1203, and Glu1210 were the residues involved in the domains of inhibitor and ATP binding ( Figure 5A). We calculated the folds of energy change for the five key residues identified (Table S3). The average folds of the relatively energy change for His1124, Lys1150, Met1199, Asp1203, Glu1210 were 5.60, 1.14, 1.14, 7.86, and 8.55, respectively ( Figure 5B).

Crizotinib and Lorlatinib Interacts with ALK in Different Modes
We generated the two-dimensional (2D) interaction diagram of ALK-crizotinib and ALK-lorlatinib to analyze the hydrogen bonds formation and hydrophobic interaction between ALK and the small molecule drug ( Figure 6). The N22 and N23 of crizotinib formed two hydrogen bonds with Glu1197 and Met1199 of ALK. The distances of the two hydrogen bonds were 2.96 Å and 3.02 Å, respectively. Lorlatinib interacted with Glu1197 and Met1199 of ALK by forming three hydrogen bonds through the N3, N17, and N24. The distances were 2.96 Å, 2.81 Å, and 3.58 Å, respectively. In terms of hydrophobic interaction, amino acid residues Leu1122, Ala1148, Leu1196, Ala1200, Gly1202, and Leu1256 were conserved between crizotinib and lorlatinib. Lys1150 and Arg1253 were only involved in ALK-crizotinib interaction, while Gly1123, Val1130, Phe1198, and Gly1269 were unique to ALK-lorlatinib interaction.

Crizotinib and Lorlatinib Interacts with ALK in Different Modes
We generated the two-dimensional (2D) interaction diagram of ALK-crizotinib and ALK-lorlatinib to analyze the hydrogen bonds formation and hydrophobic interaction between ALK and the small molecule drug ( Figure 6). The N22 and N23 of crizotinib formed two hydrogen bonds with Glu1197 and Met1199 of ALK. The distances of the two hydrogen bonds were 2.96 Å and 3.02 Å, respectively. Lorlatinib interacted with Glu1197 and Met1199 of ALK by forming three hydrogen bonds through the N3, N17, and N24. The distances were 2.96 Å, 2.81 Å, and 3.58 Å, respectively. In terms of hydrophobic interaction, amino acid residues Leu1122, Ala1148, Leu1196, Ala1200, Gly1202, and Leu1256 were conserved between crizotinib and lorlatinib. Lys1150 and Arg1253 were only involved in ALK-crizotinib interaction, while Gly1123, Val1130, Phe1198, and Gly1269 were unique to ALK-lorlatinib interaction. increasing or decreasing, were likely to affect the activity of ALK (i.e., contributing to the drug sensitivity change). According to this criterion, we identified 174 amino acid residues (Table S2), among which His1124, Lys1150, Met1199, Asp1203, and Glu1210 were the residues involved in the domains of inhibitor and ATP binding ( Figure 5A). We calculated the folds of energy change for the five key residues identified (Table S3). The average folds of the relatively energy change for His1124, Lys1150, Met1199, Asp1203, Glu1210 were 5.60, 1.14, 1.14, 7.86, and 8.55, respectively ( Figure 5B).

Crizotinib and Lorlatinib Interacts with ALK in Different Modes
We generated the two-dimensional (2D) interaction diagram of ALK-crizotinib and ALK-lorlatinib to analyze the hydrogen bonds formation and hydrophobic interaction between ALK and the small molecule drug ( Figure 6). The N22 and N23 of crizotinib formed two hydrogen bonds with Glu1197 and Met1199 of ALK. The distances of the two hydrogen bonds were 2.96 Å and 3.02 Å, respectively. Lorlatinib interacted with Glu1197 and Met1199 of ALK by forming three hydrogen bonds through the N3, N17, and N24. The distances were 2.96 Å, 2.81 Å, and 3.58 Å, respectively. In terms of hydrophobic interaction, amino acid residues Leu1122, Ala1148, Leu1196, Ala1200, Gly1202, and Leu1256 were conserved between crizotinib and lorlatinib. Lys1150 and Arg1253 were only involved in ALK-crizotinib interaction, while Gly1123, Val1130, Phe1198, and Gly1269 were unique to ALK-lorlatinib interaction.  We also calculated the distance of the inhibitor to the hydrogen bonds forming amino acid residues, Glu1197 and Met 1199, with respect to time (Table 1). During the whole process of simulation, both Glu1197 and Met1199 fluctuated significantly in the ALK-crizotinib and ALK-lorlatinib complexes. This result supported that L1198F mutation caused dramatic conformational changes to the inhibitor binding pocket, especially the amino acid residues nearby.

Root-Mean-Square Fluctuation Analysis of the Key Amino Acid Residues
We then analyzed the trends of root-mean-square fluctuation (RMSF) change for the amino acids that participated in ATP binding, inhibitor binding, proton binding, and phosphorylation ( Figure 7). The RMSF values of His1124 and Glu1210 in ALK-crizotinib were constantly higher than that of ALK-lorlatinib complex throughout the simulation (Table S4). The average folds of change were 1.80 for His1124 and 1.29 for Glu1210 (Table S4). His1124 is involved in ATP binding and Glu1210 is involved in inhibitor binding, which are the two key events affecting the physiological outcome of treatment targeting ALK. We also calculated the distance of the inhibitor to the hydrogen bonds forming amino acid residues, Glu1197 and Met 1199, with respect to time (Table 1). During the whole process of simulation, both Glu1197 and Met1199 fluctuated significantly in the ALK-crizotinib and ALK-lorlatinib complexes. This result supported that L1198F mutation caused dramatic conformational changes to the inhibitor binding pocket, especially the amino acid residues nearby.

Root-Mean-Square Fluctuation Analysis of the Key Amino Acid Residues
We then analyzed the trends of root-mean-square fluctuation (RMSF) change for the amino acids that participated in ATP binding, inhibitor binding, proton binding, and phosphorylation ( Figure 7). The RMSF values of His1124 and Glu1210 in ALK-crizotinib were constantly higher than that of ALK-lorlatinib complex throughout the simulation (Table S4). The average folds of change were 1.80 for His1124 and 1.29 for Glu1210 (Table S4). His1124 is involved in ATP binding and Glu1210 is involved in inhibitor binding, which are the two key events affecting the physiological outcome of treatment targeting ALK.

L1198F Mutation Affects the ATP Association of ALK
Both the electrostatic energy ( Figure 5) and RMSF (Figure 7) change trends analysis identified His1124 as a key amino acid involved in regulating the activity of ALK. Since ALK is an abnormally expressed kinase in many types of cancers, it is not surprising to see His1124 that is a key amino acid residue participating in ATP binding functions in the resensitization of crizotinib. A structural study by Bossi et al. suggests that the O2 oxygen of the α-phosphate is within hydrogen-bonding distance of the backbone carbonyl of His1124 of the P-loop [23]. A modeling study of the ALK Gly1123-His1124 segment indicated that mutations in this part of the protein are likely to sterically impede ATP binding [24]. We further analyzed the interaction of ALK and ATP through RMSD and RMSF. Because the inhibitor binding site and the ATP binding site are very close, we were not able

L1198F Mutation Affects the ATP Association of ALK
Both the electrostatic energy ( Figure 5) and RMSF (Figure 7) change trends analysis identified His1124 as a key amino acid involved in regulating the activity of ALK. Since ALK is an abnormally expressed kinase in many types of cancers, it is not surprising to see His1124 that is a key amino acid residue participating in ATP binding functions in the resensitization of crizotinib. A structural study by Bossi et al. suggests that the O2 oxygen of the α-phosphate is within hydrogen-bonding distance of the backbone carbonyl of His1124 of the P-loop [23]. A modeling study of the ALK Gly1123-His1124 segment indicated that mutations in this part of the protein are likely to sterically impede ATP binding [24]. We further analyzed the interaction of ALK and ATP through RMSD and RMSF. Because the inhibitor binding site and the ATP binding site are very close, we were not able to perform a MD simulation of ALK-inhibitor-ATP complex. Instead, we isolated the ALK protein from ALK-crizotinib and ALK-lorlatinib complex and combined them with ATP molecule to run the simulation. The fluctuation of ALK in the crizotinib-associated complex was slightly wider than that of lorlatinib ( Figure 8A). The ATP in the two complexes behaved differently ( Figure 8B). In the crizotinib-associated complex, the average value of ATP fluctuation was 0.21 nm, in comparison with 0.29 nm of the lorlatinib-associated complex. The largest fluctuation in the crizotinib-associated complex was 0.38 nm, in comparison with 0.42 nm in the lorlatinib-associated complexes. These differences indicated that ATP binding was less stable in the ALK-lorlatinib complex. to perform a MD simulation of ALK-inhibitor-ATP complex. Instead, we isolated the ALK protein from ALK-crizotinib and ALK-lorlatinib complex and combined them with ATP molecule to run the simulation. The fluctuation of ALK in the crizotinib-associated complex was slightly wider than that of lorlatinib ( Figure 8A). The ATP in the two complexes behaved differently ( Figure 8B). In the crizotinib-associated complex, the average value of ATP fluctuation was 0.21 nm, in comparison with 0.29 nm of the lorlatinib-associated complex. The largest fluctuation in the crizotinib-associated complex was 0.38 nm, in comparison with 0.42 nm in the lorlatinib-associated complexes. These differences indicated that ATP binding was less stable in the ALK-lorlatinib complex. We then analyzed the RMSF change between the two complexes ( Figure 9). After 60 ns of simulation, His1124 was identified as the residue that gave the largest RMSF value difference (0.27 nm) between the crizotinib-and lorlatinib-associated ALK protein (Table S5). This result confirmed the importance of His1124 in resensitizing crizotinib to C1156Y-L1198F mutated ALK protein. In addition, the RMSF analysis also found that Arg1279 gave the second largest RMSF value difference (0.24 nm) and Tyr1278 showed a 0.11 nm RMSF difference between the two complexes compared (Table S5). Both Tyr1278 and Arg1279 are in the YRASYY sequence that is present in the activation loop of the kinase domain, and Tyr1278 has been defined as the first tyrosine to be phosphorylated in this sequence [25]. Presumably, the L1198F mutation led to conformational changes in the ATP binding pockets, including His1124. Such changes decrease the binding affinity of ATP and hinder the autoactivation of ALK regulated by the phosphorylation on Tyr1278. These events are likely to affect the downstream phosphorylation of ALK substrates. The critical amino acids identified in our MD analysis, including His1124, Lys1150, Met1199, Asp1203, Glu1210, and Tyr1278, were in line with the experimental results by another group (Table 2). These amino acid residues are involved in inhibitor binding (Lys1150, Met1199, Asp1203, and Glu1210) [24,26], ATP binding (His1124) [24,25], and the phosphorylation (Tyr1278) of the ALK protein [27]. We then analyzed the RMSF change between the two complexes ( Figure 9). After 60 ns of simulation, His1124 was identified as the residue that gave the largest RMSF value difference (0.27 nm) between the crizotinib-and lorlatinib-associated ALK protein (Table S5). This result confirmed the importance of His1124 in resensitizing crizotinib to C1156Y-L1198F mutated ALK protein. In addition, the RMSF analysis also found that Arg1279 gave the second largest RMSF value difference (0.24 nm) and Tyr1278 showed a 0.11 nm RMSF difference between the two complexes compared (Table S5). Both Tyr1278 and Arg1279 are in the YRASYY sequence that is present in the activation loop of the kinase domain, and Tyr1278 has been defined as the first tyrosine to be phosphorylated in this sequence [25]. Presumably, the L1198F mutation led to conformational changes in the ATP binding pockets, including His1124. Such changes decrease the binding affinity of ATP and hinder the autoactivation of ALK regulated by the phosphorylation on Tyr1278. These events are likely to affect the downstream phosphorylation of ALK substrates. The critical amino acids identified in our MD analysis, including His1124, Lys1150, Met1199, Asp1203, Glu1210, and Tyr1278, were in line with the experimental results by another group (Table 2). These amino acid residues are involved in inhibitor binding (Lys1150, Met1199, Asp1203, and Glu1210) [24,26], ATP binding (His1124) [24,25], and the phosphorylation (Tyr1278) of the ALK protein [27].

Secondary Structure Analysis of the ATP and Inhibitor Binding Sites
The activity of a protein is often affected by its conformational changes. We analyzed the secondary structure change of the amino acids around His1124 and Glu1210. The secondary structure of the analyzed region (aa1122-1134) was relatively steady in the first 5 ns of the simulation and His1124 was located in a turn ( Figure 10A, panels a and d). At 20-25 ns, the structure of the corresponding domain in crizotinib-associated ALK was significantly different from that of the lorlatinib-associated protein ( Figure 10A, panels b and e). The change was represented by conversion from turns to coils. This difference was still observable at the end of the simulation ( Figure 10A, panels c and f). As shown in Figure 10B, the secondary structure of aa1203-1213 was mainly composed of -helixes and coils. There was no noticeable difference at residues 1203 to 1213 between ALK-crizotinib and ALK-lorlatinib complexes at the beginning of simulation (5-10 ns). During 20-25 ns, the secondary structure of ALK-crizotinib dramatically shifted from α-helixes to coils ( Figure 10B, panels b and f), which was not the case for ALK-lorlatinib. The secondary structure of ALK-lorlatinib remained to be fairly steady throughout the simulation. The predicted secondary structure difference was consistent with the RMSD of crizotinib and lorlatinib, as shown in Figure 4B, indicating dramatic conformational shifts of the ALK protein at the crucial domains to accommodate crizotinib binding.

Secondary Structure Analysis of the ATP and Inhibitor Binding Sites
The activity of a protein is often affected by its conformational changes. We analyzed the secondary structure change of the amino acids around His1124 and Glu1210. The secondary structure of the analyzed region (aa1122-1134) was relatively steady in the first 5 ns of the simulation and His1124 was located in a turn ( Figure 10A, panels a and d). At 20-25 ns, the structure of the corresponding domain in crizotinib-associated ALK was significantly different from that of the lorlatinib-associated protein ( Figure 10A, panels b and e). The change was represented by conversion from turns to coils. This difference was still observable at the end of the simulation ( Figure 10A, panels c and f). As shown in Figure 10B, the secondary structure of aa1203-1213 was mainly composed of α-helixes and coils. There was no noticeable difference at residues 1203 to 1213 between ALK-crizotinib and ALK-lorlatinib complexes at the beginning of simulation (5-10 ns). During 20-25 ns, the secondary structure of ALK-crizotinib dramatically shifted from α-helixes to coils ( Figure 10B, panels b and f), which was not the case for ALK-lorlatinib. The secondary structure of ALK-lorlatinib remained to be fairly steady throughout the simulation. The predicted secondary structure difference was consistent with the RMSD of crizotinib and lorlatinib, as shown in Figure 4B, indicating dramatic conformational shifts of the ALK protein at the crucial domains to accommodate crizotinib binding.

Molecular Docking
The molecular graphics of ALK-inhibitor complexes were prepared and analyzed with the University of California, San Francisco (UCSF) Chimera package [28]. In this process, (i) solvent and non-complexed ions were removed; and (ii) hydrogens and charges (of amber ff99sb force field) were added to the protein [28]. The docking was performed with DOCK 6.5 program (UCSF) under the

Molecular Docking
The molecular graphics of ALK-inhibitor complexes were prepared and analyzed with the University of California, San Francisco (UCSF) Chimera package [28]. In this process, (i) solvent and non-complexed ions were removed; and (ii) hydrogens and charges (of amber ff99sb force field) were added to the protein [28]. The docking was performed with DOCK 6.5 program (UCSF) under the following parameters: 5AAB, box margin = 0.5 nm, select spheres = 1.0 nm, max orientations = 2300; 5AA8, box margin = 0.5 nm, select spheres = 1.0 nm, max orientations = 5000.
AutoDock Vina analyses were performed in cubes with 1.5 nm side length (5AAB, center_x = 37.935, center_y = 47.567, center_z = 17.016; 5AA8, center_x = 38.053, center_y = 47.101, center_z = 16.983). The side chains of the ligands were allowed for flexible torsion. The maximum energy difference between the optimal energy and the highest energy was set to be 3, and the lowest predicted energies were extracted for further analysis.
The ALK-inhibitor complexes were first analyzed by grid scoring of DOCK 6.5, followed by a second round scoring with Descriptor and Hawkins GB/SA algorithms. Descriptor scoring calculates the standard energy within a system. Hawkins GB/SA scoring calculates the energy through molecular mechanics generalized born surface area (MM/GBSA).

Molecular Dynamics Simulation
Molecular dynamics simulations were performed using GROMACS 4.6.7 [29] package and Amber ff99sb force field with TIP3P water model [30]. Crizotinib and lorlatinib were under the general Amber force field (GAFF) and the charges were added by AM1-BBC method. We performed system check with parmchk program and generated additional parameters of force field. To generate the ligand topology file, we used AnteChamber PYthon Parser interfacE (ACPYPE) [31]. Particle Mesh Ewald (PME) [32] was utilized to consider the long-range electrostatic interactions and the Linear Constraint Solver (LINCS) [33] algorithm was used to constrain bonds. The receptor-ligand complexes were solvated in a dodecahedron box of water, with a distance of 1.0 between the solute and the box. All systems were neutralized by adding Na + and Cl − at 0.15 mol/L. Before MD simulations, the complexes were relaxed to <1000 kJ/mol·nm by up to 50,000 cycles of steep descent minimization. After energy minimization, temperature of the system was controlled in the constant number of particles, volume, and temperature (NVT) ensemble to 300 K over 100 ps. The 100 ps constant number of particles, pressure, and temperature (NPT) equilibration was then performed with a reference pressure of 1 bar. After that, 30 ns MD simulations were performed with a time step of 2 fs and the coordinates of the complexes were saved every 8 ps.

Root-Mean-Square Fluctuation and Root-Mean-Square Deviation Analysis
The RMSF, RMSD, and secondary structure elements were analyzed by g_rmsf, g_rmsd, and do_dssp modules of the GROMACS 4.6.7 suite [29].

Binding Energy Calculation and Energy Decomposition
Molecular dynamics trajectory visualization was performed with VMD 1.9.2 [34]. The g_mmpbsa module of GROMACS 4.6.7 and mm_pbsa.pl tool of Amber 9 (University of California, San Francisco, CA, USA) were applied to calculate the free energies in biomolecular interactions.
The MM/PBSA was used to calculate the binding free energy of ALK-ligand complexes, as described in our previous study [35]. In this approach, the binding free energy (∆G) calculation formula is [36]: The Van der Waals (∆E vdw ) and electrostatic interaction (∆E ele ) can be calculated through molecular mechanics (MM) method (Equation (2)). The bonded interactions (∆E bonded ) consisted of bond, angle, dihedral, and improper interactions. In the single trajectory approach, ∆E bonded is always taken as zero because the conformations of protein and ligand are identified as a constant before and after they bond together [37]. The free energy of salvation (∆G sol ) consists of the polar solvation free energy and the nonpolar part which can be estimated with the Poisson-Boltzmann (PB) equation and the solvent-accessible surface area (SASA). −T∆S stands for the changing of conformational entropy upon ligand binding, which usually is ignored in practice because of its high computational cost and low prediction accuracy.
In our research, the binding free energy by MM/GBSA is the sum of the gas-phase interaction energy (∆E MM ) plus the solvation free energy (∆G sol ) minus the product of the absolute temperature and entropy change (T∆S) (Equation (4)).
The gas-phase interaction energy (∆E MM ) was composed of the van der Waals (∆E vdw ) and electrostatic interaction (∆E ele ) energies between receptor and ligand, and the solvation free energy (∆G sol ) was the sum of the electrostatic contributions calculated by the generalized Born (GB) approximation model (E GB ) and the nonpolar contributions obtained using the SASA (E SUR ). The entropic contribution (T∆S) was ignored as it is computationally expensive to calculate and yields a low prediction accuracy. These terms were calculated using Equations (5) and (6).

Calculation of the Electrostatic Energy and Root-Mean-Square Fluctuation Change Trends
If the energy change from ALK-crizotinib to ALK-lorlatinib of the same amino acid residue was the same between 5-10 ns and 20-25 ns simulation, the residue received a score of 1. If the energy change trend was opposite, the residue received a score of −1. The trend between 20-25 ns and 25-30 ns was scored in the same way. The total energy change trend parameter was calculated by adding the score of the two-time periods together. A value of "2" suggested that the energy was changing in the same trend from ALK-crizotinib to ALK-lorlatinib in the investigated periods of time. The RMSF change trend of the key amino acids was analyzed in the same way as described above.

Conclusions
Acquired resistance to small molecule inhibitors is now a great challenge in the treatment of NSCLC and many other types of cancer. It is of great value to understand the molecular mechanism resulting in the occurrence of such resistance, which will in turn help the development of more specific and potent drugs to treat the disease. In this study, we exploited molecular dynamics simulation approach to understand how the lorlatinib ALK resistance mutation L1198F led to the resensitization to the first-generation inhibitor crizotinib. We found that crizotinib bound to C1156Y-L1198F ALK with a higher affinity than lorlatinib. This difference is at least partially caused by the conformational change that resulted from the L1198F mutation, as revealed by the RMSD analysis. By energy decomposition, we identified five amino acid residues that are located in the inhibitor and ATP binding domains of ALK. Root-mean-square fluctuation analysis further confirmed the importance of His1124 and Glu1210 as key amino acid residues that regulate the activity of ALK. We also found that Tyr1278 and Arg1279, located in the activation loop of the kinase domain, were also affected by the L1198 mutation. With all these results, we concluded that the L1198F mutation led to conformational changes at the inhibitor and ATP binding sites of ALK protein, which affected the downstream phosphorylation of the key residues. In clinical practice, it is very rare to see a resistance point mutation result in resensitization to the previous generation of inhibitors. In this study, we elucidated the molecular mechanism of such an intriguing drug sensitivity shift. These findings are valuable to the design of new targeted therapies for the treatment of ALK-positive cancer.