Dynamic cluster formation determines viscosity and diffusion in dense protein solutions.

We develop a detailed description of protein translational and rotational diffusion in concentrated solution on the basis of all-atom molecular dynamics simulations in explicit solvent. Our systems contain up to 540 fully flexible proteins with 3.6 million atoms. In concentrated protein solutions (100 mg/mL and higher), the proteins ubiquitin and lysozyme, as well as the protein domains third IgG-binding domain of protein G and villin headpiece, diffuse not as isolated particles, but as members of transient clusters between which they constantly exchange. A dynamic cluster model nearly quantitatively explains the increase in viscosity and the decrease in protein diffusivity with protein volume fraction, which both exceed the predictions from widely used colloid models. The Stokes-Einstein relations for translational and rotational diffusion remain valid, but the effective hydrodynamic radius grows linearly with protein volume fraction. This increase follows the observed increase in cluster size and explains the more dramatic slowdown of protein rotation compared with translation. Baxter's sticky-sphere model of colloidal suspensions captures the concentration dependence of cluster size, viscosity, and rotational and translational diffusion. The consistency between simulations and experiments for a diverse set of soluble globular proteins indicates that the cluster model applies broadly to concentrated protein solutions, with equilibrium dissociation constants for nonspecific protein-protein binding in the Kd ≈ 10-mM regime.

The interior of cells is a densely crowded medium, in which macromolecular concentrations range from 90 mg/mL in red blood cells to 300 mg/mL in the mitochondrial matrix (1, 2). Macromolecular crowding influences the stability of proteins, reaction rates, the catalytic activity of enzymes, protein–protein association, and diffusion (3–13). Excluded volume through steric repulsion (14) and attractive protein–protein interactions as well as hydrodynamic interactions affect protein diffusion (6, 15–19). To address the influence of specific protein–protein interactions on protein diffusion (20), crowded solutions with proteins serving as both agents and readout have been studied (5, 14, 16, 21–30).

Experimental techniques to study the effects of macromolecular crowding on diffusion (15) include tracer boundary spreading (14), light scattering spectroscopy (31), fluorescence recovery after photobleaching (FRAP) (32–34), electron spin resonance (35), single-particle tracking (36), fluorescence correlation spectroscopy (FCS) (37–39), quasielastic neutron backscattering (27, 40), and NMR spectroscopy (24, 41, 42). Particle-based simulations complement these experiments (15), treating the proteins as spheres or ellipsoids (20, 43, 44), with residue-level coarse graining (45–47), or as rigid all-atom models (16, 48). Hydrodynamic interactions contribute significantly to the slowdown of protein diffusion in crowded environments (19). In implicit solvent, they are ignored or approximated via the diffusion tensor (16, 19, 44, 49).

Rapid advances in computing hardware and simulation algorithms have opened up the opportunity to study macromolecular crowding using atomistic molecular dynamics (MD) simulations. Explicit solvent accounts directly for excluded volume effects and hydrodynamic interactions and mediates short-range attractive and long-range electrostatic protein–protein interactions (5, 28–30, 50–52). Here, we use atomistic MD simulations of dense protein solutions to calculate the viscosity and protein diffusion coefficients as a function of protein concentration (Fig. 1). Ubiquitin (UBQ), the third IgG-binding domain of protein G (GB3), hen egg white lysozyme (LYZ), and villin headpiece (VIL) are used as model proteins.

Fig. 1.

Representative simulation snapshots of dense UBQ solution (200 mg/mL, ) and GB3 solution (100 mg/mL, ). Blue lines indicate the periodically replicated simulation boxes. Proteins are depicted as green surface (UBQ) or green ribbons (GB3); and ions as blue and cyan balls, respectively; and TIP4P-D water as red sticks.

Soluble proteins self-associate in concentrated solution to form transient and dynamic clusters (19, 24, 53–58). Clustering has also been reported for membrane proteins (59). The influence of cluster formation on the protein translational and rotational diffusivity has recently been addressed by atomistic simulation studies (29, 30). Here, we build on these findings and put cluster formation in the framework of the Stokes–Einstein relations connecting viscosity, cluster size, and diffusion. Central questions are (i) whether the Stokes–Einstein relations remain valid in concentrated protein solutions (60), (ii) how transient protein interactions affect the diffusivity and apparent hydrodynamic radii of proteins in concentrated solutions (42), (iii) how viscosity depends on protein concentration (39, 61, 62), and (iv) whether colloid models apply to concentrated protein solutions.

To address these questions, we perform extensive MD simulations, develop a cluster model of concentration-dependent protein diffusion, and compare our results to the theoretical predictions for hard-sphere (HS) colloidal suspensions without and with attractive interactions. We put our findings in the context of a wide range of experimental and simulation studies and obtain a remarkably consistent picture of the diffusive dynamics in concentrated solutions of soluble globular proteins.

Results

Shear Viscosity Increases Strongly with Protein Volume Fraction.

