All-atom simulations of crowding effects on ubiquitin dynamics

Simulations were run with the NAMD [33] code using the amber99SB [34] force field for the protein, TIP3P [35] parameters for water and Glycam06 [32] parameters for α-D-glucose. The structure of human ubiquitin in PDB ID 1D3Z [36] (model 1) was used as a starting point for minimization, equilibration and production in both simulations. For the simulation of ubiquitin in water, the molecule was solvated in a TIP3P box (315 396 ų); whereas for the simulation under crowded conditions, the molecule was first added to 317 glucose molecules using the Packmol [37] program, and then solvated with TIP3P water molecules. After minimization and equilibration the volume of the box was 291 135 ų, resulting in an average concentration of 325 g L⁻¹.

To build the reference frame for trajectory projection, principal components analysis (PCA) was carried out on the covariance matrix of Cα positions for residues 2–70 of 72 high-resolution x-ray structures of human ubiquitin aligned to model 1 of PDB ID 1D3Z [36]. This is the same strategy used by others [22, 27], but is based on a larger set of structures. This set includes free monomers, covalent oligomers and complexes with binding partners, all solved at a resolution equal or better than 2.5 Å (the list of PDB codes and PDB titles can be found in table S1 (available from stacks.iop.org/PhysBio/10/045006/mmedia)). The projection of the trajectories on the reference frame was also preceded by an alignment to model 1 of PDB ID 1D3Z. The sum of squared PCA loadings (i.e. eigenvalues) 1 and 2 shows that the structural variability that spreads the x-ray structures on the reference frame arises from sequence segments 6–11, 33–36 and 46–49 (figure S2B (available from stacks.iop.org/PhysBio/10/045006/mmedia)), which indeed show different conformations in the set of structures (root mean squared fluctuations (RMSF) from structures in figure 1(A)).

Zoom In Zoom Out Reset image size

Free-energy landscapes were built from the trajectories projected on the reference frame by binning the projection into a grid and counting the number of frames inside each cell of the grid (N_i) to compute its free energy according to:

$$\begin{equation*} G_i = - k_B T\ln ( {N_i /N_0 }) \end{equation*}$$

where N₀ corresponds to the most populated bin, which thus sets the free energy offset.

The average number of ubiquitin–partner interactions per residue was approximated as the average (through all x-ray structures that contained a protein interacting with ubiquitin) number of Cα atoms of ubiquitin-bound proteins that are within 5 Å of each Cα atom of ubiquitin.

N–H order parameters (s²) and residual dipolar couplings (RDC) were computed from the trajectories using equations reported in the literature. The trajectory is aligned to model 1 of PDB ID 1D3Z and the normalized N–H vectors u are retrieved for each peptide bond in each frame of the aligned trajectory. The calculation of s² is done as described elsewhere [38] assuming full decay of the autocorrelation function for u during the time of the trajectory:

$$\begin{equation*} s_i^2 = \frac{1}{2}\bigg( {3\sum\limits_{\alpha = x,y,z} {\sum\limits_{\beta = x,y,z} {\langle {u_{i,\alpha } u_{i,\beta } } \rangle ^2 } } - 1} \bigg) \end{equation*}$$

where the average 〈u_i,α u_i,β〉 is computed over the whole trajectory.

The N–H RDC were computed for each alignment medium for which data are available, as described by Showalter et al [39]. Briefly, for each medium an alignment vector A = (A_xx, A_yy, A_xy, A_xz, A_yz) is fit that minimizes the difference between the experimental RDC data (vector D_exp with as many elements as N–H peptide bonds for which data are available in that alignment medium) and the back-calculated RDC values (D_back of same size) in a least-squares sense. The back-calculated data is D_exp = M × A, where M is a matrix derived from the trajectory formed by one row of 〈u_x²〉 − 〈u_y²〉, 〈u_y²〉 − 〈u_z²〉, 2〈u_xu_y〉, 2〈u_yu_z〉 and 2〈u_yu_z〉 values for every N–H pair for which data are available. The least-squares problem (M × A – D_exp → min) is solved by singular value decomposition as described in detail in that same work [39].

