Previous efforts to induce local perturbations used generic heat baths to apply the perturbation to an individual residue [33],[34]. While sufficient for inducing small-scale changes, we wished to explore applicability to larger perturbations. In our experience, the application of a high-temperature heat bath propagates most of the energy along the backbone, which invariably unfolds the backbone at the point of perturbation. In contrast, the Rotamerically Induced Perturbation (RIP) circumvents the problem of the backbone unfolding by exclusively perturbing the χ angles, generating motions orthogonal to the backbone degrees of freedom.

In the RIP method, the protein is first stripped of ligands and waters, energy minimized, and pre-equilibrated without constraints to 300 K over 10 ps using Amber with GB/SA implicit solvent [35]. To avoid having to modify the Amber source code, the actual RIP calculations are performed by running sets of constant energy 100 fs simulations, which allow the rotamer perturbations to develop, with no other restraints applied. Between each interval, the instantaneous rotational velocity is calculated for each χ angle of the residue being perturbed. The bond and bond angle vibrations are suppressed, and a pure χ rotational velocity of the desired magnitude is applied to the atomic velocities and then the next 100 fs interval begins. Because of the reduced degrees of freedom, imposing a sidechain rotamer kinetic energy of 300 K results in far greater rotational velocities than those in a standard 300 K MD simulation. The direction of rotation is maintained until it exceeds a limit of ±60° from the initial χ angle, upon which the direction of rotation is reversed. This limit restricts the sidechain to exploring the basin surrounding a single rotamer, and allows long sidechains such as MET to effectively explore the landscape of a single rotamer. The resultant motion is that of the sidechain rotating back-and-forth around the initial χ angle. If there are no collisions with other residues, no energy gets transferred. But if there is a collision, a large displacement may be induced in another part of the protein. As the protein is simulated under constant energy during each interval, the transferred energy continues to propagate through the protein, resulting in a slight increase in the overall temperature.

In RIP, the rotational velocities are calculated directly from the rotational inertia of the sidechain. Thus we expect the rotational velocities of the χ angles in different sidechains to be different. For instance, at 300 K the phenyl ring in Phe should rotate more slowly about its χ2 angle than would the methyl-group in Ile about its χ2 angle. In order to demonstrate that RIP generates plausible χ angle behavior, we simulated the 17 amino acids that possess sidechains having χ-angles using a standard MD protocol. The amino acids were capped with methyl groups, and then simulated in AMBER with GBSA for 10 ps using a standard thermal bath at 300 K. The average values of the χ angle rotational velocities are then extracted from these trajectories, providing a reference set of rotational velocities for the χ-angles of each amino acid.

We then performed the RIP method on the same 17 amino acids. The average rotational velocities were extracted from the trajectories of the RIP simulations, and compared to the standard set of rotational velocities (Figure 2A). The Pearson correlation coefficient is 0.84, which shows that the RIP protocol generates the relative differences of the rotational velocities found in the standard simulations. The fit is better for low rotational velocities, which correspond to the motion of the heavier sidechains. The higher rotational velocities, which deviate more from a straight-line fit, correspond to the small sidechains or the ends of the long sidechains.

The differences between a residue perturbed by RIP and a residue regulated with a standard thermostat at 300 K can be shown in greater detail with the results for Ile (Figure 2B–D). As intended by the design of the RIP protocol, the average kinetic energy of the residue perturbed by RIP is equivalent to that in a standard simulation (Figure 2B). The differences can be found in the frequency distributions of the χ angles (Figure 2C and 2D), where RIP induces much higher rotational velocities than does the standard distribution. In the RIP rotational-velocity distributions, χ1 is peaked around 0.21 rad ps⁻¹ (Figure 2C) while χ2 is much flatter (Figure 2B). This is due to the smaller weight of the methyl group controlled by χ2, which moves much faster and collides more often with the rest of the Ile amino acid, thus broadening the distribution.

Triosephosphate Isomerase (TIM) has a ligand-binding loop that can close over the active site of the protein. In different crystal structures, this loop is observed in both an open and closed state [36] and NMR relaxation experiments detect loop fluctuations at a rate of 3×10⁴ s⁻¹ [11]. TIM typically exists as a dimer in which inter-monomer contacts are mediated through a separate dimer-interface loop.

