Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3050011

Formats

Article sections

Authors

Related links

Phys Chem Chem Phys. Author manuscript; available in PMC 2011 March 8.

Published in final edited form as:

Phys Chem Chem Phys. 2010 April 14; 12(14): 3491–3500.

PMCID: PMC3050011

NIHMSID: NIHMS273404

See other articles in PMC that cite the published article.

A cell's interior is comprised of macromolecules that can occupy up to 40% of its available volume. Such crowded environments can influence the stability of proteins and their rates of reaction. Using discrete molecular dynamics simulations, we investigate how both the size and number of neighboring crowding reagents affect the thermodynamic and folding properties of structurally diverse proteins. We find that crowding induces higher compaction of proteins. We also find that folding becomes less cooperative with the introduction of crowders into the system. The crowders may induce alternative non-native protein conformations, thus creating barriers for protein folding in highly crowded media.

The presence of macromolecules at high concentrations within the interior of a cell creates a crowded environment in which a single macromolecule must adapt and function near regions of excluded volume.^{1} Limitations of the available volume induce conformational changes in the polypeptide to fold into a compact shape. This phenomenon has attracted both experimental and theoretical studies that have quantified the effects of macromolecular crowding on the increased stability of a protein's native state.^{2}^{,}^{3}

Various other biochemical equilibria are affected by the excluded volume environment, including oligomerization and chemical reactivity. Dimerization of protein–protein complexes may favor either the association or dissociation, depending on which final state excludes less volume.^{4} The rates of chemical reactions can vary drastically depending on the speed of the reaction measured in an ideal solution. Reactions that occur quickly at ideality become increasingly slower at higher concentrations of crowders as the rate becomes diffusion-limited. Alternatively, reactions that are normally slow under dilute conditions can be expedited by the presence of crowders as the reactants have a higher probability of re-colliding at the appropriate geometry due to caging effects that keep the reactants in close proximity.^{5}

Macromolecular confinement describes a similar phenomenon in which a protein or other macromolecule is enclosed within a region of limited volume.^{6} The GroEL-GroES chaperonin system found in *E. coli* is an important biological example that serves to increase the rate of folding and prevent off-pathway protein aggregation.^{7} Protein folding under confining conditions has been studied analytically^{8} and with simulations^{9} using a variety of techniques including Gō-like models^{10} and all-atom potentials.^{11}^{,}^{12} Confining environments can influence the directionality of protein folding by compacting the unfolded state. At extremely limited volumes, proteins can experience downhill-like folding due to destabilization of the transition state.^{10} Interestingly, other studies have shown that the loss of solvent entropy in some proteins actually destabilizes the secondary structure and favors the unfolded state.^{12}^{,}^{13} Overall, the effect of confinement differs depending on the macromolecule studied.

Theoretical treatments of macromolecular crowding have relied mostly on statistical thermodynamic models, which approximate the crowding effect as a potential of mean force acting to compress a polypeptide chain into its native state.^{2} The use of simulations becomes advantageous since crowding reagents can be modeled explicitly. Several studies have employed a variety of simulation methodologies including density functional theory,^{14} Brownian dynamics,^{15}^{,}^{16} and molecular dynamics.^{17} Due to the size and complexity of proteins, most models used in simulations have been coarse-grained.

One of the main limitations with molecular dynamics simulations is its computational cost. In a recent all-atom study of crowding Qin and Zhou calculated the chemical potential for inserting a 106 residue protein inside a configuration of spheres using a novel algorithm to determine the allowed fractional volume allocated to the protein.^{17} However, configurations generated for the protein itself were limited to 30 ns. For larger proteins that may contain multiple domains, an adequate sampling of protein conformational space becomes a limiting factor.

Discrete molecular dynamics (DMD) offers a computationally efficient alternative for simulating biological systems by utilizing a combination of simplified models in conjunction with discretized potentials. DMD differs from its traditional counterparts such that simulations proceed according to ballistic equations of motion. Since each collision is sorted via an efficient search algorithm, simulations scale on the order of *N* ln *N* as opposed to *N*^{3} with traditional molecular dynamics.^{18} Given the discontinuous nature of DMD simulations, interactions between pairs of particles are described with discretized potentials. Each particle present within a simulation may represent an entire amino acid residue, regions of a residue, or a single atom. Therefore, the level of coarse-graining may enhance the speed of calculations although at the expense of accuracy. The choice of a model is dependent on which biological timescale is relevant to the study at hand. Several models of DMD are available and include applications in protein folding,^{19} and aggregation,^{20} protein–DNA complexes,^{21} RNA folding,^{22} and lipids.^{23} In comparison to traditional molecular dynamics, there can be a speed enhancement on the order of 10^{8}–10^{10} depending on the model used to describe the system.^{18} For an all-atom model, DMD still remains a slightly faster alternative due to the discretized interactions between atom pairs.^{19} However, as the potentials become more continuous (and, presumably, more accurate), the computational cost becomes equivalent to traditional molecular dynamics.

