Protein conformational transitions explored by a morphing approach based on normal mode analysis in internal coordinates.

Large-scale conformational changes are essential for proteins to function properly. Given that these transition events rarely occur, however, it is challenging to comprehend their underlying mechanisms through experimental and theoretical approaches. In this study, we propose a new computational methodology called internal coordinate normal mode-guided elastic network interpolation (ICONGENI) to predict conformational transition pathways in proteins. Its basic approach is to sample intermediate conformations by interpolating the interatomic distance between two end-point conformations with the degrees of freedom constrained by the low-frequency dynamics afforded by normal mode analysis in internal coordinates. For validation of ICONGENI, it is applied to proteins that undergo open-closed transitions, and the simulation results (i.e., simulated transition pathways) are compared with those of another technique, to demonstrate that ICONGENI can explore highly reliable pathways in terms of thermal and chemical stability. Furthermore, we generate an ensemble of transition pathways through ICONGENI and investigate the possibility of using this method to reveal the transition mechanisms even when there are unknown metastable states on rough energy landscapes.

1 Introduction

Max Perutz and John Kendrew first determined the three-dimensional (3D) structures of hemoglobin and myoglobin in the 1960s, which laid the foundation for the field of structural biology [1-3]. Since then, numerous experiment-based studies have been performed to reveal structural information of macromolecules, resulting in more than 183,000 atomic-level structures in the Protein Data Bank (PDB) archive [4]. In addition, the vast array of information has demonstrated that regulated conformational changes are of crucial importance for proteins to perform their biological functions, which has led to increasing awareness of the need to probe these large transitions. Indeed, various experimental techniques such as nuclear magnetic resonance spectroscopy [5], small-angle X-ray scattering [6], and single-molecule spectroscopy [7] have been widely utilized to analyze the dynamic behavior of proteins. However, obtaining experimental information on the conformational changes of proteins is a longstanding challenge due to not only the intrinsic properties of the transition events with short-lived intermediate conformations, but also several technical limitations like sample preparation, system size, and time scale [8, 9].

Aside from the experimental studies, computational methods have played a key role in better understanding the functionally relevant dynamics of proteins that are difficult to capture through experimental approaches. Especially, molecular dynamics (MD) simulation, which samples conformational states in atomic detail by calculating interatomic forces using molecular mechanics force fields, has become one of the most powerful and popular tools [10, 11]. However, despite its successful applications in numerous studies, MD simulation has intrinsic limitations in exploring the large-scale conformational changes: the simulated systems easily get trapped in stable or metastable states and rarely cross high-energy barriers toward functional states, even on millisecond time scales. Recently, various MD strategies, such as development of special-purpose supercomputers [12, 13] and enhanced sampling methods [14-16], have contributed greatly to improving the performance of MD simulation, but the time-scale limitation remains to be resolved.

As an alternative approach to overcome the issue of computational complexity, normal mode analysis (NMA) has received much attention because it provides an efficient way to elucidate the intrinsic dynamics of proteins that are related to the global transitions [17-20]. NMA calculation is based on harmonic approximation of the potential energy function, and the resulting mode shapes are valid only near an equilibrium state. In other words, this method has inherent limitations in directly predicting conformational transitions that require inharmonic movements over energy barriers. Therefore, various methods combining NMA with other computational techniques have been developed to explore effective transition pathways between two end-point conformations [21-24].

In this study, we propose a new NMA-based pathway generation method called internal coordinate normal mode-guided elastic network interpolation (ICONGENI), an improved technique over the normal mode-guided elastic network interpolation (NGENI) [25]. The fundamental concept of both methods is to obtain intermediate conformations comprising a transition pathway by iteratively calculating displacement vectors to minimize error between the simulated intermediates and the targeted ones. In this process, NMA calculation is required to represent the displacement vectors as linear combinations of the lowest normal mode shapes and makes a critical difference between the two techniques: it is achieved in internal coordinates (IC-NMA) and Cartesian coordinates (CC-NMA) in ICONGENI and NGENI, respectively. CC-NMA has been widely used in studying protein dynamics due to the inherent nature of Cartesian coordinates (CCs): high computational efficiency and intuitive expression of protein dynamics but has the disadvantage of producing the mode shapes having unrealistic distortions like bond length stretching and bond angle bending [26]. On the other hand, IC-NMA has a distinctive advantage in describing the large-scale transitions in proteins. Internal coordinates (ICs) inherently facilitate the separatation of the torsion angles from the others, so that IC-NMA can be performed in torsion angle space. Given that the conformational changes are dominantly influenced by the variations in the torsion angles, not in the bond lengths and the bond angles that are nearly rigid, this strategy enables it to produce chemically relevant mode shapes preventing the unrealistic distortions of bond lengths and bond angles and extending the validity of the harmonic approximation in calculation [26-28]. In terms of computational complexity, IC-NMA is less efficient than CC-NMA because it has extra calculations involved in transformation either from CCs to ICs or from ICs to CCs (see section 2.2 for further details), but this is not a critical issue because both methods can be performed at the personal computer level. In other words, IC-NMA is more suitable than CC-NMA for describing and exploring the conformational changes in proteins, which provides an insight into the development of ICONGENI.