The RIP analysis was applied to all residues in the TIM monomer having the open state of the mobile loop (Figure 1A). The simulations were also performed on the closed state of TIM and similar results were generated (discussed below). For each perturbing residue, the response of the protein can be measured by the Cα RMSD of the residues in the protein at the end of the 10 ps simulation. As an example, applying RIP to Glu128, a residue located near the mobile loop of TIM, results in a Cα RMSD deformation response with one large peak at residue 175 (Figure 3A). By overlaying the 10th ps conformation of the RIP simulation over the crystal structures (Figure 3B), it can be seen that the peak of large conformational change at residue 175 corresponds to the ligand-binding loop at residues 167–177. There is also a displacement of the helix near Glu128.

A global map of deformability in TIM can be constructed by applying RIP to every residue along the entire length of the protein. Each column in the RIP deformation map represents the 10 ps Cα RMSD response to a perturbation of RIP on the residue with sequential numbering on the X-axis (Figure 4A). To facilitate visualization, the actual response is shown on the map only if the responding residue is perturbed by more than a defined threshold. The threshold was chosen by analyzing the probability distribution (Figure 3C) of the Cα RMSD responses in the set of RIP simulations applied to every residue in TIM. This distribution has an average of 1.5Å, a σ of 1.2 Å. The peak of Cα RMSD at 0.8 Å corresponds to the residual motion from the initial pre-equilibration of 300 K. Using a cutoff of 6Å (mean+3.75σ) results in a sparse global map consisting of mostly contiguous segments of deformation (Figure 4A). Lower thresholds result in noisier maps, showing many single-residue deformations.

The first point to note is that there is no systematic response along the diagonal in the RIP deformation map. Residues adjacent to the perturbed residue are not automatically disturbed. This demonstrates the key property of the RIP method: perturbing a sidechain does not systematically disturb the local backbone, unless there is a specific interaction of the perturbed sidechain to the backbone. Consequently, if deformations are observed then they can be directly attributed to the perturbing sidechain.

As is evident from the TIM RIP deformation map, perturbing some residues can induce large changes in the protein while others have virtually no effect. The magnitude of the perturbation inducible by a particular residue can be quantitated by counting the number of residues that respond significantly to the perturbation (above the 6Å threshold). Residues capable of inducing significant perturbations will be referred to as structural linchpins (Figure 4B). While structural linchpins are generally larger amino acids, some smaller amino acids also show up.

Another way of extracting useful information from the RIP deformation map is to quantify the susceptibility of each residue to local perturbation (Figure 4C) by summing the number of above-threshold deformations horizontally across the deformation map (Figure 4A). Unlike the perturbation strength defined above, this local flexibility is a global property based on the entire set of perturbations along the protein chain.

In TIM, there are two segments with significant local flexibility (Figure 4C). One of the segments of large local flexibility corresponds to a striking horizontal band of deformations at residue 75 in the deformation map (Figure 4A). This segment corresponds to the dimer-interface loop. As deformation of the dimer-interface loop occurs independently of the location of the applied perturbation, a reasonable interpretation is that the dimer-interface loop is an intrinsically mobile loop in the monomer state. In contrast, there is a scattered band of flexibility in the segment corresponding to the ligand-binding loop at residue 175. Flexibility in the ligand-binding loop is induced only by the cluster of structural residues that surround it (Figure 4B). This can be interpreted as the ligand-binding loop being conditionally flexible, i.e., large fluctuations of neighboring residues can readily displace the mobile loop. To illustrate the range of motion induced by the RIP method in the mobile loop of TIM, the entire ensemble of the final conformations of each perturbation is shown (Figure 5A). The range of motion generated by RIP is far greater than the differences between the open and closed conformations of the mobile loop found in the crystal structures (Figure 1A).

Dihydrofolate Reductase (DHFR) catalyzes the reduction of dihydrofolate by NADP. The protein binds both ligands through the Met20 loop. Crystallography of the enzyme with different ligands has defined three states for the Met20 loop (Figure 1B) suggesting that loop dynamics may be functionally important: i) a closed state where the Met20 loop covers both dihydrofolate and NADP; ii) an occluded state where the Met20 loop packs against the dihdyrofolate but blocks NADP; and iii) an open state where the loop may be disordered in solution (Sawaya and Kraut 1997).

