Deciphering Conformational Changes of the GDP-Bound NRAS Induced by Mutations G13D, Q61R, and C118S through Gaussian Accelerated Molecular Dynamic Simulations

The conformational changes in switch domains significantly affect the activity of NRAS. Gaussian-accelerated molecular dynamics (GaMD) simulations of three separate replicas were performed to decipher the effects of G13D, Q16R, and C118S on the conformational transformation of the GDP-bound NRAS. The analyses of root-mean-square fluctuations and dynamics cross-correlation maps indicated that the structural flexibility and motion modes of the switch domains involved in the binding of NRAS to effectors are highly altered by the G13D, Q61R, and C118Smutations. The free energy landscapes (FELs) suggested that mutations induce more energetic states in NRAS than the GDP-bound WT NRAS and lead to high disorder in the switch domains. The FELs also indicated that the different numbers of sodium ions entering the GDP binding regions compensate for the changes in electrostatic environments caused by mutations, especially for G13D. The GDP–residue interactions revealed that the disorder in the switch domains was attributable to the unstable hydrogen bonds between GDP and two residues, V29 and D30. This work is expected to provide information on the energetic basis and dynamics of conformational changes in switch domains that can aid in deeply understanding the target roles of NRAS in anticancer treatment.


Introduction
Small GTPases proteins (KRAS, HRAS, and NRAS) coded by human RAS genes regulate the levels of numerous cell signaling processes implicated in cell growth [1][2][3]. KRAS, HRAS, and NRAS have in common the feature of exchanging between the GTPbound active state and the GDP-bound inactive state [4][5][6]. According to previous reports, the hydrolysis reaction of guanosine triphosphate (GTP) into guanosine diphosphate (GDP) can be accelerated by the GTPase activating proteins (GAPs), which generate an inactive form of RAS protein [7]. In contrast, the transformation of GDP into GTP can be catalyzed by guanosine exchange factors, leading to an active state for the GTP-associated RAS protein [7,8]. Significant signaling pathways, responsible for key roles in cell proliferation and survival, are triggered by the binding of the GTP-bound RAS proteins to effectors, mainly involving RAF, phosphoinositide 3-kinase (PI3K), and Ral guanine nucleotide dissociation stimulator (RalGDS) [9]. Residue mutation-mediated active RAS proteins generally drive human cancers and tumorigenesis by hyperactivating the downstream signal pathways [10][11][12]. Therefore, exploring the effects of mutations on the activity of RAS proteins would be of great significance to for a deeper understanding of the target role of RAS proteins in human cancers.
Among human cancers induced by mutations of RAS proteins, KRAS is the most frequently mutated oncogene, mainly occurring at the residues G12 and G13 and accounting allosteric position and can affect the activity of NRAS. To achieve our aim, multiple-replica GaMD (MR-GaMD) simulations, free energy landscapes (FELs), principal component analysis (PCA) [70][71][72][73], and dynamics cross-correlation maps (DCCMs) [71,74] were adopted in the current study. This work is expected to contribute useful information on dynamics and energy basis to aid the understanding of the target role of NRAS in drug design aimed at anticancer treatment.

Structural Flexibility and Motion Modes between Structure Domainsof NRAS
To investigate the effects of G13D, Q61R, and C118S on the local structural flexibility of NRAS, the difference in the RMSF for each Cα atom between the mutated and the wild-type NRAS was estimated through the equation, ∆ = − (Figure 2A). It was found that the most obvious effect of mutations on the structural flexibility occurs at the two switch domains SWI and SWII. The mutations of G13D and Q61R weaken the structural flexibility of the SWI, but C118S strengthens it ( Figures 1A and 2A).In contrast the weak impact of C118S on the structural flexibility of the SWII, G13D greatly enhances its structural flexibility; however, Q61R strongly reduces the structural flexibility of this switch domain ( Figures 1A and 2A). It was also observed that G13D and C118S slightly strengthen the structural flexibility of loop L1. Furthermore, G13D increases the structural flexibility of loops L3 and L4, while Q61R slightly weakens that of loop L3 (Figures 1A and 2A).