For validation of ICONGENI, we first demonstrate the superiority of IC-NMA in predicting large-scale conformational transitions by comparing the performance of ICONGENI to that of NGENI where the pathway is explored using CC-NMA [25]. Both methods are applied to two proteins: E. coli adenylate kinase (ADK) and E. coli ribose-binding protein (RBP). The comparative analyses of the distributions of ICs and the potential energies of the resulting pathways show that the transition pathways simulated by ICONGENI have higher thermal and chemical stability than those by NGENI. However, its efficient manner of computing intermediate structures has intrinsic limitations in exploring large-scale transitions on complex energy landscapes. To address this issue, ICONGENI generated a pathway ensemble for ADK dependent on the number of normal modes used in these simulations and characterized the ensemble in interdomain angle space, demonstrating that ICONGENI can explore plausible pathways on complex free energy landscapes.

2 Materials and methods

2.1 Protein structural information

To explore transition pathways of proteins through ICONGENI and NGENI, two end-point structures of each protein should be used as reference information. In the ADK case, the open and closed structures are chain A in PDB entry 4AKE (4AKE:A) [29] and 1AKE:A [30], respectively. In the RBP case, the open and closed structures are 1BA2:A [31] and 2DRI:A [32], respectively. In addition, we used several experimental intermediate structures of ADK whose PDB code 1ZIN, 1ZIO, 1ZIP [33], and 1DVR [34] to experimentally evaluate the ICONGENI simulation results. Because these intermediate structures have similar conformations, but different sequences with the reference structures (i.e., 4AKE:A and 1AKE:A), homology modeling was implemented using Modeller v9.25 [35]. In detail, 10 candidate models of the intermediate structures were constructed by using their 3D conformations (as templates) and the 2D sequences of the reference structures (as target proteins), and the best models for each template were selected based on DOPE score [36]. The selected structures were refined by energy minimization (500 steps of conjugate gradient) using the CHARMM36m force field [37].

2.2 Internal coordinate normal mode-guided elastic network interpolation (ICONGENI)

The ultimate goal of ICONGENI is to predict pathways for large-scale conformational changes of macromolecules based on structural information on two end-point (i.e., initial and final) conformations. To do this, valid displacement vectors can be obtained iteratively through ICONGENI, leading to determination of the consecutive intermediate conformations that comprise the pathways. A brief explanation of its algorithm is as follows.

First, IC-NMA is required prior to any procedure because we assume that the displacement vectors are linear combinations of low-frequency normal mode vectors. Next, the cost function is constructed to compute the degree of difference in interatomic distance between the resulting and desired conformations. Using the cost function, we devise compromise solutions (i.e., a series of weights assigned to the normal modes of the displacement vectors) between minimizing the values of the cost function and constraining the degrees of freedom (DOFs) of the structural dynamics with some low-frequency mode shapes. By repeating these cycles, promising transition pathways can be predicted. A detailed explanation of ICONGENI will be introduced in the following subsections.

2.2.1 Elastic network model (ENM)

ICs characterize molecular geometry using bond lengths, bond angles, and torsion angles terms to facilitate a better understanding of the structural dynamics of molecules. While the backbone bond lengths and bond angles are nearly fixed in molecular systems, some torsion angles (phi (ϕ) and psi (ψ) angles in the protein conformation) can vary, and their dynamics exert a strong influence on the large-scale conformational changes. To effectively describe molecular systems in ICs, we use a coarse-grained modeling method, the elastic network model (ENM), wherein specific atoms of a protein backbone (i.e., N, C, and C) are sampled and linked by a unit spring constant. The spring constant k matrix is defined as where k is a binary spring constant between atoms i and j, d is an actual distance between them, and d is a cutoff distance set to be 12 Å [19].

2.2.2 Normal mode analysis in internal coordinates (IC-NMA)

For IC-NMA calculation, the second derivatives of the kinetic and potential energy functions are required. The main strategy is to calculate the functions with respect to CCs and then to transform them into those with respect to ICs. All calculations of IC-NMA will be described in this section with respect to previous works [38-41].

2.2.2.1 Transformation from Cartesian to internal coordinates. To obtain the second derivatives of the kinetic and potential energy functions, the first derivatives of CCs with respect to ICs must be defined. Only torsion angles are regarded as variables for the calculation, while the bond lengths and bond angles are considered to be fixed. The following derivation of the first derivatives will be based on some assumptions. First, the Eckart condition is assumed to separate the internal and external motions because the ICs cannot express external motions [42]. For the position vector of atom i, the condition can be satisfied by the following equations: where m and are the mass and fixed position vector of atom i, respectively, in the reference conformation.

Next, let θ be a torsion angle around chemical bond α and domains A and B be the two domains divided by bond α as shown in . Then, the two domains are regarded as rigid bodies, based on which Eq (2) can be rewritten as where , and and are the position vectors of the center of mass of domains A and B, respectively.

Schematic of a molecular system composed of two rigid bodies A and B with a chemical bond α. The relative displacement between the two bodies can be defined by the torsion angle θ around the bond α. If a bond α links atoms (i−1) to i, t(α) designates i.