Interactions between water or glucose molecules and the protein were investigated by searching for contacts between the small molecules and the protein. Two atoms were considered to be in contact when the distance between them was shorter than 1.1 times the sum of their Van der Waals radii, as optimized in recent works [40, 41].

RMSF were computed for each Cα atom as:

$$\begin{equation*} {\mathop{\rm RMSF}\nolimits} _i = \sqrt {\frac{{\sum_j {( {r_{i,j} - \overline {r_i } } )^2 } }}{N}} \end{equation*}$$

where r_i denotes the position of the Cα atom i, and j runs through the N frames of simulated time or x-ray structures. Cα atoms 2 to 70 were considered for the ensemble of x-ray structures and Cα atoms 1 to 76 for the simulations. Calculations were done after alignment to a reference structure (i.e., model 1 of 1D3Z).

The analysis proposed by Hess [42] was carried out on the basins observed in the conformational landscape of each simulation to describe them in terms of their deepness and roughness. First, ∼20 ns long sections of the trajectories were identified in which the protein remained inside each basin according to their projections on the two-dimensional reference frame built from the x-ray structures. For each sub-trajectory, PCA was carried out to obtain the most important fluctuations inside each basin; the first and second eigenvalues (λ) were used to compute the force constants k that defines the harmonic well in the direction of the two most important eigenvectors, using the formula:

$$\begin{equation*} k = \frac{{k_B T}}{\lambda } \end{equation*}$$

where k_B is Boltzmann's constant and T is 300 K. The autocorrelation functions were computed for the first principal components and their first decaying parts were fitted to obtain the reported τ parameters. The internal friction-like coefficient η is then:

$$\begin{equation*} \eta = k\tau \end{equation*}$$

which measures the roughness of the basin completing its description together with the harmonic well constant k.

The dynamic features of ubiquitin in dilute solution have been vastly studied by NMR, revealing nano- to micro-second timescale motions that are important for target recognition [20, 22, 27–29]. These experimental observations were complemented with molecular dynamics simulations in several works to further dissect the contributions of conformational selection and induced fit mechanisms to target binding [27, 28, 43]. Of particular importance, 0.5 to 1 µs long simulations performed with amber99SB, a highly accurate force field that fairly reproduces NMR dynamic data, showed that free ubiquitin is constantly exploring conformational states that resemble those observed in x-ray structures of the protein complexed with different partners [27, 28]. Based on these works, we first performed a reference 500 ns long simulation of free ubiquitin in water using the amber99SB force field. As reported for previous simulations, ours also reproduces fairly well the experimental order parameters available from relaxation and RDC data [14] (figure S1A (available from stacks.iop.org/PhysBio/10/045006/mmedia)) and the RDC measured in 39 alignment media [17, 18, 25, 44–46] (figure S1B (available from stacks.iop.org/PhysBio/10/045006/mmedia)). The RMSF profile derived from the simulation (black curve in figures 1(C) and pictorial representation (E)) follows that observed in a set of x-ray structures of ubiquitin bound to different partners (figure 1(A)), both pointing at increased dynamics of the loops involved in the recognition of ubiquitin's binding partners (figures 1(B) and (D)).

All in all, this analysis shows that the amount of flexibility observed in the trajectory is consistent with experimental data in dilute solution and that flexible regions of the protein are important for its ability to bind different partners. However, the analysis does not reveal any information about the conformational space accessible to the protein through collective motions as a result of such flexibility. In order to explore this we projected our simulation on a two-dimensional reference frame built by mapping through principal components analysis the structural variability observed in 72 experimental x-ray structures of human ubiquitin under different scenarios (figure S2A (available from stacks.iop.org/PhysBio/10/045006/mmedia), where each red point is an x-ray structure). This method has recently been introduced to analyze simulation trajectories and ensembles for which no specific reaction coordinate can be defined a priori [22, 27, 28]. Briefly, this method combines the collective fluctuations of the most variable regions of the protein into two new variables, i.e. the principal components 1 and 2, which account for as much fluctuation as possible. Projection of the Cα coordinates of our trajectory on this two-dimensional frame (figure 2(A) in blue, where only one frame was plotted every 1 ns for simplicity) shows that the protein explores the whole range of conformations sampled by the x-ray structures, covering smoothly the conformational space. Figure 2(B) is a simpler and more informative representation of the projection, in which the density of frames projected on a grid was converted into a free energy landscape relative to the deepest well. The largest barriers separating these states are around 1 kcal/mol, indicating a nearly flat energy surface consistent with the very fast exchange observed between conformational states and with the fast dynamics determined by NMR.