Here we introduce the first implementation of DMD to study macromolecular crowding with the incorporation of residue-specific information for each protein by employing replica-exchange simulations with the DMD algorithm to rapidly sample crowder and protein conformational space. The advantage of this approach is that we can observe how crowders directly interact with the protein to influence folding, as opposed to calculating the chemical potential of stabilization. Thus, by explicitly simulating the interactions between crowders and proteins at a residue-level, it is possible to study the folding pathway instead of only the initial and final states. We generate thermodynamic profiles for several proteins that give insight on how crowding effects scale with respect to sequence length. The thermodynamic profiles also provide a method for viewing crowding effects on the folding pathway enabling new mechanistic insights to emerge. Our results suggest that while crowding imposes compaction of a protein, the interactions necessary for cooperative folding become disrupted.

Traditional molecular dynamics simulate the motions of particles by solving Newton's equations of motion for a defined system using an integration algorithm. In DMD, simulations proceed according to the conservation laws of energy, momentum, and angular momentum and are evaluated as a series of two-body interactions.^{24} The efficiency of the engine is based on an algorithm that searches through an event table, where velocities are only modified as necessary. Here we classify an event as the instance in which two particles are within a defined interaction range as defined by their potential. The potentials used in DMD are discretized to accommodate the discontinuous nature of the simulations. Further details of the DMD algorithm can be found elsewhere.^{25}^{,}^{26} Here we extend its applications toward the modeling of crowding effects.

Let us consider a system of particles, including two particles *i* and *j* both of mass *m* that occupy initial positions of *r*_{i0} and *r*_{j0} with initial velocities of *v _{i}* and

$${U}_{\mathit{ij}}=\{\begin{array}{cc}\infty ,\hfill & \phantom{\rule{1em}{0ex}}r<{\sigma}_{1}\hfill \\ -\epsilon ,\hfill & \phantom{\rule{1em}{0ex}}{\sigma}_{1}<r<{\sigma}_{2}\hfill \\ 0,\hfill & \phantom{\rule{1em}{0ex}}r>{\sigma}_{2}\hfill \end{array}\phantom{\}}$$

The variable σ_{α} defines an interaction distance, where α = 1 refers to the hard sphere repulsion and α = 2 is the attractive interaction. For a system with *N* particles, the potential energy of the system is

$$E=\frac{1}{2}\sum _{i,j=1}^{N}{U}_{i,j}$$

During a simulation, a table is generated by calculating event times for each pair of particles. The time interval in which an event may occur between two particles is

$${t}_{\mathit{ij}}^{\left(\alpha \right)}=\frac{-{b}_{\mathit{ij}}\pm {[{b}_{\mathit{ij}}^{2}-{v}_{\mathit{ij}}^{2}({r}_{\mathit{ij}}^{2}-{\sigma}_{\alpha}^{2})]}^{1\u22152}}{{v}_{\mathit{ij}}^{2}}$$

where

$${\mathbf{r}}_{\mathit{ij}}={\mathbf{r}}_{i0}-{\mathbf{r}}_{j0}$$

$${\mathbf{v}}_{\mathit{ij}}={\mathbf{v}}_{i}-{\mathbf{v}}_{j}$$

$${b}_{\mathit{ij}}={\mathbf{r}}_{ij}\cdot {\mathbf{v}}_{ij}$$

and the plus or minus sign refers to two particles either approaching or receding from each other, respectively. The trajectory of particle *i* is evaluated with respect to time as

$${\mathbf{r}}_{i}(t+{t}_{u})={\mathbf{r}}_{i}\left(t\right)+{\mathbf{v}}_{i}\left(t\right){t}_{u}$$

where the time unit *t _{u}* is determined by comparing the shortest event time to the maximum allowed time interval

When particles *i* and *j* are determined to interact at ${t}_{\mathit{ij}}^{\left(\alpha \right)}$ the squared difference in the particles' positions is compared to the squared interaction distance prior to a change in the potential (i.e., before ${r}_{\mathit{ij}}^{2}-{\sigma}_{\alpha}^{2}=0$). There are three possible scenarios for our one-step potential as the two particles approach one another. During an attractive encounter, if ${r}_{\mathit{ij}}^{2}>{\sigma}_{2}^{2}$ then the magnitude of separation between the particles decreases and the change in velocities will be

$$\Delta {\mathbf{v}}_{i}=-\Delta {\mathbf{v}}_{j}=\frac{-{\mathbf{r}}_{\mathit{ij}}}{2{\sigma}_{2}^{2}}\left[{b}_{\mathit{ij}}+{\left({b}_{\mathit{ij}}^{2}-\frac{4{\sigma}_{2}^{2}\epsilon}{m}\right)}^{1\u22152}\right]$$

However if ${r}_{\mathit{ij}}^{2}<{\sigma}_{2}^{2}$, then the two particles increase their magnitude of separation and the change in velocities will be

$$\Delta {\mathbf{v}}_{i}=-\Delta {\mathbf{v}}_{j}=\frac{-{\mathbf{r}}_{\mathit{ij}}}{2{\sigma}_{2}^{2}}\left[{b}_{\mathit{ij}}-{\left({b}_{\mathit{ij}}^{2}+\frac{4{\sigma}_{2}^{2}\epsilon}{m}\right)}^{1\u22152}\right]$$

In the case of hard sphere repulsions, the change in velocities is simply

$$\Delta {\mathbf{v}}_{i}=-\Delta {\mathbf{v}}_{j}=\frac{-{\mathbf{r}}_{\mathit{ij}}{b}_{\mathit{ij}}}{2{\sigma}_{1}^{2}}$$