If and are the rotation vectors of domains A and B, respectively, their relative rotation vectors and δ are defined as follows: where .

Using Eq (5), d in domain B (corresponding to the second equation of Eq (6)) can be rewritten as

From Eqs (6) and (7),

Subsequently, using Eqs (3) and (6) and the concept of angular momentum, where the inertia tensors are given by and is the 3×3 identity matrix.

From Eqs (4), (5), and (8), where M = M+M.

Using Eqs (9) and (10), is expressed as where = +.

By substituting Eqs (10) and (11) into Eq (6), the derivative of CCs with respect to ICs is

Finally, these equations can be rewritten in the form of matrix-vector multiplication: where , and = [, ×].

2.2.2.2 Construction of the second derivative of kinetic energy. As shown in , the kinetic energy of the molecular system in ICs can be calculated using two torsion angle variables θ and θ as follows: where is the second derivative of kinetic energy for bonds α and β.

Schematic of a molecular system composed of three rigid bodies A, B, and C with two chemical bonds α and β.

If the bond α links atoms (i−1) to i, t(α) designates i.

Then, using Eq (13), can be defined as where and

2.2.2.3 Construction of the second derivative of potential energy. Basically, the potential energy of a molecular system is defined based on interatomic interactions (e.g., van der Waals bond and covalent bond potentials). With the aid of a simple assumption that one rigid-body domain is fixed instead of the Eckart condition, we can simply reformulate against Eq (13). If domain A is fixed in , the following equations are satisfied:

The second derivatives of potential energy in ICs can be derived from those in CCs as where is the potential energy function, and is the second derivative of potential energy for the ICs θ and θ (for the CCs and ).

is defined simply by the following function of interatomic distances [19]:

Then, supposing that domain B is fixed in the diagram of and using Eqs (16) and (18), can be formulated as

Using Eqs (15) and (19), we can construct the equation of motion with respect to ICs. For a molecular system having n torsion angle dynamics, where , and .

We can obtain normal modes (i.e., pairs of eigenvalues and eigenvectors) by solving a generalized eigenvalue problem for Eq (20). Next, an extra calculation is required to transform the resulting eigenvectors in ICs to those in CCs. Therefore, the final form of the normal mode vectors can be determined by the following equation [41]: where is the k normal mode vector of atom i (of bond α) with respect to CCs (ICs).

2.2.3 Construction of the cost function

The cost function is defined as a function of error in interatomic distances between the simulated conformations and the desired ones [25, 43]: where is the displacement vector of atom i describing conformational changes, and q is the desired distance between atoms i and j.

q in a target intermediate can be determined through linear interpolation between the two end-point conformations as where and represent the position vectors of atom i for the initial and final conformations, respectively. s is a proportional representation of the location of an intermediate to be simulated on the pathway when the total length of the pathway is set to 1.

In accordance with the strategy of ICONGENI, can be represented as the linear combination of a set of low-frequency normal mode vectors for reference conformations: where m denotes the number of low-frequency normal modes used in the simulation, and w is the weight of the n normal mode vector.

Using Eq (24), Eq (22) can be rewritten as where , and .

Then, to find the value of that minimizes C(), we simplify Eq (25) into the form of matrix-vector multiplication using a Taylor series expansion: where .

Then, Eq (27) can be written in the form of matrix-vector multiplication: where and

Finally, the optimal displacement vectors can be determined by solving for from the following equation:

The optimal displacement vector allows us to determine the intermediate conformation on the pathway. By repeating this calculation process using the simulated intermediate as a new reference, the transition pathways from initial to final conformations can be generated. The number of iteration steps was determined so that consecutive intermediate structures differ by a root-mean-square deviation (RMSD) of about 0.1 Å. Here, we set the number of iteration steps for the case of ADK (RBP) is set to 71 (62) because the RMSD value between the two end-point conformations is 7.16 Å (6.25 Å).

3. Results and discussion

3.1 Comparing the ICONGENI pathways to NGENI pathways

In this section, we discuss the effectiveness of ICONGENI by comparing the resulting pathways to those developed by NGENI [25] under the same conditions. Because the main difference between the two techniques is the coordinate space in which the NMA is performed (i.e., ICONGENI and NGENI are based on IC-NMA and CC-NMA, respectively), it is expected that this comparative analysis will demonstrate the superiority of IC-NMA in describing large deformations of proteins. We performed ICONGENI and NGENI for two proteins: ADK and RBP. ADK as a phosphotransferase enzyme catalyzing the reaction ATP + AMP ⇔ 2 ADP is composed of three domains: CORE, NMP, and LID, and undergoes two pairs of hinge motions of NMP and LID relative to CORE to fulfill its function [44]. RBP as one of the periplasmic binding proteins binds ribose through a hinge motion of two domains, which enables cells to sense and transport the ligand [31]. To predict the transition from open state to closed state, we obtained their 3D structures from PDB; the open and closed structures of ADK are chain A in PDB entry 4AKE:A and 1AKE:A, respectively, and those of RBP are 1BA2:A and 2DRI:A, respectively. The DOFs of these simulations were set to be the 50 lowest normal modes considered empirically sufficient to simulate the conformational changes within the experimental resolution based on our previous study [25]. For better understanding, the transition pathways explored by ICONGENI of ADK and RBP are provided in and , respectively.