Structural Flexibility and Motion Modes between Structure Domains of NRAS
To investigate the effects of G13D, Q61R, and C118S on the local structural flexibility of NRAS, the difference in the RMSF for each C α atom between the mutated and the wild-type NRAS was estimated through the equation, ∆RMSF = RMSF mutant − RMSF WT (Figure 2A). It was found that the most obvious effect of mutations on the structural flexibility occurs at the two switch domains SWI and SWII. The mutations of G13D and Q61R weaken the structural flexibility of the SWI, but C118S strengthens it ( Figures 1A and 2A). In contrast the weak impact of C118S on the structural flexibility of the SWII, G13D greatly enhances its structural flexibility; however, Q61R strongly reduces the structural flexibility of this switch domain ( Figures 1A and 2A). It was also observed that G13D and C118S slightly strengthen the structural flexibility of loop L1. Furthermore, G13D increases the structural flexibility of loops L3 and L4, while Q61R slightly weakens that of loop L3 ( Figures 1A and 2A). To reveal mutation-mediated influences on the global flexibility of NRAS, molecular surface areas (MSAs) of the GDP-bound WT, G13D, Q61R, and C118S NRAS were calculated by using the single-joined MR-GaMD trajectory (SJMT), and their frequency distributions are depicted in Figure 2B. The MSAs of the GDP-bound G13D, Q61R, and C118S NRAS were decreased by 242.5, 238.5, and 235.5 Å 2 , respectively, which indicates a lower global structural flexibility for NRAS with the binding state caused by the mutations, G13D, Q61R, and C118S.
To probe mutation-induced effects on motion modes between structural domains of NRAS, DCCMs of the Cα atoms from the GDP-bound WT, G13D, Q61R, and C118S NRAS were calculated based on the SJMT ( Figure 2C-F). It can be noted that G13D, Q61R, and C118S exert evident effects on correlated motions between structural domains of NRAS. In the GDP-bound WT NRAS ( Figure 2C), the domains D1 and D2 To reveal mutation-mediated influences on the global flexibility of NRAS, molecular surface areas (MSAs) of the GDP-bound WT, G13D, Q61R, and C118S NRAS were calculated by using the single-joined MR-GaMD trajectory (SJMT), and their frequency distributions are depicted in Figure 2B. The MSAs of the GDP-bound G13D, Q61R, and C118S NRAS were decreased by 242.5, 238.5, and 235.5 Å 2 , respectively, which indicates a lower global structural flexibility for NRAS with the binding state caused by the mutations, G13D, Q61R, and C118S.
To probe mutation-induced effects on motion modes between structural domains of NRAS, DCCMs of the C α atoms from the GDP-bound WT, G13D, Q61R, and C118S NRAS were calculated based on the SJMT ( Figure 2C-F). It can be noted that G13D, Q61R, and C118S exert evident effects on correlated motions between structural domains of NRAS. In the GDP-bound WT NRAS ( Figure 2C), the domains D1 and D2 produce slightly anti-correlated motions (cyan), while the domain D3 generates obvious anti-correlated movement (blue). Among the D1-D3 domains from the GDP-bound WT NRAS, the D1 and D2 domains, respectively, exhibit slightly anti-correlated movements for the switches SW I and SW II relative to the P-loop (residues [8][9][10][11][12][13][14][15][16], while the D3 domain shows strong anticorrelated motion for the switch SW II relative to the SW. The domains D4 and D5 from the GDP-bound WT NRAS yield strongly positive correlated motions, indicated by the yellow and red in Figure 2C, in which the D4 shows strongly positive correlated motion for the β-sheet β2 relative to β1 ( Figures 1A and 2C),while the D5 domain shows strongly positive correlated movement between the loop L4 and the P-loop ( Figures 1A and 2C). Compared to the GDP-bound WT NRAS, G13D obviously strengthens the anti-correlated motions of the SWI and SWII relative to the P-loop but slightly weakens the anti-correlated movement between the SWI and SWII ( Figure 2D);moreover, G13D strengthens the positive correlated motions of β2 relative to β1 and ofL4 relative to the P-loop ( Figure 2D). In comparison with the GDP-bound WT NRAS, Q61R weakens not only the anti-correlated motion between the SWI and the SW II but also the positive correlated motions of β2 relative to β1 and ofL4 relative to the P-loop ( Figure 2E). By interacting with the GDP-bound WT NRAS, C118S slightly strengthens the anti-correlated movement between the SWI and the SW II but softly weakens the positive correlated motion of L4 relative to the P-loop ( Figure 2F).
It was found from the above analyses that G13D, Q61R, and C118S change the structural flexibility and correlated motions of the switch domains. Structurally, the switch domains are involved in the binding of the RAS proteins to effectors. Therefore, the alterations in the conformation and internal dynamics of the switch domains certainly produce a significant effect on the binding of NRAS to effectors, which regulates the activity of NRAS. The work of Kessler et al. suggests that Q61R and C118S change the "on" or "off" state of NRAS and affect the binding modes of inhibitors to NRAS [30], which agrees with our current findings. The study from Johnson et al. verifies that G13D greatly changes the structural flexibility of the switch domains and exerts vital impacts on the activity of NRAS [31], supporting our results well.

Conformational Changes of NRAS Revealed by Free Energy Landscapes
To reveal the energetic basis of mutation-mediated conformational changes of NRAS, FELs were constructed using the distance of residue D33 from residue 61 (Q61 or R61) and root-mean-square deviations (RMSDs) of backbone atoms from NRAS as reaction coordinates. The reasons why we selected the seas reaction coordinates are as follows: (1) the distance from D33 to residue 61 efficiently reflects the transformation of the conformational state between SW I and SW II, (2) the RMSDs of backbone atoms can rationally embody the total structural fluctuations of NRAS through MR-GaMD simulations. To understand the changes in the interaction sites of GDP with NRAS, a protein−ligand interaction profiler (PLIP) server [75,76] was applied to detect interaction networks of GDP with NRAS in different energetic states.
For the GDP-bound WT NRAS, MR-GaMD simulations detected two main energy valleys EV1 and EV2 ( Figure 3A), indicating that the GDP-bound WT NRAS is distributed at two conformational sub-spaces. The depth of the potential well EV1 was much deeper than that of the EV2, implying that the probability of transition from EV2 into EV1 of NRAS is much higher than that from EV1 into EV2. Thus, the conformations of the GDP-bound WT NRAS mostly fall into EV1. The distances of D33 from Q62 in the energetic states EV1 and EV2 were25.4 and 12.4 Å, respectively, showing that the switch domains SWI and SWII of the GDP-bound WT NRAS situated at EV1 have an open conformation, while those located at EV2 exist in a closed conformation ( Figure 3B,C). In two energetic states, GDP forms a salt bridge interaction with K16 and D119 while GDP produces hydrogen bonding interactions (HBIs) with the common residues G15, K16, S17, A18, V29, D30, N116, K117, and K147 ( Figure S1A,B, see Supplementary Materials). Although GDP loses two HBIs with V14 and S146 in EV2 compared to EV1, two new HBIs also appear between GDP and the two residues G13 and A145 in EV2 ( Figure S1A,B). The representative structure at EV1 was superimposed together with that at the EV2 ( Figure 3D). It was found that the switch domains SWI and II become highly disordered by interacting with the other structural domains. In addition, the helix α3 and the loops L1 and L3 generate slight deviations. Moreover, the structures of the GDP and Mg 2+ ion situated at EV1 and EV2 are aligned together ( Figure S1C). The adenine group and the middle ring from GDP obviously slide away from each other and the Mg 2+ ions in EV1 and EV2 generate a deviation of 2.2 Å ( Figure 1C), which can strongly affect the binding of NRAS to effectors. HBIs also appear between GDP and the two residues G13 and A145 in EV2 ( Figure  S1A,B). The representative structure at EV1 was superimposed together with that at the EV2 ( Figure 3D). It was found that the switch domains SWI and II become highly disordered by interacting with the other structural domains. In addition, the helix α3 and the loops L1 and L3 generate slight deviations. Moreover, the structures of the GDP and Mg 2+ ion situated at EV1 and EV2 are aligned together ( Figure S1C). The adenine group and the middle ring from GDP obviously slide away from each other and the Mg 2+ ions in EV1 and EV2 generate a deviation of 2.2 Å ( Figure 1C), which can strongly affect the binding of NRAS to effectors. With regard to the GDP-bound G13D NRAS, three low energetic valleys EV1-EV3 were captured by MR-GaMD simulations ( Figure 4A), showing that the GDP-bound G13D NRAS is mainly clustered into three conformational sub-spaces. The potential well EV1 was found to be deeper than EV2 and EV3; as a result, the GDP-bound G13D NRAS was mostly located at the energetic valley EV1 ( Figure 4A). The distances of D33 from Q61 were15.0, 18.1, and 16.4 Å at EV1, EV2, and EV3, respectively; hence, the switch domains SW I and SWII form an open state at EV2, while the switch domains SWI and II are located at an intermediate state at EV1 and EV3 ( Figure 4B-D). Similarly to the GDP-bound WT NRAS, GDP forms salt bridge interactions with K16 and D119 in three energetic states EV1, EV2, and EV3 ( Figure S2). Meanwhile, GDP yields HBIs with the conserved residues G15, K16, S17, A18, N116, K117, S145, A146, and K147 ( Figure S2). However, G13D also induces an alteration in interaction sites of GDP with NRAS. In detail, it was found that two new HBIs appear between GDP and the residues D13 and V29 in EV1 ( Figure S2A), GDP generates two additional HBIs with V14 and D33 in EV2 (Fig- With regard to the GDP-bound G13D NRAS, three low energetic valleys EV1-EV3 were captured by MR-GaMD simulations ( Figure 4A), showing that the GDP-bound G13D NRAS is mainly clustered into three conformational sub-spaces. The potential well EV1 was found to be deeper than EV2 and EV3; as a result, the GDP-bound G13D NRAS was mostly located at the energetic valley EV1 ( Figure 4A). The distances of D33 from Q61 were15.0, 18.1, and 16.4 Å at EV1, EV2, and EV3, respectively; hence, the switch domains SW I and SWII form an open state at EV2, while the switch domains SWI and II are located at an intermediate state at EV1 and EV3 ( Figure 4B-D). Similarly to the GDP-bound WT NRAS, GDP forms salt bridge interactions with K16 and D119 in three energetic states EV1, EV2, and EV3 ( Figure S2). Meanwhile, GDP yields HBIs with the conserved residues G15, K16, S17, A18, N116, K117, S145, A146, and K147 ( Figure S2). However, G13D also induces an alteration in interaction sites of GDP with NRAS. In detail, it was found that two new HBIs appear between GDP and the residues D13 and V29 in EV1 ( Figure S2A), GDP generates two additional HBIs with V14 and D33 in EV2 ( Figure S2B), and GDP not only forms three additional HBIs with D13, V14, and D30 but also produces a new salt bridge interaction with D30 in EV3 ( Figure S2C). The alignment of the structures situated at EV1, EV2, and EV3 was carried out ( Figure 4E), and the results showed that the switch domains are in a highly disordered state. It was also found that two additional Na + ions appear at the EV1 and a Na + ion is at EV3 ( Figure 4B,D,E). The structures of GDP and Mg 2+ ions in EV1-EV3 were also superimposed together ( Figure 4F). The results indicated that the adenine group, the middle ring, and the phosphate group of GDP are aligned well. Although G13D induces more energetic states compared to the GDP-bound WT NRAS, the appearance of Na + ions caused by G13D maintains the structural stability of GDP.
Molecules 2022, 27, x FOR PEER REVIEW 7 of 20 ure S2B), and GDP not only forms three additional HBIs with D13, V14, and D30 but also produces a new salt bridge interaction with D30 in EV3 ( Figure S2C). The alignment of the structures situated at EV1, EV2, and EV3 was carried out ( Figure 4E), and the results showed that the switch domains are in a highly disordered state. It was also found that two additional Na + ions appear at the EV1 and a Na + ion is at EV3 ( Figure 4B,D,E). The structures of GDP and Mg 2+ ions in EV1-EV3 were also superimposed together (Figure 4F). The results indicated that the adenine group, the middle ring, and the phosphate group of GDP are aligned well. Although G13D induces more energetic states compared to the GDP-bound WT NRAS, the appearance of Na + ions caused by G13D maintains the structural stability of GDP. In the case of the GDP-bound Q61R NRAS, four energetic valleys EV1-EV4 were detected in the MR-GaMD simulations ( Figure 5A), suggesting that the GDP-bound Q61R NRAS is mainly classed into four conformational sub-spaces. The depths of the four potential wells EV1-EV4 were almost the same; hence, the conformations of the GDP-bound Q61R NRAS were almost equally distributed at four conformational sub-spaces. The distances of D33 from R61 were 12.6, 16.2, 19.4, and 23.3 Å in the energetic states EV1, EV2, EV3, and EV4, respectively( Figure 5B-E). Therefore, it was found that the switch domains SW I and SWII form a closed state in EV1, induce an intermediate state in EV2, and lead to two open states in EV3 and EV4. In the four energetic states EV1-EV4, GDP produces two salt bridge interactions with K16 and D119 ( Figure S3). It was also observed that GDP forms HBIs with the conserved residues V14, G15, K16, S17, N116, K117, S145, and K147 in all four energetic states ( Figure S3). Interestingly, a π-π interaction appears between F28 and the adenine group of GDP in EV1 and EV2 ( Figure  S3A,B). In EV1, EV3, and EV4, GDP generates an additional HBI with A146 ( Figure  S3A,C,D). The residue D30 forms an HBI with the phosphate group of GDP in EV1 and EV4 ( Figure S3A,D). In addition, GDP also yields an HBI with V29 from the SWI in the EV1. The representative structures situated at EV1-EV4 were superimposed together ( Figure 5F). As exhibited in Figure 5F, the switch SWI has a highly disordered state while the SWII only generates obvious deviations. In addition, the helix α3 and the loops In the case of the GDP-bound Q61R NRAS, four energetic valleys EV1-EV4 were detected in the MR-GaMD simulations ( Figure 5A), suggesting that the GDP-bound Q61R NRAS is mainly classed into four conformational sub-spaces. The depths of the four potential wells EV1-EV4 were almost the same; hence, the conformations of the GDPbound Q61R NRAS were almost equally distributed at four conformational sub-spaces. The distances of D33 from R61 were 12.6, 16.2, 19.4, and 23.3 Å in the energetic states EV1, EV2, EV3, and EV4, respectively( Figure 5B-E). Therefore, it was found that the switch domains SW I and SWII form a closed state in EV1, induce an intermediate state in EV2, and lead to two open states in EV3 and EV4. In the four energetic states EV1-EV4, GDP produces two salt bridge interactions with K16 and D119 ( Figure S3). It was also observed that GDP forms HBIs with the conserved residues V14, G15, K16, S17, N116, K117, S145, and K147 in all four energetic states ( Figure S3). Interestingly, a π-π interaction appears between F28 and the adenine group of GDP in EV1 and EV2 ( Figure S3A,B). In EV1, EV3, and EV4, GDP generates an additional HBI with A146 ( Figure S3A,C,D). The residue D30 forms an HBI with the phosphate group of GDP in EV1 and EV4 ( Figure S3A,D). In addition, GDP also yields an HBI with V29 from the SWI in the EV1. The representative structures situated at EV1-EV4 were superimposed together ( Figure 5F). As exhibited in Figure 5F, the switch SWI has a highly disordered state while the SWII only generates obvious deviations. In addition, the helix α3 and the loops L1 and L4 also produce evident deviations. It was also observed that additional sodium ions Na + respectively appear at EV2, EV3, and EV4 ( Figure 5C-E). According to the structural alignment of GDP and Mg 2+ ions in the  Figure S4), the structures of GDP were aligned well in four energetic states, and they showed high stability throughout the entirety of the MR-GaMD simulations but, due to repulsive interactions of additional Na + and Mg 2+ ions,they yielded evident deviations. L1 and L4 also produce evident deviations. It was also observed that additional sodium ions Na + respectively appear at EV2, EV3, and EV4 ( Figure 5C-E). According to the structural alignment of GDP and Mg 2+ ions in the EV1-EV4 ( Figure S4), the structures of GDP were aligned well in four energetic states, and they showed high stability throughout the entirety of the MR-GaMD simulations but, due to repulsive interactions of additional Na + and Mg 2+ ions,they yielded evident deviations. As for the GDP-bound C118S NRAS, MR-GaMD simulations captured three low energetic valleys EV1-EV3 ( Figure 6A), verifying that the GDP-bound C118S NRAS is distributed at three conformational subspaces. The potential well EV2 was found to be deeper than that of EV1 and EV3; hence, the conformations of the GDP-bound C118S NRAS were mostly situated at EV2. The distances of D33 from Q61 were 12.6, 16.2, and 19.4 Å in the energy valleys EV1, EV2, and EV3 ( Figure S5), respectively, suggesting that the switch domains SW I and SWII individually form a closed conformation in the EV1, an intermediate state in the EV2, and an open conformation in the EV3. According to Figure 6B-D, GDP not only produces salt bridge interactions with K16 and D119 but also forms HBIs with the common residues G13, V14, G15, K16, D30, N116, K117, S145, and K147. In addition, a new HBI is formed between GDP and A146 in EV1 ( Figure 6B), three additional HBIs of GDP with S17, A18, and E31 are generated in EV2 ( Figure 6C), and two new HBIs of GDP with S17 and A18 appear in EV3 ( Figure 6D), which reflect the effect of conformational changes on interaction sites of GDP with NRAS. The representative structures in EV1-EV3 were superimposed together ( Figure S6A) and the results indicated that the switch domains SW I and SWII are highly disordered. In addition, the P-loop, the helix α3, and the loop L1 were found to generate slight deviations ( Figure  S6A). Interestingly, two additional Na + ions were detected in the EV2, which possibly affect the electrostatic environment around GDP (Figures S5B and S6A). The structures of GDP and Mg 2+ ions were aligned together ( Figure S6B). Although the adenine group and the middle ring of GDP were aligned well, the phosphate group of GDP and Mg 2+ ions yielded evident deviations ( Figure S6B). As for the GDP-bound C118S NRAS, MR-GaMD simulations captured three low energetic valleys EV1-EV3 ( Figure 6A), verifying that the GDP-bound C118S NRAS is distributed at three conformational subspaces. The potential well EV2 was found to be deeper than that of EV1 and EV3; hence, the conformations of the GDP-bound C118S NRAS were mostly situated at EV2. The distances of D33 from Q61 were 12.6, 16.2, and 19.4 Å in the energy valleys EV1, EV2, and EV3 ( Figure S5), respectively, suggesting that the switch domains SW I and SWII individually form a closed conformation in the EV1, an intermediate state in the EV2, and an open conformation in the EV3. According to Figure 6B-D, GDP not only produces salt bridge interactions with K16 and D119 but also forms HBIs with the common residues G13, V14, G15, K16, D30, N116, K117, S145, and K147. In addition, a new HBI is formed between GDP and A146 in EV1 ( Figure 6B), three additional HBIs of GDP with S17, A18, and E31 are generated in EV2 ( Figure 6C), and two new HBIs of GDP with S17 and A18 appear in EV3 ( Figure 6D), which reflect the effect of conformational changes on interaction sites of GDP with NRAS. The representative structures in EV1-EV3 were superimposed together ( Figure S6A) and the results indicated that the switch domains SW I and SWII are highly disordered. In addition, the P-loop, the helix α3, and the loop L1 were found to generate slight deviations ( Figure S6A). Interestingly, two additional Na + ions were detected in the EV2, which possibly affect the electrostatic environment around GDP (Figures S5B and S6A). The structures of GDP and Mg 2+ ions were aligned together ( Figure S6B). Although the adenine group and the middle ring of GDP were aligned well, the phosphate group of GDP and Mg 2+ ions yielded evident deviations ( Figure S6B). In summary, the analyses of the FELs revealed that G13D, Q61R, and C118S induce more energetic states than the GDP-bound WT NRAS, which increases the extent of disorder in the switch domains in the mutated states of NRAS. The work of Johnson et al. has verified that G13D leads to more open SWI and II and more disordered switch domains for NRAS due to the increased entropic contributions of conformations caused by G13D [31]. The representative conformations captured by the FELs indicated that G13D, Q61R, and C118S cause different numbers of Na + ions to appear at some energetic states and compensate for the changes in the electrostatic environment induced by mutations, which certainly affect the interactions of GDP with NRAS. This interesting phenomenon was also found in previous work [77,78],supporting our current study well.