The viscosity of concentrated protein solutions and of TIP4P-D solvent at different ion concentrations was calculated from the pressure tensor fluctuations (). The quadratic function Eq. fits well over the simulation range, (Fig. 2, squares). We obtained similar fit parameters for UBQ, GB3, and VIL solutions, whereas of the LYZ solutions is significantly lower (Table 1). We calculated mPas from the average viscosity of solutions of TIP4P-D (65) water with 120 mM, 157 mM, and 200 mM NaCl. The values of the calculated LYZ viscosities are in remarkably good agreement with experimentally determined viscosities of the same protein (61). For all four proteins, the viscosities calculated from the pressure tensor fluctuations are bracketed by the experimental values for BSA (62), human serum albumin (HSA), and ribonuclease A (RNaseA) (39) solutions (Fig. 2).

Fig. 2.

Viscosity of concentrated protein solutions as a function of protein volume fraction from all-atom MD simulations (symbols), experiments (dotted lines and gray shading), and HS colloid theory (dashed lines). Blue, UBQ; orange, GB3; green, LYZ; red, VIL. The viscosity was calculated from MD using fluctuations in the pressure tensor (, squares; fit of Eq. , solid lines) and the Stokes–Einstein (SE) relations (Eq. ; small systems , circles; large systems , diamonds). Green and purple dotted lines: fits of Eq. to experimental data on concentrated LYZ (61) and BSA (62) solutions. Magenta and olive dotted lines: exponential fit [equation 4 in Zorrilla et al. (39)] to the experimental data on HSA and RNaseA solutions (39). Black dashed line: Einstein expression for HS colloidal suspensions (Eq. ). Brown dashed-dotted line: Ladd expression (63, 64) for HSs (Eq. ).

Table 1.

Viscosity coefficient , clustering propensity , Baxter parameter , and effective dissociation constant for nonspecific protein–protein binding in concentrated solutions of UBQ, GB3, LYZ, and VIL

Protein	b,mPa⋅s	ζ	τ^cl	Kdcl, mM	τ^RDF	KdRDF, mM	τ^D_t	KdD_t, mM	τ^k_off	Kdk_off, mM	τ^η	Kdη, mM
UBQ	55.3	8.2	0.12	19	0.13	21	0.13	21	0.49	∼80	0.038	6.1
GB3	63.1	8.8	0.11	26	0.19	44	0.13	30	0.34	∼80	0.033	7.6
LYZ	32.5	3.6	0.27	26	—	—	—	—	0.80	∼80	0.071	6.9
VIL	56.2	11.6	0.09	28	—	—	0.13	43	0.24	∼80	0.037	12.2

and are from cluster sizes, and are from radial distribution function g(r), and are from (ɸ) at low protein concentration, and are from binding off rates, and and are from viscosity term b.

The viscosity of dense protein solutions exceeds the Einstein prediction for noninteracting HS colloids (66) even after correction for high concentration (63). The pronounced increase in the viscosity with protein concentration beyond the nonlinear HS viscosity model (63) indicates that short-range attractive interactions between the proteins cannot be ignored. For colloids, the second-order term in Eq. for the viscosity is related to the attraction strength, as measured by the osmotic virial coefficient (67–69). In the following, we use in Eq. with values of listed in Table 1 to account for the dependence of the viscosity on the protein volume fraction.

Translational Diffusion Slows Down at High Protein Density.

As shown in Movie S1 for GB3 at 200 mg/mL with proteins, translational and rotational diffusion in concentrated solution is strongly impacted by protein interactions. For each protein in the simulation box, mean-squared displacement (MSD) curves were calculated and fitted to the Einstein relation in (). The MSD curves of the dense protein solutions averaged over starting times and proteins are linear at times exceeding 10 ns. The translational diffusion coefficients obtained by fitting the Einstein relation to the MSD from 10 ns to 30 ns are therefore long-time diffusion coefficients. The MSD curves of the dilute solutions (one protein in the simulation box) are linear at small delays and were fitted to the Einstein relation from 0 ns to 5 ns. We corrected for large finite-size effects using Eq. , where we used from the quadratic fit, Eq. . The values before finite-size correction are listed in .

After finite-size correction, the translational diffusion coefficient of dilute UBQ is consistent with results of NMR spectroscopy (71, 72) in dilute solution. Our calculated values of LYZ are bracketed by measurements in dilute and dense solutions (21, 24, 25, 71–76). The spread in the measured diffusion coefficients of LYZ is possibly due to differences in pH value, ionic strength, and temperature in the different experiments. All simulation values of in the dilute solution are very close to Hydropro (70) predictions (Fig. 3). The translational diffusion coefficients calculated for the large systems (with proteins) at 200 mg/mL concentration are similar to the values for the small systems () at the same concentration, decreasing slightly with increasing box sizes (number of ). For all proteins, our dynamic cluster model Eq. , developed below, accounts nearly quantitatively for the slowdown of translational diffusion with increasing concentration.

Fig. 3.