Zoom In Zoom Out Reset image size

In the next step we investigated a 500 ns long MD trajectory of ubiquitin in a solution crowded at 325 g L⁻¹ of glucose. For this simulation we combined the amber99SB force field [34] to describe the protein and the Glycam06 [32] force field to describe glucose molecules. Our setup leaves a random distribution of water and glucose molecules in a ratio of about 20:1 around the protein molecule; details on the system setup and simulation are given under methods.

A comparison of the RMSF in 325 g L⁻¹ glucose against the values computed in water (red versus black in figure 1(C), respectively) shows that the amplitudes of loop motions are slightly reduced under crowded conditions during the timescale of the simulations. Notably, the effect on RMSF is not as drastic as that observed in a previous simulation of myoglobin under much more crowded, near-dry conditions at almost 90% w/w sucrose [47, 48]. As a result of the mild reduction in loop flexibility, the RMSF and s² values of the loops homogenize in the crowded condition, indicating similar mobility at 325 g L⁻¹ crowder concentration. Importantly, and as shown in figure S4 (available from stacks.iop.org/PhysBio/10/045006/mmedia), this equalizing effect is also seen on the N–H order parameters, which can be obtained experimentally and compared to test our findings. One final point worth noticing is that the root mean squared velocities of Cα atoms are roughly similar in both simulations (data not shown) indicating that the effect of crowding is on the amplitudes of local atomic fluctuations but not on the atoms' speeds.

In order to explore the impact of crowding on collective motions and the conformational landscape explored by the protein, we projected the simulation on the two-dimensional reference frame built above. Projection of one frame every 1 ns (figure 2(C)) shows that the protein has not explored the same conformational space as in water during the same time length. Instead, it is restricted to small excursions slightly away from its initial conformation, always staying within the same basin. This is also evident in the time evolution of the conformational space sampled in figure 3. The finding of a strong restriction in the size of the explored conformational space with only a minor reduction in RMSF indicates that local fluctuations take place to a similar extent under both conditions but that collective motions (i.e. transitions between basins) have been compromised or more likely slowed down in the crowded condition.

Zoom In Zoom Out Reset image size

In only 40 ns of simulation time the protein simulated in water has escaped from the initial energy well and has already explored most of the conformational space that it will sample during the rest of the simulation (figure 3, left). In contrast, the simulation in 325 g L⁻¹ glucose is still trapped around its starting point at 300 ns (figure 3, right) and has made only one unsuccessful attempt to visit an alternative conformation at around 400 ns (figures 2 and S3 (available from stacks.iop.org/PhysBio/10/045006/mmedia)). We also observe that the speed at which the protein explores the conformational space is slower in the crowded condition than in water. In order to quantify this, we studied the rate at which the conformational basins are explored based on a formalism proposed by Hess [42] to model the 'diffusion' of a protein inside a conformational basin in terms of the internal friction coefficients (η) and harmonic force constants (k) characteristic of the basin. We analyzed 20 ns long sections of the trajectories in which the protein remained inside each of the two basins of the simulation in water and in the single basin observed in glucose solution. The parameters obtained for the two first principal components of motion inside each basin are given in table 1. The two basins observed in water are characterized by similar harmonic force constants, whereas the basin observed in the crowded condition features a ∼5 times stronger k meaning a deeper and sharper energy well. The later is also rougher with η being an order of magnitude larger, meaning an increased 'friction' and a much slower exploration of the basin in the presence of crowder.

In brief, our results predict that crowding slows down the internal dynamics by exerting an increased 'internal friction' and possibly by imposing higher energy barriers around the conformational basins thereby reducing the rate of exchange between them. Importantly, the effect is important for large scale motions but not that much for atomic fluctuations. Similar ideas implying that crowding slows down internal motions have been suggested by experimental observations in the last few years [49, 50]. Also, recent computational and experimental studies showed that the confinement of a protein inside a micelle reduces its internal dynamics due to interactions with the surfactant [51, 52]. With all this evidence, the notion that internal protein dynamics are affected by crowding starts to gain support.