Principal Component Analysis
To better understand the effects of G13D, Q61R, and C118S on the conformational changes of the switch domains caused by NRAS, PCA was carried out using the SJMT. In general, eigenvalues arising from PCA are applied to clarify the total motion intensity of receptors. In the current study, the functions of eigenvalues as eigenvector indexes were used, as shown in Figure S7. The first several, larger eigenvalues mainly represent concerted motions of structural domains in NRAS. It can be seen from Figure S7 that the first five eigenvalues account for 89.9, 92.4, 83.3, and 95.2% of the total movements of the GDP-bound WT, G13D, Q61R, and C118S NRAS, respectively. In comparison with the WT NRAS, the first eigenvalues of the GDP-bound G13D and C118S NRAS are increased, while those of the GDP-bound Q61R NRAS are greatly reduced. Therefore, G13D and C118S strengthen the total motion intensity of NRAS relative to the WT state but Q61R weakens it. In summary, the analyses of the FELs revealed that G13D, Q61R, and C118S induce more energetic states than the GDP-bound WT NRAS, which increases the extent of disorder in the switch domains in the mutated states of NRAS. The work of Johnson et al. has verified that G13D leads to more open SWI and II and more disordered switch domains for NRAS due to the increased entropic contributions of conformations caused by G13D [31]. The representative conformations captured by the FELs indicated that G13D, Q61R, and C118S cause different numbers of Na + ions to appear at some energetic states and compensate for the changes in the electrostatic environment induced by mutations, which certainly affect the interactions of GDP with NRAS. This interesting phenomenon was also found in previous work [77,78],supporting our current study well.