Protein diffusion and clustering. (A and B) Concentration-dependent translational (A) and rotational (B) protein diffusion. (A) Dependence of the translational diffusion coefficient on protein volume fraction . Solid circles and diamonds show the finite-size corrected diffusion coefficients of the small systems () and large systems (), respectively. Open squares show Hydropro calculations (70), and other open symbols show experimental data. UBQ data: (71), (72). LYZ data: (73), (74), (75), (21), (76), (24), (25). Lines show the predictions from the dynamic cluster model Eq. with fitted . (B) Dependence of the rotational diffusion coefficients on protein volume fraction . Solid circles and diamonds show the diffusion coefficients obtained from the anisotropic diffusion tensor and corrected for finite-size effects () for the small systems () and large systems (), respectively. Stars () and plus signs () show the finite-size corrected diffusion coefficients from integration of for the small and large systems, respectively. Triangles () represent NMR data for dilute UBQ (blue) (77), dilute GB3 (orange) (78), and dilute LYZ (green) (73) solution. Green crosses () show NMR data for LYZ (25). Open squares show results from Hydropro (70) calculations. Lines show the predictions from the dynamic cluster model Eq. with fitted . B, Inset shows a zoom-in at high protein volume fraction. (C and D) Protein-cluster model. (C) Dependence of mean cluster size on protein volume fraction [circles for small systems (), diamonds for large systems ()]. The dashed line shows a linear fit to the data at 0–100 mg/mL with binding propensity listed in Table 1. Purple circles and diamonds are results from a Monte Carlo simulation of Baxter’s sticky HSs with and particles in the box, respectively, using and finite-range attractive interactions up to 1.05 . The dotted curve shows the fit to the MC data up to . Green crosses show experimental data on LYZ cluster formation (79). (C, Inset) Representative simulation snapshot of dense UBQ solution (100 mg/mL). Blue lines indicate the periodically replicated simulation boxes. Colors and colored labels indicate transient UBQ protein clusters and cluster size , respectively (solvent not shown). (D) Dependence of reduced hydrodynamic radii cubed, , on protein volume fraction . The effective is calculated from the Stokes–Einstein relations for translation, Eq. (circles), and for rotation, Eq. (squares). The dashed lines show the prediction from C with values from Table 1.

Crowding Strongly Affects Rotational Diffusion.

Rotational diffusion coefficients and of the dense protein solutions were obtained by fitting quaternion correlations (80) and by integrating the orientational correlation function , respectively, over the time range 0–100 ns (). Fits to elements of the quaternion covariance matrix are shown in . The resulting rotational diffusion coefficients – are shown in . Fits to the orientational correlation function are shown in .

Orientationally averaged diffusion coefficients and after finite-size correction (81) decrease strongly with increasing (Fig. 3). At infinite dilution, the UBQ and GB3 results are bracketed by the rotational diffusion coefficients obtained from Hydropro (70) calculations and from NMR spectroscopy (77, 78). The experimental rotational diffusion coefficient in dilute LYZ solution reported in ref. 73 is slightly lower than the calculated values, whereas the rotational diffusion coefficients reported in ref. 25 agree well with our data at all protein concentrations. The calculated rotational diffusion coefficients of dilute LYZ and VIL are in fair agreement with Hydropro (70) predictions. As for the translational diffusion, the rotational diffusion coefficients calculated for the large systems () at 200 mg/mL concentration are close to the values of the small systems () at the same concentration, being slightly lower for UBQ and slightly higher for GB3. The dynamic cluster model Eq. predicts the rotational diffusion coefficients of UBQ, GB3, VIL, and LYZ accurately over the entire concentration range, except for the LYZ solution at 100 mg/mL concentration, where the effect of the weak clustering (Fig. 3) is somewhat overestimated.

Diffusion in Dense Protein Solutions Follows the Stokes–Einstein Relation.

Given translational and rotational diffusion coefficients, the viscosity can be estimated from the Stokes–Einstein relations (Eq. ). For all small systems (), we observe excellent agreement between the viscosity calculated from the autocorrelation of the pressure tensor fluctuations and from the diffusion coefficients (Fig. 2), indicating that the dense protein solutions show normal (Stokes–Einstein) diffusion for all protein volume fractions considered here. For the large systems (), slight deviations in and (Fig. 3) lead to strongly overestimated viscosity (diamonds in Fig. 2). Therefore, in practical calculations we advise against calculating the viscosity via the Stokes–Einstein relation, Eq. , because the results are quite sensitive to the uncertainties in and .

Hydrodynamic Radius, Cluster Size, and Diffusion Are Related.

We obtained very similar hydrodynamic radii from the Stokes–Einstein relations for translation and rotation, Eqs. and , respectively (Fig. 3). Therefore, after correcting for finite-size effects with actual shear viscosities , both translational diffusion and rotational diffusion follow Stokes–Einstein theory even at high protein volume fractions. Only values for the small systems () are shown, because only for these systems was the computationally expensive calculation of the viscosity from the pressure fluctuation autocorrelation function performed.