After each event, we remove from the table events that were calculated with the previous two particles since their velocities and positions have changed. The number of possible events with these two particles are then recalculated and sorted within the table. Simulations then proceed through the table, where particles are allowed to move between time intervals, until either ${t}_{\mathit{ij}}^{\left(\alpha \right)}<{t}_{m}$ or the table of events is depleted.

Thermodynamic comparisons of a protein in the absence or presence of crowders require adequate sampling of each system. During a simulation, various configurations of the system lead to local minima within its potential energy landscape. Higher crowder concentrations further frustrate the energy landscape between different minima by inhibiting the protein from exploring spatially larger conformations. While increasing the temperature is an effective way to cross such barriers, entropic contributions to the free energy begin to dominate over the enthalpic effects of interest. We utilize replica-exchange discrete molecular dynamics (REXDMD) as an efficient alternative to exploring the conformational space of each protein.^{19}^{,}^{27}

In REXDMD, there are *R* replicas (copies) of a given simulation. Each replica contains the system at a unique temperature *T _{i}*, where

$$p=\{\begin{array}{cc}\hfill 1,\hfill & \hfill \text{for}\phantom{\rule{1em}{0ex}}\lambda \le 0\hfill \\ \hfill \text{exp}(-\lambda ),\hfill & \hfill \text{for}\phantom{\rule{1em}{0ex}}\lambda >0\hfill \end{array}\phantom{\}}$$

where

$$\lambda =\left(\frac{1}{{k}_{\mathrm{B}}{T}_{i}}-\frac{1}{{k}_{\mathrm{B}}{T}_{j}}\right)({E}_{j}-{E}_{i})$$

and *E* is the total energy of the replica system.

We analyze the trajectories obtained from replica-exchange simulations *via* the weighted histogram analysis method (WHAM) provided by the MMTSB tool.^{28} This methodology generates thermodynamic profiles of a system by combining data obtained through multiple biased simulation trajectories. A detailed derivation describing the method^{29} and its application to replica-exchange simulations^{30} have been described previously in the literature. Here we give a brief overview of applying WHAM for replica-exchange analysis.

Given a set of *i* (where *i* = 1,…,*R*) replica-exchange simulations with *s* snapshots, we begin constructing a histogram by taking the reaction coordinate *A* and the potential energy *E* and placing them into the bins *x*_{α} (where α =1,…,Θ) and *E*_{β} (where β =1,…,Ξ) respectively. The count of the αth *x* bin and βth *E* bin within the *i*th simulation is denoted as *n _{i}*

$${p}_{i(\alpha ,\beta )}={f}_{i}{\lambda}_{i(a,\beta )}{\stackrel{\u02da}{p}}_{\alpha ,\beta}$$

where the biasing factor is

$${\lambda}_{i(\alpha ,\beta )}=\mathrm{exp}\left[-\frac{{E}_{k}}{{k}_{B}}\left(\frac{1}{{T}_{i}}-\frac{1}{{T}_{0}}\right)\right]$$

and *f _{i}* is the normalization factor such that

$$\sum _{\alpha ,\beta}{\stackrel{\u02da}{p}}_{\alpha ,\beta}=1$$

Here the unbiased probability is represented by ${\stackrel{\u02da}{p}}_{\alpha ,\beta}$ We extend our discussion of a single histogram to the generalization for WHAM, where the conditional probability observing a given value of *A* at a given energy *E* is

$$p(A\mid E)=\frac{\sum _{i=1}^{R}\sum _{\alpha =1}^{\Theta}\sum _{\beta =1}^{\Xi}{n}_{i(\alpha ,\beta )}\mathrm{exp}\left(-\frac{{E}_{i}}{{k}_{\mathrm{B}}T}\right)}{\sum _{i=1}^{R}{s}_{i}\mathrm{exp}\left({f}_{i}-\frac{{E}_{i}}{{k}_{\mathrm{B}}T}\right)}$$

The generalized density of states from the *i*th simulation is computed

$${\Omega}_{i}(A,E)={n}_{i}(A,E)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left(\frac{{E}_{i}}{{k}_{\mathrm{B}}T}-{f}_{i}\right)$$

and the total density of states is

$$\Omega (A,E)=\sum _{i=1}^{R}{\omega}_{i}\left(E\right){\Omega}_{i}(A,E)$$

The set of weights {*w*} are determined through minimization of the statistical error for the best estimate of Ω_{i}(*A*,*E*).

Once *P*(*A*|*E*) and Ω_{i}(*A*,*E*) are determined, we can calculate the potential of mean force (PMF) as a function of reaction coordinate *Avia*

$$PMF\left(A\right)=-\mathrm{ln}\left(\int P(A\mid E)\Omega (A,E)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left(-\frac{E}{{k}_{\mathrm{B}}T}\right)\mathrm{d}E\right)+C$$

where *C* is a constant in which the lowest *PMF* is normalized to zero. With the density of states known, it is possible to calculate the specific heat of folding with the partition function