Principal Component Analysis
To better understand the effects of G13D, Q61R, and C118S on the conformational changes of the switch domains caused by NRAS, PCA was carried out using the SJMT. In general, eigenvalues arising from PCA are applied to clarify the total motion intensity of receptors. In the current study, the functions of eigenvalues as eigenvector indexes were used, as shown in Figure S7. The first several, larger eigenvalues mainly represent concerted motions of structural domains in NRAS. It can be seen from Figure S7 that the first five eigenvalues account for 89.9, 92.4, 83.3, and 95.2% of the total movements of the GDP-bound WT, G13D, Q61R, and C118S NRAS, respectively. In comparison with the WT NRAS, the first eigenvalues of the GDP-bound G13D and C118S NRAS are increased, while those of the GDP-bound Q61R NRAS are greatly reduced. Therefore, G13D and C118S strengthen the total motion intensity of NRAS relative to the WT state but Q61R weakens it.
The first eigenvector corresponding to the first eigenvalue can reflect the main motion behavior of NRAS. Based on this, the first eigenvector stemming from the PCA was visualized using the software VMD [79], along with the optimized structure (Figure 7). The switch domains SW I and SW II have stronger motions compared to the other structural domains of NRAS. With regard to the GDP-bound WT NRAS, the two switches SWI and II move in an opposite direction away from each other, inducing a loose switch state; furthermore, the motion strength of SW I is stronger than that of the SW II ( Figure 7A). By interacting with the WT NRAS, G13D changes the motion direction of SW I and SW II, and these two switches move close to each other ( Figure 7B). Moreover, G13D not only alters the motion direction of the P-loop, the helix α3, and the loop L1 but also strengthens the movement of these three structural domains ( Figure 7B). In contrast to the WT NRAS, Q61R changes the motion direction of the switch domains and inhibits their motion strength; at the same time, Q61R also strengthens the motion of loop L1 ( Figure 7C). Compared to the WT NRAS, C118S slightly inhibits the movement of SW II and the partial regions in SW I, but it enhances the moving intensity of loop L1 and helix α3 ( Figure 7D). The first eigenvector corresponding to the first eigenvalue can reflect the main motion behavior of NRAS. Based on this, the first eigenvector stemming from the PCA was visualized using the software VMD [79], along with the optimized structure (Figure 7). The switch domains SW I and SW II have stronger motions compared to the other structural domains of NRAS. With regard to the GDP-bound WT NRAS, the two switches SWI and II move in an opposite direction away from each other, inducing a loose switch state; furthermore, the motion strength of SW I is stronger than that of the SW II ( Figure 7A). By interacting with the WT NRAS, G13D changes the motion direction of SW I and SW II, and these two switches move close to each other ( Figure 7B). Moreover, G13D not only alters the motion direction of the P-loop, the helix α3, and the loop L1 but also strengthens the movement of these three structural domains ( Figure  7B). In contrast to the WT NRAS, Q61R changes the motion direction of the switch domains and inhibits their motion strength; at the same time, Q61R also strengthens the motion of loop L1 ( Figure 7C). Compared to the WT NRAS, C118S slightly inhibits the movement of SW II and the partial regions in SW I, but it enhances the moving intensity of loop L1 and helix α3 ( Figure 7D). According to the above analyses, the three mutations G13D, Q61R, and C118S from NRAS generate different effects on the motion behavior of the switch domains. It is well-known that SW I is involved in the binding of NRAS to effectors; hence, G13D, Q61R, and C118S tune the activity of NRAS. Except for the obvious effect of G13D on the conformation of the P-loop, G13D, Q61R, and C118S also produce impacts on the movements of the helix α3 and the loop L1;structurally, these two regions are close to According to the above analyses, the three mutations G13D, Q61R, and C118S from NRAS generate different effects on the motion behavior of the switch domains. It is wellknown that SW I is involved in the binding of NRAS to effectors; hence, G13D, Q61R, and C118S tune the activity of NRAS. Except for the obvious effect of G13D on the conformation of the P-loop, G13D, Q61R, and C118S also produce impacts on the movements of the helix α3 and the loop L1;structurally, these two regions are close to the allosteric position of NRAS to some extent. Thus, G13D, Q61R, and C118S possibly affect the allosteric regulation of the activity of NRAS. The conformational changes in the switch domain induced by the phosphorylation of Y32 disturb the binding of the RAS protein to effectors, which agrees with our results [80].