If the increase in viscosity were to capture all factors that contribute to the concentration-dependent slowdown of protein diffusivity, then the hydrodynamic radius, calculated from the Stokes–Einstein relations, Eqs. and , should remain constant at all concentrations. Instead, we observe that the effective hydrodynamic radius cubed, , increases with protein volume fraction . This dependence is consistent with protein cluster formation (Fig. 3). Indeed, when calculating the cluster size distribution based on an -carbon distance cutoff criterion, the cluster size distribution shifts to larger clusters at increasing protein volume fraction (). For protein concentrations up to 100 mg/mL, the mean number of proteins in a cluster grows linearly as , with clustering propensity listed in Table 1 (Fig. 3). The highest concentration (200 mg/mL) was not included in the fit, because the close proximity of proteins causes a significant dependence of the calculated mean cluster size on the cutoff criterion (). Given the linear increase of cluster size with protein volume fraction , the effective hydrodynamic radius cubed should likewise increase linearly with , , where we assumed that the hydrodynamic radius cubed is proportional to the cluster volume. Remarkably, the dynamic cluster model accounts nearly quantitatively for the relative increase of (Fig. 3).

As shown in , the clustering propensity is related to an effective dissociation constant via the protein volume , . Given protein volumes , , , and (), we obtain dissociation constants of mM, mM, mM, and mM with Avogadro’s constant (Table 1).

Effective Viscosity Accounts for Hydrodynamic Interactions.

We investigated whether the effective viscosity captures the indirect (hydrodynamic) effects of increased protein concentration and cluster formation, i.e., all of the effects that are not accounted for by an increased effective hydrodynamic radius. For each protein, we recorded the times at which it is free, i.e., not in any cluster (). For these trajectory segments, we calculated MSD curves () and translational diffusion coefficients . could not be reliably determined for the proteins at the highest concentration (200 mg/mL) due to insufficient sampling of unbound proteins. For all four protein species, the product is approximately constant as a function of concentration (). We conclude that the effective viscosity indeed accounts for the hydrodynamic contributions to the diffusivity slowdown.

Displacement Pair Correlation Shows Contribution from Direct and Hydrodynamic Interactions.

We calculated the displacement pair correlation introduced by Ando and Skolnick (19) (). We analyzed the protein pair correlation for pairs at distances 0.6–3 nm. At distances corresponding to cluster formation, we observed highly correlated motion for all protein pairs at all concentrations and time delays (). At increasing pair distance (∼2–3 nm), the pair correlation decreased gradually.

Protein Binding Interfaces.

The interactions between the proteins in clusters were loose but not entirely random in their orientation (). For UBQ, the preferred binding interface coincides remarkably well with the noncovalent dimer interface reported from NMR measurements (82), more or less independent of protein concentration (Fig. 4). It includes residues 8–11, 20, 24, 25, 28, 31–42, 46–49, 54–60, and 71–76. The C-terminal tail (residues 71–76) and an adjacent relatively hydrophobic surface patch show strong involvement, in line with experimental evidence (82). In dense LYZ solutions, we observed that residues Asp48 and Arg73 contribute most to LYZ–LYZ interaction (). In a Brownian dynamics study (83), these residues were found to play crucial roles in the formation of a LYZ–LYZ encounter complex.

Fig. 4.

UBQ protein–protein interaction interfaces. Red, white, and blue indicate strong, intermediate, and weak involvement in UBQ–UBQ contacts. The NMR experimental noncovalent dimer interface is derived from Liu et al. (82). The rightmost view indicates hydrophobic (white), polar (green), acidic (red), and basic (blue) surface residues. C indicates the C terminus.

Colloidal Suspension Model.

Baxter’s attractive (sticky) HSs (68, 69, 84–86) are widely used as a model for suspensions of interacting colloidal particles. Their association constant is related to the dimensionless Baxter parameter as (derived in ), where is the HS volume. Low indicates high protein stickiness. We conducted Monte Carlo (MC) simulations of Baxter sticky HSs with and finite-range attractive interactions up to 1.05 of and particles in a simulation box. The mean cluster size grows linearly at low protein concentrations, deviating for (Fig. 3). Considering the rough approximations of this model, the cluster sizes of MC simulations of sticky HSs are in surprisingly good agreement with those of the atomistic MD simulations. We can also relate and to the second virial coefficient (69) . Integration of the radial distribution function of the UBQ centers at 30 mg/mL and GB3 centers at 39 mg/mL over the first peak gives us values of and , respectively () (no data for LYZ and VIL, as we did not conduct simulations at low concentrations for these proteins). With and and under the assumption of the sticky HS model, we thus obtain mM and mM, in good agreement with the values from the cluster model (Table 1).

For the sticky HS model, Cichocki and Felderhof (68) derived a low-density expansion of the long-time single-particle diffusion coefficient in terms of the volume fraction,We determined the slope (Fig. 5, dashed-dotted line) from a fit to our protein simulation data for UBQ and GB3 at low , which gives us a value of . From , , and , we obtain apparent dissociation constants of mM and mM that are in good agreement with the values obtained directly from the cluster model and by integration of (Table 1). Since the relative concentration-dependent slowdown of VIL diffusivity is similar to that of UBQ and GB3 (Fig. 5), we also calculated mM using the same . This does not apply for LYZ solutions.