First, the convergence issue of the pathways was addressed. To assess geometric convergence, we measured the RMSDs of the consecutive structures comprising the paths with respect to the final conformation (i.e., the closed structures) and judged that the paths satisfied the convergence condition if the RMSD values steadily decreased below certain thresholds that is selected as the smaller of the experimental resolutions of two reference structures (i.e., open and closed structures). As shown in , the ADK pathways of the two techniques had similar graphs and converged below the value of the resolution. Similarly, their RBP pathways also got close to the final conformation at a level beyond the resolution (). This confirmed that both techniques had no problem in terms of pathway convergence when generating the transition pathways based on the DOFs of the 50 lowest normal modes.

Conformational transition from open to closed states of (A) ADK and (B) RBP. Upper figures represent crystal structures of open and closed states of the proteins. ADK is composed of three domains: CORE (residues 1–29, 60–121, and 160–214), NMP (residues 30–59), and LID (residues 122–159). RBP is composed of Domain 1 (residues 1–103 and 236–264) and Domain 2 (residues 104–235 and 265–271). The lower graphs describe the convergence of simulated pathways of each protein by measuring changes in RMSD between predicted intermediates and the final structure. The results of ICONGENI and NGENI pathways are shown as red and blue lines, respectively. The black dotted lines represent the corresponding experimental resolution of the proteins.

Next, we investigated backbone bond length and bond angle distributions on the simulated pathways. Proteins undergo conformational changes mainly through variations of two types of backbone torsion angles: ϕ around the N−C bond and ψ around the C−C bond while variations of backbone bond lengths and bond angles are impractical during the transitions. The distributions of the bond lengths and bond angles provide key information to evaluate how well the proteins keep their molecular shape during large deformation. First, the backbone bond lengths are divided into three types: N−C, C−C, and C−N. According to the type, we calculated average (avg) and standard deviation (std) values over ICONGENI and NGENI pathways for two proteins (ADK and RBP) and analyzed their distributions using the experimental data as the reference avg and std values of the backbone bond length types () [45]. The avg values of the bond length in the ICONGENI paths were more concentrated around the corresponding experimental values than were those in the NGENI paths. Moreover, most std values of the bond lengths in the ICONGENI paths were distributed below the corresponding experimental values while many of those in the NGENI paths were higher than the experimental values. Subsequently, we investigated the backbone bond angle distributions (including N−C−C, C−C−N, and C−N−C) of the simulated pathways in the same manner as described above for analysis of the bond length distributions. Similar to the results in the bond length distribution graphs, the avg values of the bond angles in the ICONGENI paths were densely distributed around the corresponding experimental values, and the std values were distributed close to zero compared to the NGENI pathways (). These quantitative data commonly showed that the ICONGENI pathways were more likely to prevent unrealistic distortions in bond lengths and bond angles than were the NGENI pathways. On the other hand, there were exceptions to this principle, such as the small number of bond lengths in the ICONGENI paths with irrational avg or std values that deviated farther from the corresponding experimental ones than did those of the NGENI path (denoted by pink circles in ). We confirmed that all bond lengths that fell under these exceptions belonged to either the first or last residue of the proteins, which would imply that these unrealistic distortions were caused by the tip effect [41]. The tip effect is an inherent weakness of NMA (regardless of the type of coordinate system) in which the highly flexible parts in protein structures (e.g., hanging loops and protruding ends) exhibit abnormal behavior in some of the lowest normal modes. IC-NMA may suffer more from the tip effect than did CC-NMA because the mode shapes in IC-NMA is intrinsically limited in describing the movements of either the first or last residue where any torsion angle cannot be defined, which can explain the exceptions in . However, this limitation of ICONGENI is not a critical issue in predicting the transition pathways because the distortions of the tip parts could be considered as local vibrations that has little effect on the global motions.

Bond length distribution of the transition pathways of (A) ADK and (B) RBP. Comparison of the distributions of the avg and std values of backbone bond lengths in ICONGENI (denoted by red) and NGENI (denoted by blue) pathways. Both methods explore the transition pathways based on the DOFs of the 50 lowest normal modes. The bond length distributions are measured for three backbone bond length coordinates: N−C, C−C, and C−N. The green lines represent the corresponding experimental values of the coordinates. The pink circles represent specific cases where the ICONGENI pathway has irrational values of bond length avg or std.

Bond angle distribution of the transition pathways of (A) ADK and (B) RBP. Comparison of the distributions of the avg and std values of backbone bond angles in ICONGENI (denoted by red) and NGENI (denoted by blue) pathways. Both methods explore the transition pathways based on the DOFs of the 50 lowest normal modes. The bond angle distributions are measured for backbone bond angle coordinates: N−C−C, C−C−N, and C−N−C. The green lines represent the corresponding experimental values of the coordinates.