Interaction Networks of GDP and NRAS
The previous analysis of FELs indicates that conformational changes caused by G13D, Q61R, and C118S induce the alterations in the interaction sites for GDP and NRAS. To evaluate the stability of the GDP-residue interactions, the distances involved in the π-π interaction and the salt bridge interactions were calculated, and their frequency distributions are depicted in Figure 8. The occupancies of HBIs used to describe the stability of HBIs were analyzed using the program CPPTRAJ in Amber 20. The information on HBIs is listed in Table 1. the allosteric position of NRAS to some extent. Thus, G13D, Q61R, and C118S possibly affect the allosteric regulation of the activity of NRAS. The conformational changes in the switch domain induced by the phosphorylation of Y32 disturb the binding of the RAS protein to effectors, which agrees with our results [80].

Interaction Networks of GDP andNRAS
The previous analysis of FELs indicates that conformational changes caused by G13D, Q61R, and C118S induce the alterations in the interaction sites for GDP and NRAS. To evaluate the stability of the GDP-residue interactions, the distances involved in the π-π interaction and the salt bridge interactions were calculated, and their frequency distributions are depicted in Figure 8. The occupancies of HBIs used to describe the stability of HBIs were analyzed using the program CPPTRAJ in Amber 20. The information on HBIs is listed in Table 1.   The structure of the lowest energy taken from the SJMT of the GDP-bound WT NRAS was utilized to depict the geometric position of the π-π interaction and the salt bridge interactions ( Figure 8A). The distance of the salt bridge interaction between K16 and the phosphate group of GDP was calculated, and its frequency distribution is exhibited in Figure 8B. The distance of K16 from the phosphate group was distributed at 3.6, 3.4, 3.6, and 3.4 Å in the GDP-bound WT, G13D, Q61R, and C118S NRAS, respectively, suggesting that G13D and C118S strengthen the salt bridge interaction of K16 with GDP compared to the GDP-bound WT NRAS. The phenyl group of GDP generates the π-π interaction with the adenine group of GDP ( Figure 8A), and the frequency distribution of its distance is displayed in Figure 8C. The peak value of the distance between the phenyl of F28 and the adenine group of GDP was located at 4.7 Å in the four current systems, and only the distribution width of the distance in the GDP-bound WT NRAS was wider than that in the three other systems; thus, G13D, Q61R, and C118S slightly strengthened the π-π interaction of F28 with GDP relative to the GDP-bound WT NRAS. The carbonyl of D119 produces the salt bridge interaction with the adenine group of GDP ( Figure 8A), and the frequency distribution of this salt bridge distance is shown in Figure 8D. The peak value of the frequency distribution for this salt bridge was situated at 6.1 Å in the four current systems, implying that G13D, Q61R, and C118S hardly affect the strength of the salt bridge between D119 and the adenine group of GDP.
The structure of the lowest energy taken from the SJMT of the GDP-bound WT NRAS was also employed to describe the geometric position of HBIs (Figure 9). According to Table 1, except for V14, the residues G13(D13), G15, K16, S17, and A18, located at the P-loop, formed HBIs with the phosphate group of GDP in the four simulated systems ( Figure 9A), and their occupancy was higher than 61.2% (Table 1), verifying that these HBIs were stable throughout the entirety of the MR-GaMD simulations. Compared to the GDP-bound WT NRAS, G13D, Q61R, and C118S increased the occupancies of HBIs of the phosphate group from GDP with G13(D13), G15, and K16 but decreased that of the phosphate group of GDP with A18 (Table 1). Although G13D increased the occupancy of the hydrogen bond between S17 and the phosphate group of GDP by interacting with the GDP-bound WT NRAS, Q61R and C118S reduced the occupancy of this hydrogen bond ( Table 1). The residues N116 in loop L1 and S145, A146, and K147 in loop L4 yielded HBIs with the adenine group of GDP ( Figure 9B,D) with an occupancy higher than 57.9% (Table 1), which suggests that these hydrogen bonds were stable during the MR-GaMD simulations. In contrast to the GDP-bound WT NRAS, G13D, Q61R, and C118S enhanced the occupancies of HBIs of the adenine group in GDP with N116 and A146 but decreased those of the adenine group in GDP with S145 and K147 ( Table 1). The residues V29 and D30 from SW I and K117 in loop L1 formed HBIs with the middle ring of GDP and, in addition, D33 in SW I also produced an HBI with the middle ring of GDP in the GDP-bound G13D NRAS ( Figure 9C,D); however, the occupancies of these HBIs were lower than 25.0%. These results suggest that the HBIs of V29 and D30 in the SW I with GDP are highly unstable. The aforementioned interaction sites of GDP and NRAS agree well with the GDP binding spots identified in the work of Johnson et al. [29]. of the phosphate group of GDP with A18 (Table 1). Although G13D increased the occupancy of the hydrogen bond between S17 and the phosphate group of GDP by interacting with the GDP-bound WT NRAS, Q61R and C118S reduced the occupancy of this hydrogen bond (Table 1). The residues N116 in loop L1 and S145, A146, and K147 in loop L4 yielded HBIs with the adenine group of GDP ( Figure 9B,D) with an occupancy higher than 57.9% (Table 1), which suggests that these hydrogen bonds were stable during the MR-GaMD simulations. In contrast to the GDP-bound WT NRAS, G13D, Q61R, and C118S enhanced the occupancies of HBIs of the adenine group in GDP with N116 and A146 but decreased those of the adenine group in GDP with S145 and K147 ( Table  1). The residues V29 and D30 from SW I and K117 in loop L1 formed HBIs with the middle ring of GDP and, in addition, D33 in SW I also produced an HBI with the middle ring of GDP in the GDP-bound G13D NRAS ( Figure 9C,D); however, the occupancies of these HBIs were lower than 25.0%. These results suggest that the HBIs of V29 and D30 in the SW I with GDP are highly unstable. The aforementioned interaction sites of GDP and NRAS agree well with the GDP binding spots identified in the work of Johnson et al. [29]. Three conclusions can be drawn from the above results: (1) the salt bridge interaction of K16 and the HBIs of G13(D13), G15, K16, S17, and A18 with the phosphate group of GDP stabilize the conformation of the P-loop; (2) the salt bridge interaction of D119 and residues N116, S145, A146, and K147 with the adenine group of GDP lead to the ordering of the loops L1 and L4 compared to the switch domains; and (C) the instability in the HBIs of V29 and D30 with the middle ring of GDP is responsible for the structural disorder of the switch domains. Three conclusions can be drawn from the above results: (1) the salt bridge interaction of K16 and the HBIs of G13(D13), G15, K16, S17, and A18 with the phosphate group of GDP stabilize the conformation of the P-loop; (2) the salt bridge interaction of D119 and residues N116, S145, A146, and K147 with the adenine group of GDP lead to the ordering of the loops L1 and L4 compared to the switch domains; and (C) the instability in the HBIs of V29 and D30 with the middle ring of GDP is responsible for the structural disorder of the switch domains.