Fig. 5.

Concentration-dependent protein diffusion. (A) Dependence of the normalized translational diffusion coefficient on protein volume fraction . Solid circles show MD data from this study. Open circles show data from simulation studies (16, 26, 29, 87). Open triangles denote data from experimental studies (21–27, 75, 76, 88–93). Hb, hemoglobin; Mb, myoglobin; “Mixed,” CI2 in different dense protein solutions; Ova, ovalbumin. The solid curves show the prediction of the dynamic cluster model (Eq. , no adjustable parameters) for UBQ, GB3, LYZ, and VIL solutions, indicated by corresponding colors. The dashed and dotted curves show the slowdown of predicted from colloid theory on noninteracting HSs by van Blaaderen et al. (94) and Tokuyama and Oppenheim (95), respectively (). The dashed-dotted line shows a linear fit to of UBQ and GB3 at . A plot with experimental studies resolved by symbols is shown in . (B) Dependence of the normalized rotational diffusion coefficient on protein volume fraction . Solid circles show MD data from this study. Open circles show data from simulation studies (16, 29). Open triangles denote data from experimental studies (22, 25, 96). “Cell”: Hb or Mb in different cell types. The solid curves show the prediction of the dynamic cluster model (Eq. , no adjustable parameters) for UBQ, GB3, LYZ, and VIL solutions, indicated by corresponding colors. The dashed curve shows the slowdown of predicted from colloid theory for noninteracting HSs (97) (). A plot with experimental studies resolved by symbols is shown in .

Cichocki and Felderhof (68) also evaluated the quadratic term in the viscosity expansion,Equating with the coefficient from the term in the viscosity expansion fitted to MD data (Eq. and Table 1), we obtain values of and (Table 1), which are again close to those obtained from the cluster model, , and (Table 1).

Dissociation Constant from Off Rate of Nonspecific Complexes.

The cumulative distribution functions of the lifetimes of protein pairs () show that most pairs stay together for 1–50 ns, indicating dynamic clustering according to Liu et al.’s (24) terminology. The lifetimes of protein pairs are independent of the protein concentration, supporting the presence of dynamic protein clusters rather than protein aggregation. From the cumulative distribution function, we obtained the same median protein pair lifetime of ns for all proteins and concentrations and defined an off rate for nonspecific complexes as s−1, assuming exponential kinetics. Assuming in addition a Smoluchowski on rate, , and substituting the Stokes–Einstein relation for translational diffusion, we expect an on rate for nonspecific complexes of s−1. The resulting “kinetic” dissociation constant mM agrees well with those obtained from the other methods (Table 1).

Discussion

Relative Slowdown in Diffusion Is Consistent with Experiment and Cluster Model.

The relative slowdown of translational and rotational diffusion of UBQ, GB3, LYZ, and VIL at increasing protein volume fraction is within the range of published experimental and simulation results. Literature data on the protein-concentration dependence of were normalized to . If only protein mass concentrations were given, protein volume fractions were approximated according to the dependence of UBQ volume fraction on UBQ mass concentration.

The dependence of on protein crowding is surprisingly consistent across the different studies (Fig. 5 and ), considering that they cover both experiments and simulations and report results for proteins of different size, shape, and charge, for different model resolution and different experimental conditions (temperature, pH value). In particular, our simulation results for of dense VIL solutions are in line with a recent study by Nawrocki et al. (29). Likewise, the slowdown of of dense LYZ solutions agrees with the majority of experimental data on LYZ diffusion (green triangles in Fig. 5 and open green symbols in ). We conclude that the slowdown of translational diffusion in concentrated protein solutions is a general feature of soluble proteins and that the extent of the slowdown is governed by the clustering propensity of the protein (notwithstanding additional shape effects, which may not be captured by our choice of globular proteins). We also conclude that the colloidal models of noninteracting HSs (94, 95) ( and dashed and dotted lines in Fig. 5) significantly underestimate the slowdown of translational diffusion of most proteins studied.

We collected data on the slowdown of the rotational diffusion from experimental and simulation studies and normalized the diffusion coefficients to . The slowdown of is less consistent in the considered studies (Fig. 5 and ) than the slowdown of . We suggest that the greater variability in the slowdown of is due to the dependence of compared with the dependence of in the Stokes–Einstein relations. Up to intermediate concentration (100 mg/mL), the slowdown of in LYZ solutions, although in line with experiment (25), does not abide by the cluster model prediction, even though the model takes into account the reduced clustering propensity of the LYZ protein, possibly from its high positive charge (+9) (Fig. 3 and in Table 1). Indeed, the slowdown of at low LYZ concentrations appears to be better represented by a noninteracting HS model. In our simulations, LYZ forms very few clusters at intermediate (100 mg/mL) concentration. At higher concentration (200 mg/mL), we find a significant contribution of higher cluster numbers () (). From small-angle X-ray scattering data, Stradner et al. (79) inferred slightly larger LYZ clusters (green crosses in Fig. 3) compared with our simulations. However, the higher pH in their experiments may have reduced the net protein charge and increased cluster formation. Scattering data were subsequently found to be compatible with lysozyme being largely repulsive (98).