In the same context as investigating the bond length and bond angle distributions, the potential energies of the intermediate conformations comprising the resulting pathways were calculated. Because both simulation methods were carried out by using a coarse-grained modeling method, non-backbone atoms and O atoms in the backbone from reference structures (4AKE:A and 1BA2:A for ADK and RBP, respectively) were grafted to all intermediates and the generated all atom models were energy minimized within CHARMM36m force field for 500 steps of conjugate gradient [37] to eliminate any steric clashes and inappropriate geometries. Next, we calculated the potential energies of all intermediates of the NGENI and ICONGENI pathways by using CHARMM36m force field to quantitatively evaluate how the corresponding transitions are stable. shows the difference of the potential energy of each frame between the NGENI and ICONGENI pathways. From the results, we confirmed that the ICONGENI pathways are generally more stable (i.e., having lower potential energies) than the NGENI pathways. Furthermore, their energy gap increased as the pathways progress, which suggests that the qualitative difference between the pathways is increasingly noticeable in that the geometric errors are gradually accumulated when anharmonic transitions are explored by harmonic modes. Finally, these simulation results imply that the ICONGENI pathways are more reliable than the NGENI pathways in terms of thermal and chemical stability.

The difference between potential energies of ICONGENI and NGENI pathways for (A) ADK and (B) RBP. The difference of the potential energy ΔU = U−U. Before calculating the potential energies, all simulated intermediate structures were transformed into all atom models based on corresponding reference structures (4AKE:A (1BA2:A) for ADK (RBP)) and were energy minimized using CHARMM36m force field for 500 steps of conjugate gradient.

3.2 Predicting a transition pathway ensemble depending on a set of lowest normal modes

In the previous section, we tried to predict conformational transition pathways using ICONGENI with the 50 lowest normal modes considered sufficient to describe large deformation. Although the resulting pathways are shown to be reliable in terms of thermal and chemical stability, this study does not verify that ICONGENI can provide information on real transition trajectories. ICONGENI with the 50 lowest normal modes finds the deterministic and most effective pathways in terms of atomic displacements due to the intrinsic properties of the established cost function (see Section 2.3). In this section, we discuss the possibility that ICONGENI can predict plausible routes for conformational changes on complex energy landscapes by applying it to ADK of which transition mechanisms have been studied in numerous research works.

First, we generated an ensemble of the transition pathways of ADK through ICONGENI depending on the number of lowest normal modes (from 5 to 100), and their convergence was measured by RMSD with the closed state, as in the previous section. As shown in , the fewer are the normal modes used in the simulation, the less likely it is that the corresponding pathway reaches the final conformation. In detail, the pathways using fewer than 25 normal modes do not satisfy the convergence condition (i.e., their RMSD from the final conformation does not converge under the experimental resolution of ADK). This is not surprising given that the progression of pathways is influenced strongly by the DOFs describing the structural motions. In addition, this result suggests that it is necessary to focus on the “incomplete” pathways simulated with relatively few normal modes, as molecular systems usually explore seemingly inefficient routes of conformational transitions to arrive at functional states over several high-energy barriers.

The pathway ensemble for ADK generated by ICONGENI.

The ICONGENI transition pathways that make up the pathway ensemble are colored according to the lowest normal modes (from 5 to 100) used in the simulation (red to blue color scheme). (A) Convergence of the pathway ensemble. The RMSD values of each path relative to a final state are measured. The black dotted line represents the corresponding experimental resolution of ADK. (B) Projection of the pathway ensemble onto θ−θ space. The green points show the positions of experimental structures on θ−θ space. 4AKE:A and 1AKE:A indicate the open and closed states of ADK, respectively. 1ZIN:A, 1ZIO:A, and 1ZIP:A (1DVR:A, and 1DVR:B) indicate experimental structures at the NMPC state (the NMPO state).

Although the detailed transition mechanisms of ADK remain to be elucidated, previous experimental and theoretical studies have proposed several pathways via the NMP-closing/LID-opening (NMPC) state or the NMP-opening/LID-closing (NMPO) state [33, 34, 46–52]. In other words, the large-scale transition of ADK is characterized by interdomain hinge motions of NMP and LID relative to CORE. To delineate the NMP and LID movements on the ICONGENI pathways, we projected them onto interdomain angle space with the NMP-CORE angle (θ) and the LID-CORE angle (θ). θ (θ) is defined by the centers of mass of the backbone including N, Cα, and C in residues 115–125, 90–100, and 35–55 (179–185, 115–125, and 125–153) in the notation used in a previous study [51]. In addition, we used some experimental structures for cross-validation with our simulation results. 4AKE:A and 1AKE:A defines the two end-point conformations of the transition. The crystal structures whose PDB codes are 1ZIN, 1ZIO, and 1ZIP [33] and those whose PDB code is 1DVR [34] approximate the NMPC state and the NMPO state, respectively. On θ−θ space, the pathway ensemble has a tendency: the larger is the number of normal modes used in the simulation, the straighter are the resulting pathways to the closed state (). This is not surprising given that the established cost function is designed to produce the most efficient and direct paths within the given DOFs. The straight path on the interdomain angle space refers to the trajectory at which NMP and LID simultaneously open during the transition but is not favorable in terms of the free energy landscapes [47, 48]. As the number of modes used in the simulation decreased to less than 35, the resulting pathways tended to closely approach the NMPC state. When using significantly fewer normal modes (less than 10) for simulations, the corresponding pathways described the transition toward the more extreme NMPC state than did the crystal structures approximating the NMPC state (i.e., 1ZIN:A, 1ZIO:A, and 1ZIP:A). This result implies that the vibrational features describing the dynamics of θ were preferentially arranged in the lowest normal modes, resulting from the flexibility between NMP and CORE is higher than that between LID and CORE. Therefore, we suggest that the open-to-closed transition via the NMPC state is more plausible and reliable than that via the NMPO state in terms of the vibrational characteristics of ADK. However, this result does not mean ICONGENI always returns a single candidate (i.e., paths via the NMPC state) of the transition paths. If the normal mode set as the system DOFs is determined under certain conditions, ICONGENI could explore transition pathways via the NMPO state, which demonstrate that ICONGENI can explore multiple transition pathways compatible to several metastable states if information of the states is given (See more details in and ).

4 Conclusion

In this study, we introduced internal coordinate normal mode-guided elastic network interpolation (ICONGENI) as a theoretical method to explore the conformational transition pathways of proteins. By linearly interpolating the coarse-grained models of the two end-point states, ICONGENI defines virtual intermediate conformations of which the transition pathway is composed. Based on structural information, ICONGENI explores the optimal transition pathway (i.e., the pathway minimizing a cost function showing the error between the simulated intermediates and the virtual ones). When iteratively obtaining the consecutive conformations describing the transition pathway, the key idea of the method is to represent the displacement vectors as a linear combination of lowest normal mode vectors produced by normal mode analysis in internal coordinates (IC-NMA). Given that IC-NMA can describe chemically relevant dynamics (suitable for describing large-scale transitions) compared to NMA in Cartesian coordinates (CC-NMA), this strategy enables the proposed method to explore reliable transition pathways in an efficient manner.

To evaluate the superiority of ICONGENI, we performed comparative studies of ICONGENI with our previous method based on CC-NMA (named NGENI). For two proteins: adenylate kinase (ADK) and ribose-binding protein (RBP), we predicted transition pathways through the two methods under the same conditions (using the 50 lowest normal modes as the system degrees of freedom). The distribution data of the bond lengths and bond angles of the resulting pathways confirmed that these coordinates remained highly stable in the ICONGENI pathways compared to those in the NGENI pathways (Figs and ). Furthermore, we also calculated the potential energies of the simulated pathways and identified the energies of the ICONGENI pathways were lower overall than those of the NGENI pathways (). In conclusion, these results suggest that IC-NMA is suitable for representing realistic dynamics of the proteins, by extension, that ICONGENI could explore more reliable transition pathways than NGENI in terms of thermal and chemical stability.

Although ICONGENI using the degrees of freedom (DOFs) of the 50 or more lowest normal modes can provide a spatial understanding of conformational transitions, this approach is insufficient to explain the actual transition events on complex energy landscapes. To address this issue, we focused on a pathway ensemble for ADK simulated by ICONGENI. First, we confirmed that the more is the number of normal modes used in the simulation, the closer the initial structure is to the final one, which is not surprising because the number of normal modes directly indicates the DOFs to describe structural dynamics (). Next, we characterized the pathway ensemble by interdomain angles of ADK (i.e., the NMP-CORE angle (θ) and the LID-CORE angle (θ)) and found that the deficient pathways (using less than 50 lowest normal modes) provided meaningful insights into the conformational transitions of ADK. When projecting the ensemble onto θ−θ space, the deficient pathways showed the conformational transitions toward a metastable intermediate state (i.e., the NMP-closing/LID-opening state) while the sufficient pathways (using more than 50 lowest normal modes) showed those directly to the final state (i.e., the closed state) with unrealistic deformation (). Therefore, it is concluded that ICONGENI can explore meaningful transition pathways on complex energy landscapes.

The key role of computational approaches in investigating conformational transitions of proteins is to predict the trajectories that are beyond experimental capabilities. Our technique outlined here can shed light on the transition mechanisms in an efficient manner using only information on experimentally observed end-point structures. Furthermore, the simulation results strongly depend on a set of low-frequency normal modes as the system DOFs, enabling the method to generate a pathway ensemble based on dynamic characteristics and to provide low-energy paths. In this regard, our technique has the potential to find good candidates of unknown intermediate states on complex energy landscapes.

The ICONGENI simulations to explore ADK transition mechanisms via the NMPO state.

(DOCX)

Click here for additional data file.

The transition pathways for ADK via the NMPO state generated by ICONGENI.

Seven pathways in the vicinity of the NMPO state are shown on θ−θ space. The pathway named “diff. a degrees” means that it was explored by ICONGENI using the lowest normal modes that satisfy the condition: Δθ−Δθ>a (see details in ). The ADK crystal structures are taken as the references (indicated by green circles). 4AKE:A and 1AKE:A indicate the open and closed states of ADK, respectively. 1ZIN:A, 1ZIO:A, and 1ZIP:A (1DVR:A, and 1DVR:B) indicate experimental structures at the NMPC state (the NMPO state). The pathway ensemble data () is also included in this figure for comparison.

(TIF)

Click here for additional data file.

The transition pathway for ADK simulated by ICONGENI based on the DOFs of the 50 lowest normal modes.

(WMV)

Click here for additional data file.

The transition pathway for RBP simulated by ICONGENI based on the DOFs of the 50 lowest normal modes.

(WMV)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

7 Jul 2021

PONE-D-21-19500

Protein Conformational Transitions Explored by a Morphing Approach Based on Normal Mode Analysis in Internal Coordinates

PLOS ONE

Dear Dr. Kim,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

Please submit your revised manuscript by Aug 21 2021 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Wenfei Li, Ph.D.

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

Reviewer #3: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: I Don't Know

Reviewer #3: N/A

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors proposed an interesting path-sampling method to when end-point structures are given a priori. This method was developed based on a morphing algorithm and the normal mode analysis (NMA) and allows us to generate intermediate states between the end-point structures with low computational costs. The theoretical part was organized well by referring to their previous studies for this extension study. However, there might be some issues to be addressed prior to actual applications. Therefore, I would like to ask the following issues before publishing.

1) This method provides a set of (or an ensemble) of transition paths between given end-point structures. However, it is difficult to directly calculate free energy profiles based on this method when evaluating generated paths. In actual applications, users want to evaluate generated transition paths quantitatively on their free energy profiles. Herein, how to evaluate the generated paths quantitively using some physical quantities? Are there additional treatments for the quantitative evaluations? I would like to know a way to validate the generated transition paths.

2) In this study, the authors specified the 50 lowest modes as DOFs for the linier combinations. Why did you specify these modes? Are there any conditions to specify optimal DOFs to generate reliable transition paths? I think that users would like to know a way to specify the optimal DOFs.

3) Does this method work if there are multiple transition paths between given end-point structures? Does the current strategy return one of the multiple transition paths? Does this method always return a single candidate of transition paths? Or can this method always return the minimum free energy path (MFEP)?

4) For Eq (23), a target intermediate is determined by the linear interpolation between given end-point structures. Is this assumption enough to describe complicated transition paths of proteins? My concern is that this simple linear interpolation might include some artifact when we must consider complicated transitions of proteins between given end-point structures.

5) Related to previous issue (the 3rd issue), did this method successfully detect multiple transition paths of AK? In previous studies, several enhanced sampling methods detected multiple transitions path between the open and closed states of AK (for example, J. Chem. Theory Comput, 2019, 15, 5199-5208, Fig. 5). I would like to a correspondence of your result with previous study’s results.

6) Why did you adopt the ENM potential as a coarse-grained (CG) potential? Are there any possibilities to specify the other CG potentials? Generally, I think that the ENM potential is used to describe harmonic modes around a relaxed configurations using a reference structure such as an experimental (X-ray crystal or NMR) structure. Therefore, I think that this ENM potential might be unsuitable for describing large-amplitude domain motions. Why can this ENM potential treat or induce anharmonic motions of AK or RBP?

7) I think there might be a typo on Eq (17). It seems that the chain rule on Eq (17) does not work correctly in the present description.

Reviewer #2: Comments for the Authors

Lumping protein structures are undergoing large-scale conformational changes, which are essential to associate with their function at the cell and organism level, but have been elusive both experimentally and computationally. In this study, Lee et al. proposed a new computational methodology called internal coordinate normal mode-guided elastic network interpolation (ICONGENI) to sample intermediate conformations by interpolating the interatomic distance between two end-point conformations with the degrees of freedom constrained by the low-frequency dynamics afforded by normal mode analysis (NMA) in internal coordinates. The applications to two case studies (adenylate kinase (ADK) and ribose-binding protein (RBP)) by analyzing their open-to-closed transitions show that ICONGENI can explore highly reliable pathways in terms of chemical and thermal stability. Take together, the method illustrated in this study exhibits a certain degree of potentiality in terms of intermediated conformations predictions for these proteins with large-scale conformational changes during functioning. This study is a continuous work to the previous one, basically by changing the position descriptor of cartesian coordinates to that of internal coordinates. However, it’s superiority was not emphasized enough by the way of presenting/writing in this manuscript.

There are several points described below,

1. What are the computational costs for ICONGENI and NGENI methods? Is ICONGENI method much more computationally cheaper than NGENI? I would suggest to list limits and merits of two methods in a table.

2. The energy minimization details should be provided. Such as what is the minimization method, how many steps.

3. What is the cutoff distance used in the elastic network model (ENM)? What is the correlation between this parameter and the predicted results?

4. “For both proteins, the avg (std) values in the ICONGENI pathways were more concentrated around (below) the experimental values than were those in the NGENI pathways (Fig. 4).” Here the authors should provide more elaborate discussion regarding the reason behind this distribution difference.

5. “While the simulation result for RBP showed that ICONGENI suffered more from the tip effect than did NGENI, this was not a critical issue because the distortions of the tip parts could be considered local vibrations that barely contributed to the conformational transitions” Why did ICONGENI suffer more from the tip effect than did NGENI? Is this a limitation for this method potentially?