Initialization of Simulated Systems
The crystal structure taken from the Protein Data Bank (PDBID: 6ZIO) was used as the initial coordinates of the GDP-bound C118S NRAS [30]. To maintain the consistency of the atomic coordinates, the residue S118 in the 6ZIO was directly revised as C118 to generate the GDP-bound WT NRAS. Then, the G13 and Q61 in the GDP-bound WT NRAS were separately changed into D13 and R61 to yield the GDP-bound G13D and Q61R NRAS by using the Leap module in Amber 20 [81,82]. The Mg 2+ ion in the crystal structure was kept in the starting models of the four NRAS-related systems. The program H++ [83] was used to check and assign rational protonation states to the residues in NRAS. The chemical bonds between the missing hydrogen atoms in the crystal structure and the heavy atoms were constructed using the Leap module. The force field parameters of NRAS were obtained using ff19SB force field [84]. The force field parameters of GDP were taken from the work of Meagher et al. [85]. An octahedral periodic box of water with a buffer of 12.0 Å was adopted to solve each NRAS-related complex and the force field parameters of water molecules were assigned through the TIP3P model [86]. The appropriate numbers of sodium ions (Na + ) were placed around each complex in the 0.15 M NaCl salt environment to produce a neutral simulation system, in which the parameters of Na + , Cl − , and Mg 2+ were obtained from the Aqvist force field [87].