We stress that our dynamic cluster model (Eqs. and and solid lines in Figs. 3 and 5) predicts the relative slowdown in diffusion based on protein cluster size and viscosity. In light of this, the prediction of the model represents the observed and calculated slowdown in diffusion strikingly well. Our cluster model relates differences in the relative diffusivities of different proteins to differences in their nonspecific interactions. This observation provides a physical basis for the protein-dependent concentration scaling factor introduced by ref. 21 to establish consistency with a HS diffusion model (95).

Separating the Effects of Clustering and Hydrodynamics.

Using an elegant MD simulation setup, Nawrocki et al. (30) found that direct protein interactions are the dominant contributors to the slowdown of rotational diffusivity, whereas hydrodynamics play only a minor role. Here, we could show that the effective viscosity accounts for the indirect, hydrodynamic effects of dense solutions on the diffusivity slowdown (). In dense UBQ, GB3, and VIL solutions at 200 mg/mL, the viscosity is ∼2.5-fold increased compared to the solvent viscosity (Fig. 2). The translational diffusion decreases by a factor of ∼4, whereas the rotational diffusion decreases by a factor of ∼6 (Fig. 3 ). The direct effect of protein clustering (i.e., the increase in the effective hydrodynamic radius) accounts for an additional factor ∼ on translational diffusion and a factor ∼ on rotational diffusion. We therefore find that direct clustering interactions contribute to ∼40% and 50% of the total slowdown in translational and rotational diffusivity, respectively, of dense UBQ, GB3, and VIL solutions.

The displacement pair correlation function assesses concerted protein motion and is used to distinguish short-range and long-range interactions in simulations of crowded systems (19). At short distances, the proteins form clusters, and their motion is highly correlated (). The pair correlation decreases at larger protein pair distances, because the short-range interactions fade out and only contributions from long-range effects (electrostatics, hydrodynamics) remain. Interestingly, at low protein concentrations (30–50 mg/mL), the protein motions are correlated up to high distances, whereas the correlation decreases quickly with increasing distance at high protein concentration. This suggests that the motion of proteins at low concentration is more effectively coupled by hydrodynamic and electrostatic forces. At high concentrations, the coupled motion of a protein pair at short distances appears to be effectively quenched by interference of competing proteins. This trend is surprising, because proteins may serve as bridging intermediates at high concentrations, which should effectively increase the pair correlation for the distances considered.

Protein Solutions as Colloidal Suspensions with Attractions.

The effective dissociation constants and Baxter parameters were derived according to the low-concentration behavior of the sticky-HS model for attractive colloidal particles. The values for nonspecific protein–protein binding derived from the structure, dynamics, and viscosity of the protein solutions are fairly consistent with each other (Table 1). As listed in Table 1, estimates from , , , , and vary within about one order of magnitude for a given protein. The corresponding values of also agree well with each other. Interestingly, the decreased clustering propensity of LYZ solutions compared with UBQ, GB3, and VIL solutions is reflected only by a higher , but not by an increase in the effective dissociation constant , because the increased value of LYZ is compensated by its larger volume . In turn, this results in a relatively small for LYZ. For our limited set of four proteins, protein stickiness decreases with size, thereby ameliorating the aggregation problem for large proteins pointed out by Ando and Skolnick (19).

The formation of 3D clusters explains the strong increase of the mean cluster size in sticky-HS solutions at large (Fig. 3). Sticky spheres have no orientational preference, which allows the formation of compact clusters. By contrast, protein–protein interactions are directional and only some protein orientations result in favorable interactions (). The orientational preference disfavors the formation of compact clusters, as can be seen in the formation of only protein connections for a cluster of size up to intermediate (100 mg/mL) protein concentration (). In consideration of these differences, should be considered an effective parameter whose value depends on the property that is fitted and on the fit range. Therefore, we emphasize that the sticky HS model is only a rough approximation of a complex protein solution, which nonetheless explains the behavior of the complex protein systems surprisingly well up to intermediate protein concentration.

From light-scattering experiments, Scherer et al. (57) estimated for concentrated monoclonal antibody mAB2 solution (150 mM NaCl concentration), which is of the same order of magnitude as our simulation results (Table 1). They observed a strong dependence of on the ionic strength, ranging from ∼1,000 to ∼0.5 when increasing the ion concentration from 40 mM to 600 mM. It will be interesting to study the effect of varying ionic strengths on by MD simulations. In contrast to mAB2 solutions, the scattering data for antibody mAB1 solutions could not be represented by Baxter’s sticky HS model due to stronger (and possibly long-range) interactions, whereas an alternative model based on two oligomeric species (57, 58) fitted the data well. This suggests that only weakly interacting proteins are well described by the attractive HS model (57).

Limitations.