Using the RIP analysis, the results of the perturbations can be used to generate a RIP deformation map of DHFR (Figure 6A). From this map, potential structural linchpins can be identified (Figure 6B), and the local flexibility calculated (Figure 6C). The RIP deformation map displays many regions of large conformational changes, dispersed in sequence space. Regions of significant local flexibility (Figure 6C) map to the Met20 loop, the adenosine-binding loop, the F–G loop, and the G–H loop. Experiments suggests that all these loops couple to the activity of the Met20 loop [37],[38].

α-Lytic protease (αLP) (Figure 1C) is a kinetically stable protein shown by experiment to be highly resistant to unfolding [39] and to have extraordinarily high hydrogen-exchange protection factors (>10¹⁰) throughout its large hydrophobic core. The reduced dynamics and absence of mobile segments is believed to be functionally important for minimizing proteolytic destruction [40]. Reflecting these properties, αLP is also found to be extremely stable in standard MD calculations [33], where loops separating secondary structural elements appear to be firmly held by interaction with the body of the protein. This provides a good test for the RIP method, as any large-scale predictions of flexibility would likely be erroneous. In agreement with the known dynamics of αLP, the deformation map for αLP (Figure 7) is remarkably empty; there are only scattered segments of limited deformation, indicating that αLP is particularly insensitive to single residue perturbations and thus, conformationally very stable.

The Estrogen Receptor belongs to a family of nuclear receptors that are ligand-inducible transcription factors [15]. The ligand binds in a hydrophobic pocket in the ligand-binding domain of the Estrogen Receptor (ER) that is covered by Helix-12, so that in the bound state the ligand is completely buried. Not only does Helix-12 need to be displaced for ligand entry and exit, but it is also the key allosteric transducer. Crystal structures of ER with bound agonists and antagonists (Figure 1D) reveal that Helix-12 responds to the bound ligand by occupying one of several conformations that serve to either support or block binding of the downstream transcriptional co-activator protein [16]. Clearly, Helix-12 is a functionally-critical mobile segment of the protein although there is no information about the timescale of its motion.

The RIP deformation map of ER (Figure 8A) calculated in the absence of ligand, shows a very limited set of responding distortions with the exception of a horizontal band at residue 162. The horizontal band corresponds to a short loop between the α-helices (Figure 8C). Apart from a large conformational change of a long segment at the C-terminus (residues 232–250), ER appears to be very stable.

The large conformational change in the C-terminus is induced by Trp83, which is the only residue that would qualify as a structural linchpin (Figure 8B). Trp83 has been found to be a highly conserved residue in the nuclear receptor family with the greatest change in the solvent accessible surface area between the apo and holo structures [41], indicating a possible key role in binding Helix-12. The conformation induced by RIP on Trp83 reveals a dramatic cooperative motion in Helix-12, where the α-helix remains intact but moves by 13 Å (red in Figure 9). Although it doesn't move to the same location found in the crystal structure (green in Figure 9), perhaps due to the short simulation, we do see a hinge motion about the same pivot point (blue and green in Figure 9) observed in the crystal structures. Notably, this conformational change does not result in a significant distortion along Helix-12. The ability to generate this kind of large cooperative motion is due in part to the removal of the ligand normally found in the crystal structure for the simulations. In the absence of ligand, Helix-12 covers a cavity and is only weakly bound to the body of the protein. Thus the collisions generated by the perturbation on Trp83 are sufficient to induce the entire helix to move cooperatively.

Th N-terminal ligand-binding domain of the chaperone HSP90 [42] consists of a active-site lid that is 30 amino acids in length (Figure 1E). The lid is involved in ATP/ADP binding that is associated with the conformational changes of HSP90. The motion of the lid corresponds to an intermediate-scale motion of the lid independent of the body of the protein, as revealed by crystal structures of HSP90 in different ligand-bound states. The RIP deformation map of the HSP90 ligand-binding domain (Figure 10A) in the ADP-bound state, calculated in the absence of ligand, clearly identifies the lid as a highly flexible region. The conformational ensemble of the perturbed structures show large conformational changes in the lid comparable with the scale of motion found in the crystal structures (Figure 5D).