MR-GaMD Simulations
To obtain data with regard to the conformational changes of NRAS, MR-GaMD simulations are executed with the GDP-bound WT, G13D, Q61R, and C118S NRAS. All of the systems firstly underwent40,000step steepest descent minimization and then minimization using a3000step conjugate gradient, which was utilized to remove bad atom-atom contacts. Subsequently, the temperature of the four systems was enhanced from 0 to 310 K through a soft heating procedure within 1 ns in the NVT condition, in which non-hydrogen atoms were constrained by means of a weak harmonic restriction of 2 kcal·mol −1 ·Å −2 .After that, a 2ns equilibrium of 310 K was implemented to deeply relax the four systems in the NPT ensemble and a restriction similar to that for the heating process was also adopted. Finally, a cMD simulation was run for 100 ns to analyze the time evolution of the four systems at a constant temperature of 310 K and pressure of 1 bar, employing periodic boundary conditions and the particle mesh Ewald (PME) method [88,89]. Two new conformations randomly extracted from the previous 100ns cMD simulation were utilized as the initial coordinates to restart two new cMD simulations by randomly assigning initial atomic velocities for each conformation with the Maxwell distribution. The final structures obtained from the above three replica cMD simulations were adopted to run MR-GaMD simulations.
In this work, the GaMD simulations employed the harmonic boost potential to reduce the free energy barriers of the systems [60], which greatly enhanced the conformational sampling of the NRAS-related systems. The system potential V r was updated as V * r if V r was lower than a threshold energy E, which is clarified in Equation (1): in which k denotes the harmonic force constant, and the two parameters E and k may be adjusted with the help of the three enhanced sampling principles described in Equations (3) and (4): Ifa lower bound E = V max is assigned to E, then the parameter k 0 is generatedusing Equation (5): Inversely, if an upper bound E = V min + 1 k is given to E,then the parameter k 0 is obtainedthrough Equation (6): in which V max , V min , and V avg represent the minimum, maximum, and averaged potential energies of the systems, respectively. The parameter σ V indicates the standard deviation of the system potential energies, while the parameter σ 0 is a user-specified upper limit for accurately reweighting. For this work, 3.6 µs MR-GaMD simulations, composed of three separate GaMD simulations of 1.2 µs, were implemented with the GDP-bound WT, G13D, Q61R, and C118S NRAS. Three replica GaMD trajectories were connected into a singlejoined MR-GaMD trajectory (SJMT) so as to be convenient for the post-process analysis. The PyReweighting toolkit provided by Miao et al. was used to reweight the data from the post-processing analysis and to detect the original energetic profiles of the four NRASrelated systems [90]. As for all cMD and GaMD simulations, all chemical bonds between hydrogen atoms and heavy ones obeyed a restriction set by the SHAKE algorithm [91]. The temperatures of the four systems were regulated with a Langevin thermostat, in which a collision frequency of 2.0 ps −1 was adopted [92]. The PME method with an appropriate cutoff value of 10 Å was adopted for the estimation of non-bond interactions. The module pmemd.cuda included in Amber 20 was employed to run all simulations [93,94].