The generality of our findings on the connection of cluster size, viscosity, and diffusion is subject to several limitations.

All four proteins studied here are small and globular. This justifies the use of averaged rotational diffusion coefficients and facilitates the comparison with colloid theory, but does not consider the effects of different shapes on the concentration dependence of diffusion. Dense solutions of larger proteins and proteins with varying degree of anisotropy will have to be simulated atomistically to address these issues.

Slight finite-size effects on protein clustering and diffusivity were observed. In the smallest simulation of dense protein solutions (UBQ at 200 mg/mL with proteins), occasional formation of a box-spanning cluster was observed, as seen in the contribution of cluster sizes to the cluster size distribution for this system (). This artifact was not seen in any of the other simulations and does not appear to strongly affect the diffusivity and viscosity of the UBQ solution (Figs. 2 and 3). The translational and rotational diffusion coefficients at high concentration (200 mg/mL) vary slightly for higher protein numbers in the simulation box. The absolute deviations in and between small and large systems are small. However, the relative deviations in and are substantial, overestimating the viscosity calculated from the Stokes–Einstein relations (Fig. 2, circles). On the one hand, we attribute these deviations to incomplete sampling of the extremely large simulation systems. We expect that the simulations with large system size () take more time to equilibrate, thus introducing a bias on the diffusion coefficients, which is not accounted for by the current finite-size corrections (Eq. and ). On the other hand, we observe a significant tail to high cluster numbers in the cluster distribution in these systems (). Whereas the mean cluster sizes in the small and large systems are similar (Fig. 3), occasional large clusters in the large systems appear to slightly suppress diffusion in a way not captured by the smaller simulation systems. Assessing and quantifying both effects in detail would merit a separate study.

Based on the clustering data for UBQ and GB3 solutions, we assume that the cluster size increases linearly with protein concentration up to intermediate (100 mg/mL) protein concentration. Although this simple model works surprisingly well in explaining our calculated diffusion data even at higher concentrations (Fig. 3 ), the actual functional dependence of cluster size on protein concentration may be nonlinear and depend on the specifics of the system (protein type, pH, ionic strength, temperature). Much longer simulations would be needed to precisely determine the cluster distribution and mean cluster size at high protein concentration.

The representation of dense protein solutions by Baxter’s sticky HS suspensions is limited to weak, short-range protein–protein interactions. As scattering data show (57), the sticky HS model does not represent experimental data well for strongly interacting particles or particles with significant attractive long-range (electrostatic) interactions. In our simulations and in cellular conditions, these interactions are effectively shielded by the ions in the solution, making this limitation less relevant in vivo. Nevertheless, it would be interesting to test the applicability of the cluster model for these cases.

Conclusions

By performing all-atom molecular dynamics simulations of dense protein solutions, we found an increase in the viscosity of the solutions at higher protein volume fractions, consistent with experimental results (39, 61, 62). This increase is considerably higher than predicted by colloidal models of noninteracting HSs, stressing the importance of measuring or calculating rather than approximating the viscosity at protein volume fractions approaching cellular crowding conditions. We calculated translational and rotational diffusion coefficients and corrected them for finite-size effects using the respective viscosity of the solution. Translational diffusion and rotational diffusion are strongly affected by protein crowding. For LYZ solutions, experimentally measured diffusion coefficients are available also at high concentration and are in excellent agreement with our simulation data (Fig. 3 ). We calculated effective hydrodynamic radii using the Stokes–Einstein relations and found that a similar increase in the effective hydrodynamic radius can be inferred from the slowdown of translational and rotational diffusion caused by the formation of dynamic protein clusters. Indeed, establishing consistency with the Stokes–Einstein relations requires accounting for protein cluster formation (42) as a result of attractive interactions (60). Overall, we conclude that the concentration dependence of protein cluster size, the translational and rotational diffusion coefficient, and viscosity are consistent with each other (exception: LYZ at 100 mg/mL) and—for the proteins studied here—are explained well by Baxter’s sticky HS model of colloidal suspensions.

Representing the diffusion data in reduced form as a function of protein volume fraction showed that the relative slowdown in translational diffusion is consistent with results from previous studies. The relative slowdown in rotational diffusion shows a larger spread, consistent with the notion that rotational diffusivity depends more sensibly on clustering propensity and thus on the specifics of the protein interactions. Dynamic cluster formation has recently been observed also for membrane proteins (59) and shown to slow down rotational diffusion. In light of our analysis in terms of the cluster and colloidal models, we would expect similar affinities mM for other abundant, soluble proteins with similar size and shape. It will be interesting to explore the limits of our cluster model, i.e., if the cluster model remains valid for large, anisotropic, and more sticky proteins. In this context, the effect of liquid–liquid phase separation (99) on protein diffusivity is an exciting question to address.

We find that the proteins favor certain orientations for interactions and our findings on UBQ contact interfaces are consistent with experiments (82). The protein interactions lead to highly correlated motion at short distances and the correlation is sustained up to larger distances at low concentration. At high concentration, despite increased protein cluster formation, the pair correlation (19) at similar distances is decreased.