There are other systems that calculate protein flexibility from a crystal structure. To compare RIP to these methods (Figure 11), we have calculated the flexibility using the Anisotropic Network Model (ANM) [43] and CONCOORD [32]. We chose ANM as a representative of contact-based network models that generate domain-level motions, where ANM generates a theoretical B-factor for each residue from the sum of the contribution of the lowest modes of oscillations. We chose CONCOORD as a representative of local instability methods, where CONCOORD calculates an RMSDf fluctuation (RMSDf is the averaged Cα RMSD over different conformations of the same structure) at each residue from an ensemble generated from a monte-carlo simulation of the protein using distance constraints derived from a standard force-field. Both methods were chosen largely due to the easy availability of the software implementation. Finally, we compare these flexibilities to various experimental data. For the four proteins with large conformational changes. In the case of TIM, for example, the RMSDf was calculated from the position of the Cα atoms between the open and closed structures (Figure 11A) after they were aligned with the weighted superposition program Theseus [44].

The flexibilities calculated from the open conformation of TIM provide different results with respect to the RMSDf. All three measures of flexibility identify the dimer-interface loop (near residue 65) as flexible even though the RMSDf of the dimer-interface loop is negligible, as both the open and closed conformations exist in the same dimer arrangement in the crystal. For the ligand-binding loop (near residue 165), both RIP and CONCOORD identify elevated flexibility whereas ANM does not. CONCOORD also identifies several other regions as flexible, where there is no corresponding elevated RMSDf values. In contrast, RIP identifies as flexible only the ligand-binding loop and dimer-interface loop. As a further comparison, we show the flexibilities calculated from the closed conformation of TIM (Figure 11B). For the ligand-binding loop, only RIP identifies the loop as intrinsically flexible whereas both ANM and CONCOORD do not.

There exists a rich set of NMR measurements of DHFR that can be used to evaluate the calculated flexibilities of DHFR. From the averaged RMSDf between the open/closed/occluded conformations in crystal structures, only the ligand-binding Met20 displays any large conformation change. We find that both CONCOORD and RIP identify the Met20 loop as flexible whereas ANM does not. This can be contrasted with the S² parameters [22], which measures the backbone ¹⁵N mobility of DHFR on the picosecond-nanosecond timescale. We find that CONCOORD provides the best agreement to the S² parameters, as would be expected from previous studies [32] (Figure 11D). However, RIP also provides a reasonable match, whereas ANM does not. The similarity between RIP and CONCOORD is surprising as they do not agree in the other proteins. One possibility is that the nanosecond motions of the loops in DHFR correspond to microsecond fluctuations. This can be measured by NMR relaxation exchange factors R_ex factors [38] that detect significant backbone fluctuations between distinct states on the microsecond/millisecond timescale. R_ex measurements for DHFR identify microsecond dynamics in the Met20 loop, the F–G loop, the adenosine-binding loop and the G–H loop (Figure 11E), suggesting that in the case of DHFR, there is significant overlap between the nanosecond and microsecond dynamics. Coupling between these loops have also been identified in nanosecond MD simulations through cross-correlation and quasi-harmonic analysis [21]–[24].

The large conformational change in Helix-12 of ER, as evident in the large RMSDf values (Figure 11F), provide a difficult challenge for these models of flexibility. There are elevated values of RMSDf in short loops near residue 35 and 165 (Figure 11F), where all three flexibilities finds elevated flexibility. None of the flexibilities identify Helix 12 as a highly flexible region. However, with RIP we can identify one large conformational change by examining the deformability map (Figure 8A), where the conformational change is consistent with the crystal structures (Figure 9). Both ANM and CONCOORD identify elevated flexibility at residue 195, which do not correspond to elevated RMSD values.

The conformational changes in HSP90 is dominated by the motion of the lid, indicated by the large values of RMSDf at residue 110 (Figure 11G). This large conformational change is not identified by ANM but both RIP and CONCOORD identify the lid region as flexible. One advantage with RIP over CONCOORD is that by examining the conformational ensemble of perturbed structures generated by RIP (Figure 5E), we can see that the generated motions are of the order of the RMSDf observed in the crystal structures. Both CONCOORD and ANM identifies a region of elevated flexibility at residue 60 that does not correspond to any elevated region of RMSDf and the large RMSDf in the N-terminal that is involved with inter-domain interactions does not correspond to any calculated flexibility.