$$Z=\int \Omega \left(E\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left(-\frac{E}{{k}_{\mathrm{B}}T}\right)\mathrm{d}E$$

From thermodynamic principles, it follows that the variance in energy can be computed from the partition function as

$$\delta E=\frac{{\partial}^{2}\mathrm{ln}Z}{\partial {\beta}^{2}}$$

where β is 1/*k*_{B}*T*. Thus the specific heat of folding can be evaluated using the following relation

$${C}_{V}=\frac{{\delta}^{2}E}{{T}^{2}}$$

The radius of gyration for a given protein is evaluated as

$${R}_{g}^{2}=\frac{\sum _{i=1}^{N}{m}_{i}{r}_{i}^{2}}{\sum _{i=1}^{N}{m}_{i}}$$

where *m _{i}* is the individual mass of a bead and

Cooperative processes such as protein folding contain a peak in the heat capacity ${C}_{P}^{\mathrm{Max}}$ at the transition point and occurs at the folding transition temperature *T _{f}*. Although the heat capacity at constant pressure is generally derived from the enthalpy as opposed to the total energy, here we follow the approximation that

$$\Delta {H}_{vH}=2{\left({k}_{\mathrm{B}}{T}_{f}^{2}{C}_{P}^{Max}\right)}^{1\u22152}$$

We can measure the degree of cooperativity for protein folding as simply the ratio between the van't Hoff enthalpy to the enthalpy derived from calorimetry experiments^{32}

$${\kappa}_{2}=\frac{\Delta {H}_{vH}}{\Delta {H}_{\mathrm{cal}}}$$

We estimate the calorimetric enthalpy by evaluating the integral across the transition region^{32}^{,}^{34}

$$\Delta {H}_{\mathrm{cal}}=\int {C}_{P}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}T$$

Since WHAM enables determination of the partition function using the density of states, the thermodynamic functions and properties derived from the method are smoothed.^{29} For our analysis we use a total of 8 × 10^{6} trajectories (per REXDMD simulation) to ensure an accurate estimation of the density of states. We estimate the errors from the calculated radii of gyration by considering its fluctuations throughout a simulation at a given temperature. However, the calculation of standard errors for the specific heat of folding is a non-trivial task since it was not included within the original WHAM publication.^{29}^{,}^{35} Nonetheless, we have estimated the errors from the specific heat of folding by performing multiple WHAM calculations using independently derived trajectories from the simulations. Small estimates of error indicate adequate sampling of the system to study its thermodynamic properties.

We perform simulations on four different proteins whose choice was based on the number of amino acids present on the peptide chain: Trp-cage, Villin headpiece subdomain, Src SH3 domain, and hemoglobin α-chain (PDB codes 1L2Y, 1YRF, 1NLO, and 1FN3, respectively) (Fig. 1). Proteins were represented by a two-bead Gō model, where the two-beads corresponded to the *C*-α and sidechain of each residue (Fig. 2A).^{36} The Gō energy of a system is defined as the sum of Gō interaction energies for all residue pairs excluding *N* near-neighbor residue interactions:

$${E}^{G\stackrel{\u2012}{o}}=\sum _{I}\sum _{J\ge I+N}{E}_{IJ}^{G\stackrel{\u2012}{o}}$$

The Gō interaction energy between two individual residues is given as

$${E}_{IJ}^{G\stackrel{\u2012}{o}}=\in {\Delta}_{IJ}^{\mathrm{Nat}}$$

where

$${\Delta}_{IJ}^{\mathrm{Nat}}=\{\begin{array}{cc}1,\hfill & \mid {r}_{i}^{NS}-{r}_{j}^{NS}\mid >\sigma \hfill \\ -1,\hfill & \mid {r}_{i}^{NS}-{r}_{j}^{NS}\mid \le \sigma \hfill \end{array}\phantom{\}}$$

Thus if the distance between two residues in the native state reside within a predetermined cutoff, σ, then an attractive potential is assigned. For our Gō model we defined our parameters as σ = 7.5 Å and *N* = 3 (Fig. 2B). Crowders were modeled as hard spheres with no attractive interactions in all simulations.

Proteins simulated under crowded conditions. These proteins are comprised of a variety of topologies and sequence lengths. (A) Trp-cage, (B) villin subdomain, (C) SH3 domain, (D) hemoglobin alpha chain.

Coarse-grained modeling of proteins. (A) Proteins are represented as beads on a string with each residue comprised of two beads, Cα and Cβ. Bonds are denoted by solid lines and non-bonding interactions are denoted by dashed lines. (B) **...**

Three different conditions were the focus of our simulations: the crowding effect with respect to increasing the number of crowders (*N*), the crowding effect with respect to increasing the size of crowders (*r _{c}*) at constant fractional volume (

Crowders were represented as 10 Å particles within this subset of simulations. We defined the number of crowders inserted within a given simulation as *N* = 0, 100, 250, 500, 750, or 1000 particles. Since the volume of the box is 10^{6} Å^{3}, the corresponding fractional volumes of crowder are 0%, 5%, 13%, 26%, 39%, and 52%.

In this set of simulations, we define *N* = 250 and vary the crowder diameter as 3, 5, 7, 10, or 15 Å. The respective fractional volumes of crowder are 0.35%, 1.6%, 4.45%, 13%, and 44%.

The diameter of the crowder was varied in these set of simulations, where *d _{c}* was set to 5, 10, 15, or 20 Å. Crowders were inserted into the box until

We model the protein Trp-cage (*PDB* code 1L2Y, Fig. 2) using an all-atom model that is previously published.^{19} Interactions between atoms in this model are described using a discretized version of the CHARMM force field, as well as new parameters accounting for explicit hydrogen bond formation. All heavy atoms and polar hydrogens are included in the simulation. Crowders are defined as 10 Å particles that undergo hardcore collisions with all atoms. Each atom's hard sphere radius cutoff with the crowder is approximated to be 1.5 Å. Similar to the coarse-grained simulations, we vary the concentration of crowders by assigning N to be equal to 0, 500, 750, or 1000 particles. The simulation box volume is 10^{6} Å^{3} thus the corresponding fractional volumes occupied by crowders are 0%, 26%, 39%, and 52%.

All simulations start from the native conformation for each protein. After random insertion of crowders into a system, we initially relax the system to relieve any potential steric clashes by assigning temporary soft potentials to the crowders. The system then undergoes replica-exchange DMD for 10^{6} time units. We utilize 8 replicas where we increment the temperature by 0.1 ε/*k*_{B} such that reduced temperatures range from 0.5–1.2 ε/*k*_{B}. The temperature between any two replicas is varied every 1000 time units during the course of the simulation.

Depending on cell types, the number of macromolecules present in the cytoplasm varies. We first investigate how the number of particles within a system affects protein folding equilibria. For all simulated proteins, there is a shift in equilibrium (Fig. 3) corresponding to a stabilized native state with respect to increased crowder concentration. Using temperature as a reaction coordinate, the larger proteins (SH3 domain and hemoglobin alpha chain) exhibited sharp transitions between the folded and unfolded states. The peak that separates the two states represents the folding transition temperature, and is equivalent to the midpoint temperature required to unfold a protein in thermal denaturing experiments. The folding transition temperature can therefore be considered a method for studying folding stability of a given protein. Increasing the crowder concentration not only increases the folding transition temperature, but additionally decreases the slope of the folding transition. Even the unfolded states at high crowder concentrations exhibit an overall lower radius of gyration compared to those with lower crowder concentrations.

Insertion of crowders shifts protein folding equilibria toward the native state. The native state (as measured using its radius of gyration) becomes increasingly resistant to thermal denaturation as a function of increasing crowder concentration. (A) **...**

Changes in the slope of the folding transition justify a closer look at the protein folding equilibrium itself. We compute the specific heat as a function of temperature and concentration of crowders (Fig. 4). Not only does the folding temperature peak increase as a function of increasing crowders, but the specific heat required for the folding transition decreases as well. Most interestingly, we find that the folding transition between the two states becomes less cooperative at higher crowder concentrations.

Macromolecular crowding influences specific heat of folding. As the concentration of crowders increase, the specific heat of folding decreases. The transition between the folded and unfolded state becomes less cooperative, denoted by peak broadening. **...**

We perform all-atom simulations of Trp-cage to test whether our coarse-grained models can accurately predict crowding effects despite the loss in atomic accuracy. The range in folding transition temperatures is narrower in the all-atom simulations since at 0% fractional volume there is a peak at 0.63 ε/*k*_{B} compared to 0.51 ε/*k*_{B} in the two-bead Gō model (Fig. 5A). However, the folding temperatures at higher fractional volumes are much similar with values of 0.7 for both models. These differences are attributed to the parameters describing the nature of interactions themselves. Nonetheless the qualitative trend is similar. We quantify the degree of cooperativity by comparing the ratio κ_{2} (defined as the quotient between the van't Hoff enthalpy and the calorimetric enthalpy, see Methods section for further details) across differing crowder concentrations. The decrease in κ_{2} over increasing crowder concentration is consistent with the conclusions derived from the specific heat plots (Fig. 5B). Regardless of resolution, both models of macromolecular crowding display loss of cooperative folding as the concentration of crowder increases.

The composition of macromolecules within a cell consists of proteins with various sizes. Next we evaluate how the size of crowders can affect protein folding equilibria at constant concentration. From a theoretical perspective, larger crowders will occupy more space and exclude more volume. Thus the effect of increasing crowder size should be similar to increasing the concentration of crowders. We compare the average radii of gyration as either a function of crowder diameter or number density (Fig. 6). Our simulation results agree with theory, where the average radius of gyration decreases with respect to increasing crowder diameter or crowder concentration.

Modulation of crowder size or number density affects protein compactness. At a constant concentration of crowders, as the diameter of the crowder increases the fractional volume increases thus leaving limited volume available to the polypeptide. The average **...**

For the smaller proteins Trp-cage and villin headpiece, there is only a slight decrease in R_{G} upon insertion of crowders. For the larger polypeptides, we observe an asymptotic decay of polypeptide compaction when occupying lower fractional volumes. The results suggest that the behavior of protein compaction as a function of available fractional volume is sigmoidal. In the absence of neighboring crowders, the native state of a protein adopts and fluctuates around a certain value of R_{G}. As the occupational volume decreases for the protein, there is a shift (akin to the folding transition temperatures) in equilibria to adopt a compressed state. Within this compressed state the protein is less likely to explore larger conformations, as indicated by the decreasing fluctuations. Larger proteins may possess a critical crowder concentration that forces the protein to shifts its equilibria toward a folded state. The lack of sigmoidal behavior for Trp-cage and villin headpiece may be due to the low initial values of R_{G} or simply sampling only of its asymptotic region. Scaled particle theory predicts there is a critical crowder size associated with stabilization.^{37} However, whether there is a critical crowder concentration remains unknown but our simulations indicate that the effect of increasing crowder size is analogous to increasing crowder concentration although to a lesser degree.

To further characterize macromolecular crowding as a function of size, we consider a scenario where we vary the diameter of the crowder while maintaining a constant fractional volume of crowders by adjusting the concentration of crowders accordingly. From our simulations, we find that smaller crowders stabilize the native state (Fig. 7). While at constant concentration of crowders larger crowders are more effective at excluding volume by virtue of occupying more space, equating the fractional volume between two different types of crowders suggests the opposite effect. The physical explanation is that the void volume available between crowders diminishes with decreasing diameter.^{37} Thus smaller crowders tend to pack more effectively and exclude more volume near the vicinity of a protein, creating a stronger crowding effect.

Effective packing of smaller crowders. We simulated the proteins (A) Trp-cage, (B) villin subdomain, (C) SH3 domain, and (D) hemoglobin alpha chain, with crowders at a constant fractional volume of 26% and varying diameter. By measuring the radius of **...**

In our crowding simulations above, we noted that the cooperativity of folding decreases at higher crowding concentrations (Fig. 4). To examine whether this effect is a general consequence of macromolecular crowding or an effect specific to high concentrations of neighboring solute, we compute the specific heats as a function of temperature for our crowding simulations (Fig. 8A). The strong peaks indicate that the SH3 domain exhibits two-state folding behavior. As the crowder size diminishes, the specific heat required for overcoming the barrier between the two states decreases. Additionally, the temperature peak shifts toward higher temperatures with smaller crowder sizes. We compare the trends in degree of cooperativity as a function of crowder concentration and the degree of cooperativity as a function of crowder size at constant concentration (Fig. 8B). Both conditions exhibit a loss of folding cooperativity which further indicates that this effect is a general feature of crowding.

Loss of folding cooperativity as a function of crowder size. (A) At constant fractional volume, both the specific heat and folding cooperativity of SH3 domain decrease with respect to decreasing crowder diameter. Legend: () 5 Å, ()10Å, **...**

An important question to address is whether or not crowders are merely compacting a protein into the native state more frequently or additionally compressing the protein into other compact non-native states. We examined the distribution of states for the SH3 domain at different crowder sizes (Fig. 8C). In absence of crowders, the SH3 domain exhibits two-state behaviour in our set of simulations. Addition of crowders diminishes the higher energy unfolded peak. As the crowder size decreases, the unfolded state becomes less populated and somewhat more compacted as indicated by a slight shift toward lower energy states. The folded peak becomes more populated and broader at small crowder sizes which suggests that crowders can compress the protein into higher energy non-native states.

The simulations presented here are unique in that they are the first to utilize the replica-exchange method in conjunction with discrete molecular dynamics to study the effects of macromolecular crowding. The combination of these two algorithms has allowed us to explore a large sample of configurations enabling us to conduct thermodynamic analyses.

Modeling proteins with a Gō potential is an effective way to explore conformational space for fast-folding two-state proteins^{9}^{,}^{38} and has been utilized extensively to study protein dynamics including under macromolecular crowding conditions.^{32}^{,}^{39} With the explicit modeling of crowders as hard spheres and proteins at a residue-specific level using a two-bead Gō model, our simulation results are in line with those previously published.^{4} In accordance with theory, increasing the fractional volume occupied by crowders within a system (either by concentration or size) compresses a polypeptide chain into a compact state thereby stabilizing the native state of a protein. Both our results and those previously reported by others suggest that macromolecular crowding can vary a protein's ensemble of states in a non-monotonic fashion. Using an off-lattice model, Cheung *et al*. found that the refolding rates of the WW domain increased nonmonotonically with increasing fractional volumes occupied by the crowders.^{40} By altering either the concentration or the width of crowders, we observe a nonmonotonic shift in protein equilibria towards the native state (Fig. 6).

Our simulations are also able to recapitulate the biophysical phenomena where smaller crowders reduces the void volume near the proximity of a protein and consequently stabilizes a compact conformation.^{37} We find that there is not only a shift in equilibria to favor the native state, but also compacted unfolded states under more aggressive crowding conditions. The shifts seen primarily in the unfolded region are similar to those reported under conditions of macromolecular confinenment.^{10} Overall, these crowding effects become more prominent for polypeptide chains of higher sequence length.

The most surprising finding within our set of simulations is that macromolecular crowding induces a loss of cooperativity in folding. This effect is observed irrespective of how macromolecular crowding is achieved. The implications of specific heat peak broadening are most easily understood in terms of two-state folding kinetics. We have simulated several proteins, in which most are known to exhibit two-state folding dynamics. Trp-cage has been characterized as a two-state folding protein both experimentally^{41} and computationally.^{42} The SH3 domain has also been shown to follow two-state folding kinetics in both computation^{43} and through experiments,^{44} consistent with our current set of simulations that SH3 domain (Fig. 8C). The villin headpiece folding pathway has also been demonstrated to be second-order.^{45}^{,}^{46} It is unclear whether or not the α-domain of hemoglobin follows a two-state folding model though the single sharp transition in the specific heat of folding suggests it does. Nonetheless, we have limited our cooperativity analyses to SH3 domain and Trp-cage.

The Gō model itself does not include actual physicochemical interactions. Rather, the Gō potential is based on the probability of residues being in close proximity within a protein's structure. Hence loss of folding cooperativity may potentially be an artifact due to the coarse-graining of our protein models or the use of a Gō potential. To test the robustness of our findings, we conducted all-atom simulations of Trp-cage using a force field that includes van der Waals packing, solvation, and environment-dependent hydrogen bond interactions. We have chosen Trp-cage since the protein is small enough to ensure adequate sampling of its conformational space despite the increase in atomic resolution. Our results with the all-atom model are consistent with those obtained using the two-bead Gō model (Fig. 4 and and5).5). As the concentration of crowders increase, the two-state folding cooperativity decreases as represented by peak broadening.

Loss of folding cooperativity presents an interesting scenario. Macromolecular crowding effects have been mostly considered a beneficial stabilizing effect for protein folding when regarded as a purely volume exclusion effect. However, the loss of cooperative folding reveals a detrimental consequence of macromolecular crowding. By reducing the two-state folding properties of a protein, the probability of folding into a metastable/alternative state increases. There are biological implications for such behavior since the existence of non-native structures may promote the formation of potentially toxic protein aggregates. High crowder concentrations may also potentially trap proteins into structural intermediates that are present within the folding pathway.

Loss of folding cooperativity leads to alternative structures since crowders may force the protein to adopt a compact state regardless whether it is the native conformation or not. The idea that crowders may induce alternative states has been a topic of recent discussion.^{47} Formation of an alternative structure due to macromolecular crowding was recently discovered experimentally for the VlsE protein found in *Borrelia burgdorferi*. The authors found that if the protein is denatured then refolded under the presence of high crowder concentrations, a change in the overall topology and secondary structure is observed. The conformational change found in this particular protein is biologically significant since crowding may be a method in which antibodies are generated against the bacteria by exposing key residues in its alternative state. Thus macromolecular crowding may modulate changes in conformations by trapping proteins in metastable states.

The Gō model has been previously employed in several studies to examine the cooperativity of folding for numerous proteins.^{32}^{,}^{48}^{,}^{49} Although it is successful in describing the cooperativity of two-state folders as well as identifying key intermediates in higher-order protein folding, the Gō model is biased toward the native state.^{38}^{,}^{50} It is possible to find novel alternative structures, but the enthalpy of those states will be higher by default since the native structure in the Gō model is assigned the lowest energy state. Additionally, the native state bias imposes restrictions on the sampling of non-native structures and thus these alternate structures will occur at a lower probability. More sophisticated force fields are necessary in order to thoroughly characterize alternative structures akin to the previous VlsE example.

Our simulations demonstrate the utility of DMD for investigating questions in macromolecular crowding. Results obtained from weighted histogram analysis are in agreement with theoretical expectations. The loss of cooperativity in two-state folding equilibria presents an exciting new facet of macromolecular crowding. Future applications of DMD aim to explore the effect of crowding on folding equilibria for multi-domain proteins.

We would like to thank Dr Allen P. Minton for helpful discussions. This work was supported by the National Institutes of Health grant R01GM080742.

1. Hall D, Minton AP. Biochim. Biophys. Acta, Proteins Proteomics. 2003;1649:127–139. [PubMed]

2. Minton AP. Biophys. J. 2005;88:971–985. DOI: 10.1529/biophysj.104.050351. [PubMed]

3. van den Berg B, Wain R, Dobson CM, Ellis RJ. EMBO J. 2000;19:3870–3875. DOI: 10.1093/emboj/19.15.3870. [PubMed]

4. Hu Z, Jiang J, Rajagopalan R. Biophys. J. 2007;93:1464–1473. DOI: 10.1529/biophysj.107.104646. [PubMed]

5. Kim JS, Yethiraj A. Biophys. J. 2009;96:1333–1340. DOI: 10.1016/j.bpj.2008.11.030. [PubMed]

6. Minton AP. J. Biol. Chem. 2001;276:10577. [PubMed]

7. Brinker A, Pfeifer G, Kerner MJ, Naylor DJ, Hartl FU, Hayer-Hartl M. Cell. 2001;107:223. [PubMed]

8. Cifra P. J. Chem. Phys. 2009;131:224903. [PubMed]

9. Wang J, Gao H. J. Chem. Phys. 2005;123:084906. [PubMed]

10. Mittal J, Best RB. Proc. Natl. Acad. Sci. U. S. A. 2008;105:20233. [PubMed]

11. Sorin EJ, Pande VS. J. Am. Chem. Soc. 2006;128:6316. [PubMed]

12. Lucent D, Vishal V, Pande VS. Proc. Natl. Acad. Sci. U. S. A. 2007;104:10430. [PubMed]

13. Zhou HX. J. Chem. Phys. 2007;127:245101. [PubMed]

14. Kinjo AR, Takada S. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2002;66:051902. [PubMed]

15. Minh DD, Chang CE, Trylska J, Tozzini V, McCammon JA. J. Am. Chem. Soc. 2006;128:6006–6007. DOI: 10.1021/ja060483s. [PMC free article] [PubMed]

16. Wieczorek G, Zielenkiewicz P. Biophys. J. 2008;95:5030–5036. DOI: 10.1529/biophysj.108.136291. [PubMed]

17. Qin S, Zhou HX. Biophys. J. 2009;97:12–19. DOI: 10.1016/j.bpj.2009.03.066. [PubMed]

18. Ding F, Dokholyan NV. Trends Biotechnol. 2005;23:450–455. DOI: 10.1016/j.tibtech.2005.07.001. [PubMed]

19. Ding F, Tsao D, Nie H, Dokholyan NV. Structure. 2008;16:1010–1018. [PMC free article] [PubMed]

20. Sharma S, Ding F, Dokholyan NV. Front. Biosci. 2008:4795–4808. [PMC free article] [PubMed]

21. Sharma S, Ding F, Dokholyan NV. Biophys. J. 2007;92:1457–1470. DOI: 10.1529/biophysj.106.094805. [PubMed]

22. Ding F, Sharma S, Chalasani P, Demidov VV, Broude NE, Dokholyan NV. RNA. 2008;14:1164–1173. DOI: 10.1261/rna.894608. [PubMed]

23. Davis CH, Nie H, Dokholyan NV. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2007;75:051921. [PubMed]

24. Emperador A, Carrillo O, Rueda M, Orozco M. Biophys. J. 2008;95:2127–2138. DOI: 10.1529/biophysj.107.119115. [PubMed]

25. Alder BJ, Wainwright TE. J. Chem. Phys. 1959;31:459–466.

26. Dokholyan NV, Buldyrev SV, Stanley HE, Shakhnovich EI. Fold. Des. 1998;3:577–587. [PubMed]

27. Sugita Y, Okamoto Y. Chem. Phys. Lett. 1999;314:141.

28. Feig M, Karanicolas J, Brooks CL., 3rd J. Mol. Graphics Modell. 2004;22:377–395. DOI: 10.1016/j.jmgm.2003.12.005. [PubMed]

29. Kumar S, Bouzida D, Swendsen RH, Kollman PA, Rosenberg JM. J. Comput. Chem. 1992;13:1011.

30. Gallicchio E, Andrec M, Felts AK, Levy RM. J. Phys. Chem. B. 2005;109:6722–6731. DOI: 10.1021/jp045294f. [PubMed]

31. Lobanov MY, Bogatyreva NS, Galzitskaya OV. Mol. Biol. 2008;42:623.

32. Knott M, Kaya H, Chan HS. Polymer. 2004;45:623.

33. Tiktopulo EI, Bychkova VE, Ricka J, Ptitsyn OB. Macromolecules. 1994;27:2879.

34. Knott M, Chan HS. Proteins: Struct., Funct., Bioinf. 2006;65:373. [PubMed]

35. Chodera JD, Swope WC, Pitera JW, Seok C, Dill KA. J. Chem. Theory Comput. 2007;3:26. [PubMed]

36. Serohijos AW, Tsygankov D, Liu S, Elston TC, Dokholyan NV. Phys. Chem. Chem. Phys. 2009;11:4840–4850. DOI: 10.1039/b902028d. [PMC free article] [PubMed]

37. Davis-Searles PR, Saunders AJ, Erie DA, Winzor DJ, Pielak GJ. Annu. Rev. Biophys. Biomol. Struct. 2001;30:271. [PubMed]

38. Paci E, Vendruscolo M, Karplus M. Biophys. J. 2002;83:3032–3038. DOI: 10.1016/S0006-3495(02)75308-3. [PubMed]

39. Zhang Z, Chan HS. Biophys. J. 2009;96:L25. [PubMed]

40. Cheung MS, Klimov D, Thirumalai D. Proc. Natl. Acad. Sci. U. S. A. 2005;102:4753. [PubMed]

41. Neidigh JW, Fesinmeyer RM, Andersen NH. Nat. Struct. Biol. 2002;9:425. [PubMed]

42. Snow CD, Zagrovic B, Pande VS. J. Am. Chem. Soc. 2002;124:14548. [PubMed]

43. Ding F, Dokholyan NV, Buldyrev SV, Stanley HE, Shakhnovich EI. Biophys. J. 2002;83:3525–3532. [PubMed]

44. Grantcharova VP, Baker D. Biochemistry. 1997;36:15685–15692. [PubMed]

45. Kubelka J, Eaton WA, Hofrichter J. J. Mol. Biol. 2003;329:625–630. [PubMed]

46. Gong H, Isom DG, Srinivasan R, Rose GD. J. Mol. Biol. 2003;327:1149–1154. [PubMed]

47. Homouz D, Perham M, Samiotakis A, Cheung MS, Wittung-Stafshede P. Proc. Natl. Acad. Sci. U. S. A. 2008;105:11754–11759. DOI: 10.1073/pnas.0803672105. [PubMed]

48. Chan HS, Bromberg S, Dill KA. Philos. Trans. R. Soc. London, Ser. B. 1995;348:61. [PubMed]

49. Dokholyan NV, Buldyrev SV, Stanley HE, Shakhnovich EI. J. Mol. Biol. 2000;296:1183–1188. [PubMed]

50. Guerois R, Serrano L. Curr. Opin. Struct. Biol. 2001;11:101. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |