In recent years, there has been a growing interest to quantify the energy landscape that governs ribosome dynamics. However, in order to quantitatively integrate theoretical predictions and experimental measurements, it is essential that one has a detailed understanding of the associated diffusive properties. Here, all-atom explicit-solvent simulations (50 μs of aggregate sampling) predict that the diffusion coefficient of a tRNA molecule will depend on its position within the ribosome. Specifically, during aa-tRNA accommodation (i.e., the process by which tRNA enters the ribosome), the apparent diffusion coefficient decreases by approximately an order of magnitude. By comparing these to values obtained with an energetically “smooth” model, we show that the observed nonuniform behavior likely arises from electrostatic and solvation interactions between the tRNA and ribosome. These calculations also reveal the hierarchical character of ribosomal energetics, where steric interactions induce a large-scale free-energy barrier, and short-scale roughness determines the rate of diffusive movement across the landscape.
I. INTRODUCTION
Energy landscape principles have been central to developing a quantitative understanding of biomolecular dynamics.1–4 Classic examples include protein and RNA folding, which are now commonly described in terms of diffusion across a multidimensional free-energy surface.5–7 Within this framework, it has been shown that energetic roughness can manifest in the form of large-scale free-energy barriers, as well as configuration-dependent diffusion coefficients.8–10 The effectiveness of energy landscape techniques to understand folding suggests that applying these strategies to describe molecular assemblies11–13 may facilitate similar insights into the physical-chemical factors that govern biological dynamics.
One of the most well-studied biomolecular assemblies is the ribosome, which is composed of three large RNA molecules (∼100 to 2800 residues, each) and approximately 50 proteins. During the elongation cycle of protein synthesis, aminoacyl-tRNA (aa-tRNA) molecules are delivered to the ribosome. Initial association is followed by a large-scale conformational change called accommodation (Fig. 1), where the aa-tRNA molecule enters the ribosome. Accommodation allows the ribosome to increase the accuracy of gene expression by “proofreading” the incoming aa-tRNA molecule.14–16 That is, if an incorrect aa-tRNA molecule initially binds the mRNA-ribosome, accommodation provides a second opportunity for the molecule to be rejected prior to incorporating the amino acid into the nascent protein chain.
Aminoacyl-tRNA (aa-tRNA) accommodation on the ribosome. (a) The bacterial 70S ribosome is composed of 23S (white) and 16S (cyan) rRNA, as well as approximately 50 proteins (blue). During the elongation cycle, tRNA molecules (yellow, red) bind and transit through the ribosome. (b) When aa-tRNA is delivered to the ribosome by EF-Tu, it adopts the A/T conformation. (c) Simulated snapshot of A/T-like conformation of aa-tRNA, shown with the P-site tRNA (red), mRNA (green), and ribosomal helix 89 (white). In the current study, tRNA movement is described by the coordinate Relbow (distance between O3′ atoms of U8 in the P-site tRNA and U60 in the aa-tRNA, purple spheres). (d) During the first substep of accommodation, the aa-tRNA molecule adopts an elbow-accommodated conformation, where the aa-tRNA arm and 3′-CCA end are not yet bound to the A site.25,34 Molecular representations generated with Visual Molecular Dynamics (VMD).37
Aminoacyl-tRNA (aa-tRNA) accommodation on the ribosome. (a) The bacterial 70S ribosome is composed of 23S (white) and 16S (cyan) rRNA, as well as approximately 50 proteins (blue). During the elongation cycle, tRNA molecules (yellow, red) bind and transit through the ribosome. (b) When aa-tRNA is delivered to the ribosome by EF-Tu, it adopts the A/T conformation. (c) Simulated snapshot of A/T-like conformation of aa-tRNA, shown with the P-site tRNA (red), mRNA (green), and ribosomal helix 89 (white). In the current study, tRNA movement is described by the coordinate Relbow (distance between O3′ atoms of U8 in the P-site tRNA and U60 in the aa-tRNA, purple spheres). (d) During the first substep of accommodation, the aa-tRNA molecule adopts an elbow-accommodated conformation, where the aa-tRNA arm and 3′-CCA end are not yet bound to the A site.25,34 Molecular representations generated with Visual Molecular Dynamics (VMD).37
With the critical need for accurate protein synthesis, there have been many techniques developed to study accommodation. Kinetic measurements have isolated the rates of interconversion between biochemically distinct states,16–18 while cryo-EM and X-ray crystallography have provided structural snapshots of the endpoints of accommodation.19–23 In addition, the dynamics of individual aa-tRNA molecules have been monitored using single-molecule techniques.24 From a theoretical perspective, simulations have predicted that the structure of the ribosome requires accommodation to be a multistep process that begins with movement of the aa-tRNA elbow region.25,26 Subsequent cryo-EM reconstructions have corroborated the ability of aa-tRNA to adopt partially accommodated conformations.27 Together, previous investigations provide a wealth of information on the overall character of accommodation and proofreading, upon which rigorous theoretical analysis may be developed.
A long-term goal of ribosome studies is to dissect the numerous factors that govern biological function. To this end, simulations provide a suitable approach for establishing relationships between molecular structure, energetic interactions, and biological dynamics. As an example, the application of all-atom models that utilize simplified energetic representations (i.e., “structure-based” models) has demonstrated how elongation factor Tu may act as a crowding agent that reduces the free-energy barrier associated with aa-tRNA accommodation.28 Subsequent simulations showed how this effect may be further amplified by changes in flexibility of the L11 stalk.29 Similarly, simplified models were used to demonstrate how the size of each tRNA molecule can lead to differential hybrid-state formation rates.30 While a range of models may be applied to explore the many possible contributors to the free-energy landscape of elongation (e.g., changes in ionic concentrations and temperature, or introduction of post-transcriptional modifications), there will be a continuous need to compare theoretical predictions to experimental kinetics. As described below, these comparisons can be facilitated by knowledge of the diffusive properties of the ribosome.
By describing aa-tRNA accommodation in terms of diffusion along a one-dimensional surface, the free energy and rate of interconversion between any two points may be related according to4,31
where ρ is a reaction coordinate that describes aa-tRNA movement and D(ρ) is the diffusion coefficient along that coordinate. For the current study, the A/T and elbow-accommodated (EA) configurations are given by ρA/T and ρEA.32 However, in order to apply this relation, it is necessary to identify a suitable reaction coordinate and have estimates of the position-dependent diffusion coefficient D(ρ). With regard to the choice of coordinate, previous simulations have been used to identify interatomic distances that can measure the rate-limiting free-energy barrier.33 While earlier simulations aimed to quantify diffusive properties of tRNA,34,35 those studies described motion along a putative experimentally-accessible coordinate36 that was subsequently found to underestimate the free-energy barrier.33 In addition, previous investigations focused primarily on diffusion about the endpoints of accommodation and not the position-dependence of the diffusive properties.
Here, explicit-solvent simulations implicate a strong position-dependence of the diffusion coefficient of aa-tRNA inside of the ribosome. This is found to be in contrast to the relatively uniform values obtained when using a simplified energetic model. Based on this comparison, we infer that the position-dependence arises from energetic roughness associated with solvent and electrostatic effects. We further correlate changes in diffusion with the formation of tRNA-ribosome interactions, which provides a molecular basis for interpreting the position-dependent diffusive properties. In addition to providing a quantitative bridge between thermodynamics and kinetics, these measurements have additional practical utility, such as allowing for the estimation of effective time scales in coarse-grained/simplified simulations. Together, the presented analysis establishes a quantitative framework for describing the hierarchical energy landscape associated with ribosome function.
II. RESULTS
A. Position-dependent diffusion of tRNA
To estimate the position-dependent diffusion coefficient of aa-tRNA during accommodation, we employed molecular dynamics simulations with an all-atom explicit-solvent model. The motivation for using a highly detailed model is that the apparent diffusion coefficient is governed by the scale and distribution of short length-scale energetic roughness, as well as the viscosity of the solvent.4 Since these models explicitly describe the solvent, in addition to nonspecific roughness that is due to electrostatics and vdW interactions, they provide a suitable representation for exploring the potential position-dependence of tRNA diffusion.
To obtain sampling over a broad range of aa-tRNA configurations, in the context of the accommodation corridor,25 numerous simulations were initiated from distinct points along the accommodation process. For each initial configuration, the system was first equilibrated and then 305 simulations were performed for 15 ns (50 μs of aggregate sampling, Fig. S1). Each simulation was initialized with random velocities drawn from a Maxwell-Boltzmann distribution. To describe the dynamics of aa-tRNA, we used the coordinate Relbow: the distance between the O3′ atoms of U60 in the aa-tRNA and U8 in the P-site tRNA (Fig. 1). This coordinate has previously been shown to capture the elbow-accommodation barrier33 that is sterically imposed by ribosomal Helix 89 (H89).28 While the current short-time scale simulations are not sufficient to directly observe large-scale displacements, this data set allows us to quantify the statistics of local relaxation processes that are characteristic of diffusive motion.
Several different analysis strategies reveal a robust position-dependence of the diffusion coefficient, where D decreases sharply as the aa-tRNA initially enters the ribosome [Fig. 2(c)]. The primary reason for applying multiple methods to calculate D is that the diffusive description is an approximation, by definition. Accordingly, without an analytical expression for D, a variety of kinetic properties may be used to infer its value. In addition, with the finite size of the simulated data set, it is possible that each estimate may be influenced differently by statistical noise. Fortunately, we find consistent trends across methods, as described below.
Diffusion coefficients determined from explicit-solvent simulations. (a) Displacement-squared as a function of lag time (Δt), for different ranges of initial Relbow values [Eq. (3)]. Each line corresponds to a different range of initial values of Relbow(t) (increments and widths of 1 Å). The slope was calculated by applying a linear fit to Δt = 1–5 ns (linear correlation coefficient 0.95). The simulated data set contained over 50 μs of aggregate sampling (Fig. S1). (b) Variance of Relbow, as a function of lag time, calculated based on different initial values of Relbow [Eq. (4)]. An average linear correlation coefficient is ∼0.95. (c) Values of the diffusion coefficient, as a function of Relbow based on Eq. (2) (black), Eq. (4) (red), or a Bayesian Inference approach8 (blue). All three methods reveal an abrupt decrease in D between 45 Å < Relbow < 55 Å.
Diffusion coefficients determined from explicit-solvent simulations. (a) Displacement-squared as a function of lag time (Δt), for different ranges of initial Relbow values [Eq. (3)]. Each line corresponds to a different range of initial values of Relbow(t) (increments and widths of 1 Å). The slope was calculated by applying a linear fit to Δt = 1–5 ns (linear correlation coefficient 0.95). The simulated data set contained over 50 μs of aggregate sampling (Fig. S1). (b) Variance of Relbow, as a function of lag time, calculated based on different initial values of Relbow [Eq. (4)]. An average linear correlation coefficient is ∼0.95. (c) Values of the diffusion coefficient, as a function of Relbow based on Eq. (2) (black), Eq. (4) (red), or a Bayesian Inference approach8 (blue). All three methods reveal an abrupt decrease in D between 45 Å < Relbow < 55 Å.
In the first approach, we used the relation
where
and ⟨⟩ denotes an average over all t for which |Relbow(t) − R| < 0.5 Å. The first 2 ns of each simulation were excluded from analysis. The mean-squared displacement was found to be linear over the interval 1 ns < Δt < 5 ns [Fig. 2(a)], from which the diffusion coefficient was estimated. We did not fit to Δt < 1 ns since these rapid motions likely reflect intramolecular fluctuations, such as deviations of individual dihedrals, rather than collective displacements of aa-tRNA. A striking feature of D(R) is that it displays a sharp decrease as the aa-tRNA enters the ribosome [Fig. 2(c), black]. Specifically, while the values of D for the A/A ensemble (Relbow ≈ 28 Å) are comparable to previous estimates for a different interatomic distance,34 the current analysis reveals an order-of-magnitude decrease in D as the aa-tRNA moves from 56 Å (A/T ensemble) to 46 Å (partially accommodated).
We next estimated the diffusion coefficient using the short-time solution to the Fokker-Planck equation, as employed previously to study peptide folding.38 In the short-time regime, D may be approximated according to
where σ2(t) and σ2(t + Δt) represent the variance in Relbow taken over all t for which |Relbow(t) − R| < 0.5 Å. Consistent with the estimates obtained using Eq. (2), the variances also predict nonconstant behavior, where D decreases from ∼10 to 1 μm2/s during accommodation [Fig. 2(c), red].
The third method we applied is based on Bayesian Inference to estimate a transition matrix that describes the system in terms of discrete transitions along the reaction coordinate.8 Consistent with the estimates obtained using Eqs. (2) and (4), this approach also implicates an order-of-magnitude decrease in D as the tRNA enters the accommodation corridor [Fig. 2(c), blue].
B. tRNA confinement does not influence diffusion
To explore the origins of position-dependent values of D, we asked whether this behavior can arise from ribosome-induced confinement of the aa-tRNA molecule. For this, we compared the values obtained with the explicit-solvent model (Fig. 2) to those calculated from a simplified model in which explicit energetic roughness is not included (Fig. 3). Specifically, we analyzed the simulated dynamics of aa-tRNA when using an all-atom structure-based model (see methods). In the structure-based model, all nonhydrogen atoms are explicitly represented, and the geometry of each molecule is maintained through bonds, bond angles, and improper dihedral angles. Attractive interactions are only defined between atoms that are in close proximity in the A/A conformation, and the potential energy minimum of each interaction is defined to stabilize the fully accommodated configuration of the aa-tRNA. For a detailed review on the construction of these models, see the work of Levi et al.39 While the structure-based model lacks nonspecific energetic roughness (e.g., electrostatics), the steric composition of the ribosome introduces a free-energy barrier of approximately 6kBT, which is centered around Relbow = 38 Å [Fig. 3(a)].
Diffusion coefficient determined from an energetically “smooth” structure-based model. (a) While the structure-based model (described in Ref. 29) lacks explicit energetic roughness arising from non-A/A interactions, there is a pronounced free-energy barrier centered about Relbow = 38 Å. (b) Diffusion coefficient as a function of position based on Eq. (2) (black), Eq. (4) (red), or a Bayesian Inference method (blue). While there are minor variations, all methods implicate a relatively uniform value of D, which is in stark contrast to the values obtained from explicit-solvent simulations (Fig. 2).
Diffusion coefficient determined from an energetically “smooth” structure-based model. (a) While the structure-based model (described in Ref. 29) lacks explicit energetic roughness arising from non-A/A interactions, there is a pronounced free-energy barrier centered about Relbow = 38 Å. (b) Diffusion coefficient as a function of position based on Eq. (2) (black), Eq. (4) (red), or a Bayesian Inference method (blue). While there are minor variations, all methods implicate a relatively uniform value of D, which is in stark contrast to the values obtained from explicit-solvent simulations (Fig. 2).
Since steric interactions can lead to a pronounced free-energy barrier, it is also possible for excluded volume interactions to introduce energetic roughness that would manifest in the form of a position-dependent diffusion coefficient. However, we find that the diffusion coefficients obtained with the structure-based model are relatively uniform (∼20% variation) for all values of Relbow [Fig. 3(b)]. The lack of significant changes in D with the energetically smooth model supports the interpretation that position-dependent diffusion arises from detailed electrostatic and solvation interactions, rather than confinement of the tRNA molecule.
C. tRNA-ribosome contacts rationalize diffusive properties
Since the diffusion coefficient is position-dependent, one can infer that the short-scale energetic roughness depends on the position of the aa-tRNA molecule. According to this interpretation, one may anticipate that changes in the number and distribution of tRNA-ribosome interactions will correlate with changes in D. To address this, we first calculated the average solvent accessible surface area (SASA) of the aa-tRNA, as a function of Relbow, for the explicit-solvent simulations. We focused our SASA calculations on the elbow, arm, and 3′-CCA end [shown in orange in Fig. 4(a)] and not the anticodon stemloop (ASL) since the ASL undergoes minimal displacements and remains bound to the ribosomal A site throughout all simulations. As expected, due to an increase in tRNA-ribosome interactions during accommodation, there is a decrease in the mean SASA value (ΔSASA ∼ −400 Å2) as the aa-tRNA accommodates [Fig. 4(b)]. Interestingly, the sharpest decrease in SASA values nearly coincides with the sharp decrease in D, which is consistent with an increase in energetic roughness arising from tRNA-ribosome contacts.
Solvent Accessible Surface Area (SASA) of aa-tRNA correlates with the diffusion coefficient. (a) The surface area of each atom was calculated based on all tRNA and ribosomal atoms. From this, we evaluated the average SASA associated with atoms in the aa-tRNA elbow, arm, and 3′-CCA end (residues 1–24, 44–76 shown in orange). (b) The average SASA value (standard deviation shown as error bars), as a function of Relbow, reveals a decrease as the aa-tRNA initially accommodates (Relbow = 46–48 Å), which is consistent with the observed decrease in D [Fig. 2(b)].
Solvent Accessible Surface Area (SASA) of aa-tRNA correlates with the diffusion coefficient. (a) The surface area of each atom was calculated based on all tRNA and ribosomal atoms. From this, we evaluated the average SASA associated with atoms in the aa-tRNA elbow, arm, and 3′-CCA end (residues 1–24, 44–76 shown in orange). (b) The average SASA value (standard deviation shown as error bars), as a function of Relbow, reveals a decrease as the aa-tRNA initially accommodates (Relbow = 46–48 Å), which is consistent with the observed decrease in D [Fig. 2(b)].
To identify the specific interactions that are associated with decreases in SASA and D values, we next quantified the timing and distribution of tRNA-ribosome contact formation. Specifically, we evaluated Pi(R): the probability that the aa-tRNA is in contact with residue i of the ribosome, as a function of the tRNA position. Here, we consider a contact to be formed if the interatomic distance between the tRNA and any atom in residue i is less than 4 Å. Consistent with a sharp decrease in D as the tRNA enters the ribosome, there are concomitant increases in Pi(R) for specific residues. When the ribosome is in the A/T ensemble (Relbow = 56 ± 1 Å), there are minimal interactions between the aa-tRNA and ribosome, and there is a relatively small probability (Pi < 0.5) of forming contacts with H90 and proteins S12 and L14 [Fig. 5(a)]. Since these interactions are distal to the tRNA elbow, it is not surprising that the diffusion coefficient along Relbow is rather large for this ensemble. However, as the aa-tRNA enters the ribosome (Relbow = 48 ± 1 Å), there is a high probability (>0.95) of forming contacts between the aa-tRNA elbow and H89 [circled in Fig. 5(b)]. An increase in contact probability suggests these interactions are primarily responsible for the initial decrease in D as the aa-tRNA enters the accommodation corridor. As the aa-tRNA molecule continues to further accommodate into the ribosome, the elbow region maintains a high frequency of contacts with the ribosome (Fig. S2), consistent with the low values of D for accommodated configurations.
Formation of tRNA-ribosome contacts rationalizes position-dependent diffusion. (a) Representative simulated snapshot (aa-tRNA, yellow; 23S rRNA, white tube; proteins, blue tubes) for Relbow ∼ 56 Å. Residues shown as spheres form contacts (4 Å cutoff) with the tRNA in at least 1% of the simulated frames for which Relbow = 56 ± 1 Å (A/T ensemble). Spheres are colored by the probability of each residue forming contacts with aa-tRNA. Proteins S12 and L14, as well as helix 90, show moderate probabilities (<0.5) although the tRNA elbow makes few contacts. (b) Same as panel (a), for Relbow = 48 ± 1 Å. There is a high probability (>0.95) of forming contacts with ribosomal helix 89 (circled), consistent with lower values of D. As accommodation continues and the tRNA enters the A site (i.e., lower values of Relbow), a high frequency of contacts is maintained (Fig. S2), consistent with lower values of D.
Formation of tRNA-ribosome contacts rationalizes position-dependent diffusion. (a) Representative simulated snapshot (aa-tRNA, yellow; 23S rRNA, white tube; proteins, blue tubes) for Relbow ∼ 56 Å. Residues shown as spheres form contacts (4 Å cutoff) with the tRNA in at least 1% of the simulated frames for which Relbow = 56 ± 1 Å (A/T ensemble). Spheres are colored by the probability of each residue forming contacts with aa-tRNA. Proteins S12 and L14, as well as helix 90, show moderate probabilities (<0.5) although the tRNA elbow makes few contacts. (b) Same as panel (a), for Relbow = 48 ± 1 Å. There is a high probability (>0.95) of forming contacts with ribosomal helix 89 (circled), consistent with lower values of D. As accommodation continues and the tRNA enters the A site (i.e., lower values of Relbow), a high frequency of contacts is maintained (Fig. S2), consistent with lower values of D.
D. Quantifying a hierarchical energy landscape
Calculating the position-dependent diffusion coefficient also allows one to quantify the magnitude of short-length-scale energetic roughness. While simplified structure-based models of the ribosome implicate a large-scale free-energy barrier (Fig. 3 and Refs. 28, 29, and 33), the current explicit-solvent simulations implicate inhomogeneous short length-scale energetic roughness. If one assumes the short length-scale energetic roughness is Gaussian distributed, then the apparent/effective diffusion coefficient may be given by40
where ϵ is the width of the distribution and Dfree is the diffusion coefficient in the absence of roughness (i.e., in solution). Scattering experiments have implicated translational diffusion coefficients of free tRNA in solution to be ∼50 μm2/s.41 Since motion along Relbow is associated with translational displacements that are correlated with the center-of-mass position,28 we will use this experimental value as an estimate of Dfree. In the current study, we find D ∼ 10–20 μm2/s for the A/T ensemble (large values of Relbow). According to Eq. (5), this would indicate a minimal level of energetic roughness (ϵ = 1–1.3kBT) in the A/T ensemble, which is consistent with a small number of intermolecular contacts [Fig. 5(a)] and an elevated SASA value [Fig. 4(b)]. As the tRNA enters the ribosome, the diffusion coefficient decreases to ∼0.5-4, indicating an increase in the roughness to ∼1.6-2.1kBT. In addition, this increase is predicted to occur prior to the aa-tRNA molecule reaching the peak in the free-energy barrier [Relbow ∼ 38 Å, Fig. 3(a)].
From a qualitative standpoint, these observations begin to reveal a hierarchical energy landscape, where the large-scale sterically induced barrier is on the order of 6–10kBT,28,29,33 while the short-length roughness in on the scale of 1–2kBT. Furthermore, these calculations indicate that, as the aa-tRNA attempts to overcome the large-scale free-energy barrier, it will encounter increasing energetic roughness that decreases D. This decrease will then amplify the effect of the free-energy barrier on kinetics. Finally, the presence of inhomogeneous roughness indicates that the rate of accommodation will be sensitive to the position of the rate-limiting barrier, as well as the height. That is, the rate will be slower if the barrier is centered about smaller Relbow values (i.e., later in the accommodation process). Accordingly, when attempting to quantify the relationship between molecular structure and kinetics, it will be necessary to identify the precise position of the thermodynamic free-energy barrier.
E. Estimating effective time scales in models
To illustrate a practical utility of diffusion calculations, we will use the above analysis to demonstrate how one may estimate time scales when using simplified and/or coarse-grained molecular models. While simpler models can allow for longer effective time scales to be simulated, one must often rely on indirect methods to estimate the simulated time. To this end, it is sometimes possible to calculate thermodynamic free-energy barriers28,29,33 and then numerically evaluate the rates via Eq. (1). However, calculating free-energy barriers is not always tractable for large assemblies, and it can be difficult to determine which reaction coordinates are kinetically relevant. An alternate strategy is to compare mean-first passage times between a specific model and experiments, in order to obtain a correction factor for the simulated time.42 While this approach is useful for predicting time scales within a class of related molecules/motions, a limitation of this strategy is that the correction factor can mask any disagreement between experimental and theoretical barrier heights. To further complicate these estimates, the correction factor will also depend on simulation parameters, such as the applied drag coefficient when employing Langevin Dynamics protocols. To help circumvent these challenges, we suggest a more direct approach for estimating time scales, where knowledge of experimental barrier-crossing rates is not required. Specifically, we will infer effective time scales through a comparison of the apparent diffusion coefficients.
Before discussing the comparison of diffusion coefficients, it is valuable to briefly describe reduced units (ru). When using any model, after defining the units of energy (e), mass (m), and length (l), the corresponding time unit will be given by
In an all-atom structure-based model, experimental measurements of base-stacking interactions suggest the modeled energy scale e is approximately 1 kCal/mol (discussed in Ref. 26). This particular variant of a structure-based model defines the mass of each nonhydrogen atom to be 1, implying a mass scale of ∼12 amu, and the employed length scale is 1 nm. According to Eq. (6), these values implicate a reduced time unit of ∼3 ps. While it is common to directly report predicted time scales in these units, it is important to recognize that the long-time simulated dynamics will be associated with an elevated diffusion coefficient since solvent and energetic roughness is not modeled. Thus, when studying long time dynamics (e.g., large-scale barrier-crossing events), it is appropriate to expect a general “speed up” of the diffusion that will result in a longer effective simulated time scale. In the context of protein folding, Kouza et al. used lattice models and off-lattice coarse-grained G-like models to estimate that one reduced unit corresponds to an effective time scale of approximately 3–10 ns.42 While this is a useful estimate when employing that particular coarse-grained model (and simulation protocols) to study folding, there is no justification for applying this factor when simulating the ribosome with an all-atom model.
To estimate the effective simulated time unit in a coarse-grained/simplified model, τru, we will directly compare the calculated diffusion coefficient in a structure-based model with that obtained from an explicit-solvent model. As described above, we have calculated the diffusion coefficient as a function of aa-tRNA position using an explicit-solvent model [Fig. 2(c)], and it was found to be in the range of 1–20 μm2/s. In comparison, when using the structure-based model, the calculated diffusion coefficient was approximately 0.1 Å2/. Equating these values of the diffusion coefficient implicates a reduced time unit of ≈50 ps to 1 ns. These are 1–2 orders of magnitude smaller than values based on a coarse-grained model for folding,42 which highlights the context-dependence of effective time scale estimates. Unfortunately, since previous studies of the ribosome utilized the folding-inspired factor, those reported effective time scales were likely significantly overestimated.26
In closing, one should note that this approach to estimating time scales in simplified/coarse-grained simulations may be extended to other models and processes. As other variations of models are developed and applied to assemblies, the comparison of diffusion coefficients may serve as a direct approach for estimating effective time scales that accounts for the effects of solvent and energetic roughness.
III. SIMULATION DETAILS
A. Explicit-solvent model
Simulations were performed using Gromacs v5.1.4.43,44 The AMBER99SB-ILDN force field45 was used with the TIP3P water model.46 We performed explicit-solvent simulations initialized from 11 configurations that represent different values of Relbow (27, 29, 31, 33, 35, 37, 39, 43, 47, 51, 55 Å). The initial structures were randomly chosen from previously reported simulations of accommodation that employed an all-atom structure-based model.29 Since spontaneous and reversible accommodation rearrangements are accessible with the all-atom structure-based model, the selected configurations describe a sterically preferred route that also preserves proper stereochemistry (i.e., excluded volume, hybridized orbital geometries, and base pairing). For each configuration, the following preparation and equilibration protocol was applied: First, the system was solvated in a cubic box, where the shortest distance between the solute and box edge was set to 10 Å (dimension of approximately 170 Å, for a total of ∼500 000 atoms). K+ ions were then added to neutralize the system, and excess K+, Cl−, and Mg2+ were added, such that and . Steepest descent energy minimization was first applied with position restraints on all heavy atoms. Steepest descent energy minimization was then performed with restraints only on heavy nonwater atoms. The system was then equilibrated in the NVT ensemble at 300 K for 1 ns. The system was then sequentially simulated for 5 ns in the NPT ensemble using the Berendsen47 thermostat at 300 K and Berendsen barostat47 at 1 atm, where nonhydrogen atoms of the ribosome were given restraints of strength 1000, 100, 10, 1, and 0 kJ/mol. The system was then subjected to an additional 5 ns of equilibration under the NPT ensemble using the Nose-Hoover48,49 thermostat set at 300 K and the Parrinello-Rahman barostat50 with pressure coupling at 1 atm. During the final equilibration simulation, as in production simulations, position restraints were only included on residues at the boundary of the simulated system, as employed in previous studies.28,29,33
After equilibration of each system, we initialized 305 independent 15 ns simulations (3355 simulations, 50 μs aggregate sampling). The first 2 ns of each simulation were not analyzed to reduce potential bias of the initial configurations. The NPT ensemble was used for all production simulations.
B. All-atom structure-based model
We performed short-length scale molecular dynamics simulations of tRNA and the accommodation corridor, where the tRNA was initialized from various positions along the elbow accommodation process (Relbow = 27 + n, for 1 ≤ n ≤ 28). The structures were extracted from the structure-based simulation described in the work of Yang et al.29 Starting from each configuration, we performed 100 simulations with random initial velocities, for a total of 2900 simulations. All simulations employed an all-atom structure-based model,51 where force field files were generated using the smog-server.org tool.52 The classical A/A conformation of cognate Phe-tRNAPhe in complex with the ribosome (PDB:3I8F53) was used to define the global energetic minimum. Since the interactions in this model describe effective energetics, and tRNA-ribosome interactions form transiently, relative to intraribosome interactions, the strengths of the tRNA-ribosome interactions were reduced. For full force field details, see Ref. 29.
The simulations were performed using a modified version of Gromacs 4.6.7,43 which included source code modifications to introduce anisotropic position restraints. We simulated a subset of the ribosomal atoms with anisotropic position restraints on the boundary atoms. For details, see Ref. 29. The temperature was held constant at 0.42 (reduced units) using Langevin dynamics protocols (temperature coupling constant of 1 time unit). Each simulation was performed for 500 000 time steps of size 0.002. In addition, we repeated all diffusion calculations using long trajectories in which barrier-crossing events were spontaneously observed (data set described in Ref. 29). Both data sets yielded comparable values of D that exhibit minimal position-dependence (not shown).
C. Estimating D
In this study, we applied three distinct methods to estimate the diffusion coefficient as a function of tRNA position.
As described in the results, the first approach was to calculate the average displacement-squared along Relbow, only including initial frames that were within a ±0.5 Å range of values. A linear fit was then applied to estimate D according to Eq. (2). The fit was applied to Δt = 1–5 ns for the explicit-solvent model [Fig. 2(b)] and Δt = 1–4 reduced units for the structure-based model. The fits to explicit-solvent data were repeated for Δt = 1–8 ns, and no notable differences were observed (not shown). The calculation of D was performed for Relbow = 27 + n1 ± 0.5 Å (n = 0–29). All linear fits yielded linear correlation coefficients that were greater than 0.95.
The second approach was to calculate the variance of Relbow values, calculated over the same time and spatial intervals used for Eq. (2). However, in this approach, a linear fit was applied to the variance, in accordance with Eq. (4).
The third approach involved estimating a transition matrix for displacements along Relbow and from the rate matrix estimating the diffusion coefficient, as described previously.8 For completeness, we will summarize the approach. First, one discretizes the configuration space along a reaction coordinate ρ. For consistency with earlier calculations, we used Relbow as the coordinate and we discretized the dynamics in 1 Å intervals. Next, one defines a fixed time interval tα, over which transitions between cells are calculated. This yields a coarse-grained trajectory that is discretized in space and time. Finally, one defines a likelihood function
Consistent with Ref. 8, we also multiplied the likelihood function by a smoothing function, which favors more modest position-dependence of D,
where D(ρi) is the estimated values of D for cell i. D is estimated according to the relation
where Δρ is the cell spacing and Ri,j represents the transition matrix elements. Using these relations, one applies a Monte Carlo search in Ri,j, where moves are accepted if the likelihood function is increased. If the likelihood function is decreased (ΔL < 0), the move was accepted with probability exp(ΔL). After excluding initial “equilibration” moves, the average D was calculated for each value of the reaction coordinate. For analysis of the explicit-solvent simulations, we used tα = 3 ns and γ = 4 * 10−6μm2/s. We repeated the calculations for various tα (1–5 ns) and γ (10−6 to 10−5) values and found that the strong position-dependence was robust to the choice of parameters (not shown). For analysis of the structure-based model, we used tα = 4–32 and γ = 10−4. The lack of position-dependence was found to be robust to choice of parameters. However, the overall scale of D did exhibit a dependence on tα. tα = 10 is shown.
IV. CONCLUSION AND DISCUSSION
As the field continues to study the energy landscapes that govern biological dynamics, a comprehensive understanding of diffusion will be essential. The most pressing need for knowledge of the diffusion coefficient is that it enables one to interconvert between theoretical free-energy landscapes and experimental kinetics. By identifying the characteristics of the diffusion coefficient, here, we provide the first analysis of position-dependent diffusion along a suitably chosen reaction coordinate for the ribosome. Specifically, we calculated the diffusion along a coordinate that can precisely identify the large-scale free-energy barrier associated with tRNA-elbow accommodation.33 With this quantitative framework, it is now possible to explore any range of theoretical representations, where experimental rates may be used as a benchmark for determining whether a model appropriately describes the physical-chemical factors that govern kinetics.
In the context of a bacterial ribosome, the observed position-dependence of the diffusion coefficient begins to reveal the hierarchical character of the ribosomal energy landscape, where large-length-scale barriers are ∼10kBT and short-scale barriers are ∼1 to 2kBT. Together, these initial insights into the diffusive properties open up many opportunities for systematically investigating the relationships between structure, energetics, and dynamics of complex biomolecular assemblies. For example, this foundation may be applied to quantitatively investigate the effects of chemical perturbations, such as posttranscriptional modifications, deacylation, and small-molecule binding. Alternately, analysis of diffusive properties may be extended to describe the energetics and kinetics of other conformational rearrangements, including hybrid-state formation and translocation, or to study the ribosomes of higher-level organisms.
SUPPLEMENTARY MATERIAL
See supplementary material for time traces of simulated system and image depicting contact formation in the elbow-accommodated ensemble.
ACKNOWLEDGMENTS
We would like to thank the Clementi group (Rice University) for sharing Matlab scripts that were used for Bayesian Inference-based estimates of the diffusion coefficient. This work was supported by the National Science Foundation (Grant Nos. MCB-1350312 and MCB-1915843). R.H. was supported by NSF Grant No. HRD 16-19629 Northeast LSAMPAlliance: 2016-2021. V.B.P.L. and J.C. were supported by FAPESP Grant No. 2014/50739-5. V.B.P.L. was also supported by FAPESP Grant Nos. 2018/18668-1 and 2016/19766-1. We acknowledge the support from the Northeastern University Discovery cluster. The authors declare no competing financial interest.
REFERENCES
The bounds on the inner integral assumes that ρA/T > ρEA.
The distance between the U8 residue of the P-site tRNA and U47 of the aa-tRNA.