Principal Component Analysis
PCA was conducted to explore concerted motion couplings with functional significance in this work. The first step in the PCA was to build a covariance matrix C of the atoms Cα by using Equation (7): where q i and q j , respectively, signify the Cartesian coordinates of the ith and jth C α atoms fromNRAS, while q i and q j indicate their averaged positions over conformational ensembles. In general, the average is estimated by superimposing the SJMT with a defined referenced structure to remove overall translations and rotations through a least-square fit procedure [95]. The second step of the PCA was to perform diagonalization on the symmetric matrix C to yield a diagonal one A by means of an orthogonal coordinate transformation matrix T based on the following equation: in which the diagonal elements of A are the eigenvalues λ i and the columns of A correspond to the eigenvectors reflecting the motion direction relative to q i . The third step of the PCA was to explore the concerted movements of the structural domains in a multidimensional space using the eigenvector and to clarify the motion strength along an eigenvector, which could rationally reflect the conformational changes in structural domains from NRAS. The PCA was implemented by utilizing the CPPTRAJ module in Amber 20 [96]. The software VMD was employed to realize the visualization of the data generated by the PCA [79], plot pictures, and unveil the impacts of G13D, Q61R, and C118S on the conformational alterations of NRAS.

Dynamics Cross-Correlation Maps
To check the internal dynamics of the structural domains from NRAS, the elements C ij of DCCMs were calculated by using the x , y, and z coordinates of the C α atoms from NRAS with Equation (9): where ∆r i and ∆r j respectively indicate the displacement of atoms i and j relative to their averaged positions, while the angle brackets describe ensemble averages over the conformations saved at the SJMT. The element values of DCCMs fluctuated in a range from −1 to 1. The C ij of the positive and negative values characterized the positively correlated motions and the anti-correlated movements between atoms i and j, respectively. Color-coded patterns were employed to represent the extent of the correlated motions. The module CPPTRAJ in Amber 20 was utilized to perform the calculations of DCCMs.

Conclusions
NRAS is regarded as an important target in drug design for anticancer treatment. Mutation-mediated conformational transformation of the switch domains can induce alterations in the activity of NRAS. In this work, 3.6 µs MR-GaMD simulations, consisting of MR-GaMD simulations withthreeseparate1.2 µs replications, were performed on the GDP-bound WT, G13D, Q61R, and C118S NRAS to enhance conformational samplings and decipher the effects of mutations on the activity of NRAS. The analyses of the RMSFs and MSAs indicated that G13D, Q61R, and C118S not only change the structural flexibility of the switch domains but also reduce the global structural flexibility of NRAS. The DCMM calculations suggest that three mutations alter the correlated motions of the switch domains relative to the P-loop from the GDP-associated NRAS. PCA performed with the SJMT showed that the motion behavior of the switch domains from the GDP-bound NRAS was strongly disturbed. The switch domains are involved in binding domains of NRAS to effectors; hence, the alterations in the structural flexibility and motion behavior of the switch domains in NRAS produce certain influences on the binding of NRAS to effectors.
The FELs constructed by using the distance of D33 from residue 61 and the RMSDs of backbone atoms verified that G13D, Q61R, and C118S induce more energetic states in the GDP-associated mutated NRAS relative to the GDP-bound WT NRAS, which results in the switch domains in the GDP-bound mutated NRAS attaining a more disordered state than in the GDP-bound WT one. The switch domains are responsible for the NRASeffector binding; thus, the changes in the degrees of order of the switch domains certainly regulate the activity of NRAS. The analyses of the GDP-residue interactions verified that the instabilities in the HBIs of GDP with V29 and D30 in SWI drive the high disorder of the switch domains in NRAS, which plays an important role in regulations of the activity of NRAS. We expect that this study can contribute useful information for understanding the target role of NRAS in drug development for anticancer treatment.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27175596/s1, Figure S1: Interaction sites of GDP with NRAS and structural superimposition in the GDP-bound WT NRAS; Figure S2. Interaction sites of GDP with NRAS in the GDP-bound G13D NRAS located at the different energetic states; Figure S3. Interaction sites of GDP with NRAS in the GDP-bound Q61R NRAS located at the energetic states EV1-EV4; Figure S4. Superimposition of GDP and magnesium ions located at the energetic states EV1-EV4 of the GDP-bound Q61R NRAS; Figure S5. Representative structures of the GDP-bound C118S NRAS falling into the energy valleys EV1-EV3; Figure S6. Superimposition of structures situated at the energy valleys EV1-EV3; Figure S7. The functions of eigenvalues vs. eigenvector indexes from principal component analysis, used to describe the motion intensity of NRAS.