Considering the high concentration of glucose molecules in the simulation box, encounters between these molecules and the protein are very likely to occur. Indeed, visual inspection reveals extensive associations of several glucose molecules with the protein, resembling the multimolecular complexes observed in a recent all-atom simulation of the cytoplasm of E. coli [3]. In the absence of crowders, the protein interacts with an average of 289.6 water molecules at any given moment; whereas in 325 g L⁻¹ glucose, this decreases to an average of 120 water molecules due to an average of 73.1 glucose molecules that attach to the protein surface. This means that every ∼169.6 water molecules that are excluded from the surface, ∼73.1 glucose molecules are attached, resulting in a ratio of ∼2.3 water/glucose present at the protein surface. This ratio is consistent with the larger area of glucose and shows how interactions between the crowder and the protein result in desolvation of its surface, as observed through simulations and experiments for near-dry conditions [47, 48, 53] but to a lesser extent. Notably, five glucose molecules never leave the protein surface along the whole trajectory. This does not imply constant interaction of the glucose molecules with the same protein residues; rather, the glucose molecules shuffle interactions to different atoms of the protein wandering around its surface. For example, there is a glucose molecule that interacts with the first loop of the protein, switching hydrogen bonds between its hanging alcohol groups and different atoms of the protein during the whole simulation (figures 4(D) and (E)).

Zoom In Zoom Out Reset image size

Analyses per residue (figures 4(A)–(C)) reveal homogeneous desolvation and glucose binding along the exposed parts of the protein, suggesting extensive enthalpic perturbations that average smoothly on the protein surface. This picture is in line with the recent finding that unspecific interactions are an important component driving crowding effects on top of the already known steric effects. Regarding desolvation, perturbations of the hydration structure and dynamics were predicted in a recent simulation meant to assess the effect of protein rather than small-molecule crowders [10]. In the case of glucose molecules, however, the perturbation of hydration structure and dynamics on the first water shell seems stronger, which can be attributed to the much smaller size of glucose and its higher density of polar groups. In our interpretation, glucose molecules bound to the protein and water molecules trapped in-between drag the exposed residues. Since the interactions are not specific, the effect is generalized and affects all collective motions making them less likely to happen; in other words, slowing them down. Further work is underway to better understand these effects of crowding on the internal protein dynamics, and to evaluate the role of protein–glucose, protein–water and glucose–water interactions.

Although the constituents of the cellular medium are chemically varied and not simply glucose rings, they contain exposed polar and charged groups capable of establishing electrostatic interactions, salt bridges and hydrogen bonds with protein surfaces, which are largely polar and charged. In our view, based on these and others' results, all these transient, unspecific interactions average out to a constant enthalpic and entropic perturbation of the protein surface. This perturbation would imply both a thermodynamic stabilization of the protein fold and a slow down in all large-scale kinetic processes thus slowing down unfolding rates too. In this regard, our observations are in agreement with the experimental report that crowding enhances the native structure of proteins [54, 55] and that part of this stabilizing mechanism acts on the folding pathways themselves [12]. On top of this, our finding that collective motions are slowed down allows us to propose that crowding agents can stabilize proteins by simply trapping the native conformations, i.e. through a kinetic contribution. Similar ideas were suggested by different experimental and computational evidence [56, 57]; thus this possibility deserves further exploration.

Notably, these effects would also stabilize protein complexes, as well as complexes composed of different kinds of macromolecules as reported in a recent all-atom simulation of E. coli's cytosol. This might be of little importance for strongly interacting proteins but would be important for proteins involved in interactions defined as 'transient' by experiments. In particular this can explain why some interactions predicted by genetic evidence are not confirmed by in vitro experiments. Stabilization of protein–protein interactions by crowders would have profound implications for many proteins, including ubiquitin itself. Moreover, crowding could have an effect on the relative weight of induced-fit and conformational selection mechanisms that guide protein–protein recognition.

Finally, we can speculate that the marginal stability observed for proteins in solution arises from the fact that they have evolved to work in crowded, stabilizing environments. Under this scenario, potential proteins that are too stable in dilute solution would be unnecessarily stable in the cellular medium and possibly too rigid to function properly.