In conclusion, we find that the flexibility of RIP identifies only loop motions that correspond to large conformational changes of intermediate-scale motions. In contrast, CONCOORD identifies more regions as flexible, where there is some overlap with the regions identified as flexible by RIP. ANM performs poorly for intermediate-scale motions.

The RIP method is implemented as a PYTHON wrapper around the Sander package of AMBER [35] and all analysis code was written in PYTHON. The simulations of the RIP method are run in AMBER, using the PARM96 force-field with an GB/SA implicit-solvent term. To prepare for the simulation, ligands and crystallographic waters are removed from the crystal structure. The structure is then minimized and a Langevin thermometer is applied for a short equilibration at 300 K for 10 ps with a friction constant of 5 ps⁻¹.

The standard protocol for a RIP method lasts for 10 ps, which is long enough for large motions to be generated. At the beginning of the RIP method, the equilibrium value of the χ angles of the residue is stored. The run is then broken up into 100 fs intervals where each interval is simulated at constant energy.

Between each interval: (1) the direction of the rotational velocity of each χ angle is stored; (2) the atomic velocities of the residue is set to zero; (3) if the value of the χ angle exceeds 60° of the equilibrium χ value, the direction of the rotational velocity is reversed; (4) the magnitude of each χ rotational velocity is calculated from the sidechain conformation; (5) the χ rotational velocities are transformed into into atomic velocities and added to each atom; (6) the kinetic energy of the residue is scaled to the rotation temperature of 300 K. By scaling the atomic velocities, the kinetic energy of the residue is effectively transfered into the rotational modes of motion. This guarantees that even though the motion is artificially large, the amount of energy in the rotation is not more than would be available for the sidechain at equilibrium, even though this is unlikely to happen.

Between the intervals, a Python module translates the AMBER restart files into a Python object, from which the RIP protocol is used to generate new AMBER restart files for the next interval. Finally, the trajectories of all the intervals are spliced into a single trajectory. Since the modifications are made on the velocities, the coordinate trajectories are continuous.

In the RIP method, a rotational velocity for each χ angle of a sidechain is calculated at the beginning of every interval. From this rotational velocity, the atomic velocities are generated. To generate the the rotational velocities of the χ angles, each χ angle is assumed to be an independent degree of freedom. Based on the equipartition theorem, each independent χ angle can be assigned an energy E derived from the temperature T. This E is drawn randomly from a Gaussian distribution with mean energy ½kT and standard deviation √(½kT).

To convert a rotational velocity into an atomic velocity, a frame of reference for the axis of rotation must be chosen. As rotational velocities are only defined relative to the axis of rotation; rotations can occur on either end of the axis, and still give the same rotational velocity. Since the purpose of the RIP method is to minimize the motion of the backbone, only the sidechain atoms on the side of the rotation axis away from the backbone are rotated. Consequently, the rotational inertia of each χ angle, I = Σ mr², is calculated as the sum of the moment of inertia of these sidechain atoms, where r is the perpendicular radius of each atom from the χ angle axis of rotation.

To convert E into a rotational velocity ω, the equation of rotational energy E = Iω² is used. This is converted to a tangential velocity v through v = rω. This velocity is applied to the atom along the direction of the tangent to the axis of rotation. The atomic velocities due to each χ angle are then added cumulatively to each atom.

However, the different χ angles of the same sidechain do not represent completely independent degrees of freedom. As such, the final atomic velocities are re-scaled such that the total kinetic energy of the sidechain is E = 3/2 nkT where T = 300 K. This scaling only changes the magnitudes of the rotations and preserves the pure rotation around the χ angles.

In the analysis of the RIP simulations, rotational velocities of the χ angles need to be extracted from the trajectories. In the generation of rotational velocities, only atoms that are on the side of the rotation axis of the χ angle away from the backbone contribute to the rotational velocity. Therefore, in the extraction of the rotation velocities, only these atoms are considered. For each atom that fits the criteria, the tangential velocity v to the axis is calculated. This v is converted to a rotational velocity by ω = v/r where r is the perpendicular radius from the axis. As the contributions of each atom to the total rotational velocity of a χ angle depends on its moment of inertia, a weighting (w) for each atom is calculated from the moment of inertia I = mr² of the atom. The weighting is given by w = I / I_total where I_total is the sum of the I for all the atoms involved in the χ angle. The overall rotational velocity is then given by ω_total = Σ wω.