6.Are the structures of the intermediated states along the predicted pathways extractable? If so, their structures should be presented (if too many intermediated structures, the authors can make a movie). If not, why?

7. In Fig.3, “The black dotted lines represent the corresponding experimental resolution of the proteins.” Sorry, I don’t understand the description here. As the vertical axis is RMSD. I also didn’t quite understand the meaning of frame in the horizontal axis. Are they predicted intermediated conformations? Why did authors use ~ 70 frames for ADK but ~ 60 frames for RBP?

8. Page 22, line 425, a “than” should be added after “reliable”.

9.The color bar used in the Fig.6 is very difficult to differentiate, the rainbow color scheme might be better.

10. The superiority was not emphasized enough by the way of presenting/writing in this manuscript. More analyses and comparisons should be conducted.

11.Is the link of code available? If so, I suggest to provide the link in Abstract.

Reviewer #3: This manuscript describes a new methodology based on Normal Mode Analysis framework to decipher protein conformational transitions by quantitatively capturing often elusive intermediate conformational states. While the methodology has merit, I have several major comments that the authors must address to render the manuscript suitable for publication in PLOS ONE.

1. The authors describe in the Introduction (line 51-53) section that "MD simulation has some problems (like those of experimental methods)". This sentence is written rather vaguely because the context of comparing MD methodology to numerous experimental techniques has to be clearly laid out, which has not been done by the authors. The authors must provide context and provide clear contextual justification when they are comparing the technical details (which are distinct) in experimental techniques and computational techniques, including MD simulation.

2. The authors must elucidate (line 73-74) regarding the novelty of NMA obtained in Internal Coordinates (IC). How does the choice/assignment of IC affect the underlying analysis framework? The authors should comment on this.

3. Line 111: the authors must clearly state what they mean by conformational energy and how it is assigned to a conformational state.

4. Line 116: it is needed to justify why the authors chose DOPE score to select the best models corresponding to each template? Were any other scoring functions used? If not, then the authors must clearly weigh in on the merit of using the DOPE score as a selection criteria of their models.

5. Lines 127-128: The authors should comment on why curvilinear coordinates do not have a more significant impact on the relative spatial orientation of intermediate states.

6. The sentence in lines 153-154 should be re-written for clarity.

7. The authors should comment on why considering bond lengths and bond angles to be fixed does not affect the framework. While it is an approximation in their model, more quantitative justification is needed.

8. In Figure 1, the authors should comment on the validity of considering the two domains, A and B, as rigid bodies, and not consider local/distal fluctuations.

9. How does comparing the two approaches (previously developed, NGENI) and the current model (ICONGENI) ensure that the pathways considered are the same? Even setting all physical quantities in the respective calculations to be the same, how is it inherently ensured that they sample the same pathway? The authors must comment on/clarify this.

10. Lines 410-411 : The authors should clearly state what experimental structures were chosen for their experimental check? This sentence is ambiguous.

11. The authors should comment on how do they ensure that given two end point conformations, they are sampling the same pathway in multiple iterations. How will the framework apply if they are considering open/closed transition of a protein that has metastable intermediate states which may be sampled at certain physical conditions?

12. The Conclusion section seems to have a broad overlap with the Materials and Methods section, and as such it needs to be re-organized The authors must attempt to provide a comprehensive evaluation of the merits and limitations of their method when summarizing their models and its applicability to different types of systems.

In addition, the figures must be re-formatted as they are of rather poor resolution.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: Yes: Prithviraj Nandigrami

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

28 Aug 2021

First of all, we thank the referees for providing thoughtful and helpful suggestions and comments. Our responses to the comments are included in an uploaded file named "Response to Reviewers".

Submitted filename: Response to Reviewers.docx

Click here for additional data file.

6 Oct 2021

Protein Conformational Transitions Explored by a Morphing Approach Based on Normal Mode Analysis in Internal Coordinates

PONE-D-21-19500R1

Dear Dr. Kim,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Wenfei Li, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

This work reports a computational method for sampling the transition pathways of proteins between two functional structures, which can be useful in the studies of protein functional dynamics. Based on my reading and the recommendations by three referees, I am happy to recommend its acceptance for publication.

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: N/A

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors have carefully addressed all the issues pointed by the reviewers. The revised manuscript is publishable.

Reviewer #2: All of my previous comments have been addressed throughly. Thus, I would like to recommend the direct publication of the current version.

Reviewer #3: The authors have adequately addressed my concerns and comments. I would like to recommend publication of the manuscript.

One minor (optional) point I would like to mention: the authors should attempt to provide a comprehensive account of the merits and limitations of their current method, and compare it against similar method(s) available in the literature. This would most certainly help the readers weigh the applicability of the account presented in this manuscript more thoroughly.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: Yes: Wei Liu

Reviewer #3: No

27 Oct 2021

PONE-D-21-19500R1

Protein Conformational Transitions Explored by a Morphing Approach Based on Normal Mode Analysis in Internal Coordinates

Dear Dr. Kim:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Wenfei Li

Academic Editor

PLOS ONE