In the cellular environment, the situation is complicated by molecular heterogeneity, reactions, partitioning in microenvironments by phase separation, interactions with membranes and structural proteins, and other factors (99–102). Nevertheless, the findings here and in earlier work (103), as well as the observation that diffusion in cell lysates is similar to diffusion in crowded protein solutions (22), suggest that both in concentrated solution and in cells, proteins appear to diffuse not as isolated particles, but as members of dynamic clusters between which they constantly exchange. From the consistency of our diffusivity results with experiments in solution and on the basis of our cluster model, we conclude that—in crowded conditions corresponding to the cellular concentration—the strength of nonspecific protein–protein interactions for abundant proteins such as UBQ should correspond to low-millimolar binding.

We can now carry out atomistic MD simulations of crowded simulations at an unprecedented scale (29), here with up to 540 proteins and 3.6 million atoms in the box simulated over microseconds. Atomistic simulations of solutions of protein mixtures, possibly reflecting the distribution of proteins in the cell, no longer seem out of reach (52). Developments in nucleic acid force fields (104–109) will make it attractive to test the above findings on dense nucleic acid solutions and dense protein–nucleic acid mixtures. Ultimately, the macromolecular diversity in the cell will have to be considered (101) to predict passive diffusion in vivo.

SI Appendix.

contains supplementary text; ; ; Movie S1; and SI references.

Movie S1.

Shown is an atomistic MD simulation of 540 GB3 proteins in concentrated solution (200 mg/mL) at simulation time 0–500 ns. The fully flexible proteins are shown in surface representation and differentiated by color. For clarity, water and ions are omitted. Proteins that seem to appear and disappear traverse the periodic boundaries.

Materials and Methods

MD Simulations of Dense Protein Solutions.

We performed all-atom MD simulations of solutions of human UBQ [PDB code 1UBQ (110)], GB3 [PDB code 1P7F (111)], LYZ [PDB code 1E8L (112)], and VIL [PDB code 1VII (113)] at up to five different densities with (UBQ) and (GB3, LYZ, VIL) proteins in the simulation box. In addition, the most concentrated UBQ and GB3 solutions (200 mg/mL) were simulated with large simulation boxes containing from to proteins (). To mimic an infinitely dilute system, MD simulations with a single protein copy were carried out. The simulation procedures are detailed in .

Viscosity Calculation and Approximations.

The low-frequency, low-shear viscosity of dense protein solutions differs from the viscosity of the pure solvent consisting of only water and ions. We determined and from MD simulations by integration of the autocorrelation functions of the pressure tensor fluctuations (114), as detailed in . The dependence of on the protein volume fraction is well captured by a quadratic function,with parameter fitted to the calculated viscosities, the solvent viscosity averaged over the NaCl concentrations used in this study, and the coefficient 2.5 adopted from Einstein’s colloid theory (66). The experiments of Woutersen and De Kruif (69) have shown that the coefficient increases with the strength of attractions between colloidal particles.

We compared the calculated viscosities to predictions from colloid theory. At low solute volume fractions , Einstein (66) predicted a linear dependence of the viscosity of HS suspensions on ,For higher solute volume fractions, this expression was modified by Ladd (63, 64) towith .

In addition, we also estimated the viscosity by assuming that the Stokes–Einstein relations for rotational diffusion, , and for translational diffusion, , are exactly satisfied, such thatwhere the bar indicates that we average over orientational asymmetries.

Translational Diffusion.

Translational diffusion coefficients were obtained for each protein density by fitting the Einstein relation to the MSD, as detailed in . The diffusion coefficients were corrected for finite-size effects using (115)with and the edge length of the cubic simulation box. In the simulations of dilute proteins, where the protein radius is comparable to the box dimension , we used the expanded correction, (115).

We compared the calculated translational diffusion coefficients to the predictions of a dynamic cluster model without any free parameters,where is the mean cluster size, as defined in and calculated directly from the MD structures. This model is based on the assumption that the slowdown in translational diffusion is linked to the increase in the effective hydrodynamic radius and in the viscosity via the Stokes–Einstein relation for translational diffusion. We also compared the reduced translational diffusion coefficients to the approximate form for monodisperse noninteracting HS colloidal suspensions (94, 95) ().

Rotational Diffusion.

Rotational diffusion coefficients were calculated following the procedure by Linke et al. (80, 116), as detailed in . Additionally, an effective rotational diffusion coefficient was obtained from fits to the orientational correlation function (117, 118), as detailed in .

We compared the calculated rotational diffusion coefficients to the predictions of the dynamic cluster model,where we again assumed the Stokes–Einstein relation to apply with a hydrodynamic radius cubed proportional to the mean cluster size . We also compared the reduced rotational diffusion coefficients to predictions from colloidal models of noninteracting HSs (97) ().

Hydrodynamic Radius.

We solved the Stokes–Einstein relations for translational and rotational diffusion to define effective hydrodynamic radii for translation,and for rotation,Gaussian error propagation of , , and was employed to estimate the errors of and .