Diffusion in crowded biological environments: applications of Brownian dynamics
© D¿ugosz and Trylska; licensee BioMed Central Ltd. 2011
Received: 9 February 2011
Accepted: 2 March 2011
Published: 2 March 2011
Skip to main content
© D¿ugosz and Trylska; licensee BioMed Central Ltd. 2011
Received: 9 February 2011
Accepted: 2 March 2011
Published: 2 March 2011
Biochemical reactions in living systems occur in complex, heterogeneous media with total concentrations of macromolecules in the range of 50 - 400 . Molecular species occupy a significant fraction of the immersing medium, up to 40% of volume. Such complex and volume-occupied environments are generally termed 'crowded' and/or 'confined'. In crowded conditions non-specific interactions between macromolecules may hinder diffusion - a major process determining metabolism, transport, and signaling. Also, the crowded media can alter, both qualitatively and quantitatively, the reactions in vivo in comparison with their in vitro counterparts. This review focuses on recent developments in particle-based Brownian dynamics algorithms, their applications to model diffusive transport in crowded systems, and their abilities to reproduce and predict the behavior of macromolecules under in vivo conditions.
Intracellular organelles are packed with small solutes, macromolecules, membranes and skeletal proteins. Such crowding is characteristic not only of the cell's interior but also for extracellular tissues . The kinetics and thermodynamics of macromolecular processes and biochemical reactions taking place in vivo are known to be affected by complex, volume-occupied environments [2–5].
More and more experimental methods have become available to investigate macromolecules under in vivo conditions. These include nuclear magnetic resonance [6–8] and fluorescence spectroscopies . Techniques such as SPT (single-particle tracking) [10, 11], FRAP (fluorescence recovery after photobleaching) [12, 13], FCS (fluorescence correlation spectroscopy)  have been applied to measure diffusion constants of macromolecules in the cytoplasm and membranes. All experiments show that diffusion of proteins in vivo is significantly reduced compared to dilute conditions. In the cytoplasm of eukaryotic cells, diffusion of both large and small molecules is slowed down three to four times . FRAP measurements of the diffusion coefficient of GFP in the Escherichia coli (E. coli) cytoplasm [15, 16] yield about 10 times smaller values than at infinite dilution in water.
Theories of diffusion for colloidal soft matter are well established. Within their framework, the behavior of colloidal systems (both dilute and concentrated), consisting of mesoscopically large colloidal particles dispersed in low-molecular-weight solvent, can be simulated, reproduced, and more importantly predicted [17–21]. In cell biology, however, an appropriate theory of diffusion phenomena is still lacking. Complications arise due to the cellular heterogeneity. Moreover, there are different cell types whose compositions vary over space and time on the scales associated with diffusive motion and the knowledge of the surroundings of diffusing species is essential for successful applications of theoretical approaches. Unfortunately, even for the best-studied systems, experimental data is too scarce to describe elementary processes involving macromolecules, their timing and spatial distribution.
Simulation methods used to study the in vivo diffusion of macromolecules should take into account their biological localization, sizes and shapes, and positions in three dimensional space. Intermolecular interactions (both specific and non-specific, repulsive and attractive) should be also included . Ideally, the modeling should be performed at atomic resolution. However, due to the μm size of an average cell and a large number of components occupying its volume, reductionist (coarse-grained) approaches are unavoidable. The simulation methods  that can be applied to directly account for the crowded conditions, sizes and shapes of molecules and their interactions, as well as other features of the cellular environment, include different variants of particle-based Brownian dynamics (BD) [23–26]. In this review, we describe recent applications of BD to study diffusion in the crowded, in-vivo-like systems. Examples of different modeling approaches, either coarse-grained or more elaborate, are presented and discussed.
with the D trans as a proportionality constant.
In biological setup, solutes interact with their environments and other molecules which affects their Brownian motion, thus single-particle diffusion tensors are no longer the unique determinants of diffusion. Nevertheless, the relation between the MSD and time t (not necessarily linear) can still be used to measure the apparent translational diffusion of molecules [4, 11].
Brownian dynamics is a stochastic simulation approach with continuum space and time. In BD the molecules exhibit noise as they are propagated according to the Langevin equation . The equation contains random and frictional forces that are intended to represent the interaction between the diffusing molecules and solvent (treated implicitly). High-friction limit is assumed; momentum relaxation is much faster than position relaxation, thus the equations of motions do not contain particles' velocities. Below, we briefly describe two common BD algorithms.
Here, D is the precomputed and constant in the simulation 6×6 diffusion tensor of the i-th molecule at its diffusion center. Random displacements of particles can be computed via Cholesky decomposition of D . Intermolecular hydrodynamic interactions between diffusing bodies are neglected.
A simplified version of the above algorithm , which assumes isotropic properties of molecules and replaces their diffusion tensors with translational and rotational diffusion coefficients, has been widely used. Apart from treating the diffusing molecules as spheres, this simplification also ignores the coupling between translational and rotational diffusion. It is however questionable, whether using pre-averaged diffusion tensors is justified for nonglobular molecules (e. g., see the discussion on the influence of anisotropy on transport and association rates of model molecules  or on the diffusion of an ellipsoid near planar surfaces ).
When the Rotne-Prager-Yamakwa [33, 34] form of the diffusion tensor is used (see below) the gradient term in the above equation (∇D) vanishes. The 6N-dimensional random displacement vector has zero mean and fulfills < (Δt) (Δt) > = 2DΔt. Again, random displacements of particles can be computed via the Cholesky decomposition of D . However, the decomposition of D must be performed at each BD step due to its dependence on the positions of particles.
The influence of environments on the diffusion of macromolecules manifests through intermolecular forces acting on individual particles. The types and functional forms of interactions included in BD simulations vary depending on the level of detail used to describe the studied systems. Different types of intermolecular interactions that can be modeled in BD simulations are briefly described below.
Where ϵ o is the vacuum permittivity, ϵ is the dielectric constant of the immersing medium, κ is the inverse of the Debye screening length, and r ij is the separation between molecules. However, the above description is not valid for highly charged systems such as nucleic acids where the linearized form of the Poisson-Boltzmann equation  is not applicable and there is strong electrostatic coupling between the molecular and ionic species. In such cases, the net charges can be replaced with renormalized charges .
A more sophisticated approach to treat electrostatic interactions in BD simulations of rigid molecules with atomic detail [42, 43] is based on the model of effective charges developed by Gabdoulline and Wade . In this approach, each molecule is described with a limited set of optimized, discrete, effective charges immersed in an uniform dielectric. These effective charges are derived by fitting the electrostatic potential resulting from the Debye-Hückel approximation to the external molecular potential obtained from the numerical solution of the Poisson-Boltzmann equation. Electrostatic energies (and thus forces) are evaluated using the values of the electrostatic potential (derived from the solution of the Poisson-Boltzmann equation on a three-dimensional grid) generated by one molecule at positions of effective charges of the second molecule. The electrostatic potential grids and effective charges are precomputed so when multiple molecules are simulated, the electrostatic potential grids need to be translated and reorientated at each simulation step along with translating and rotating molecules. This approach ignores the dielectric polarization effects that occur when two molecules come close to each other. These effects may be particularly important in the crowded systems, where intermolecular distances are comparable with the sizes of molecules. Polarization can be accounted for by using an analytical formula proposed by Gabdoulline and Wade . The performance of the effective charges approach is worsened for highly charged molecules in cases where the nonlinear form of the Poisson-Boltzmann equation should be used to evaluate electrostatic potentials. The Debye-Hückel approximation is not valid in these cases, especially in the vicinity of the electrostatic field sources .
Non-specific interactions between molecules can be modeled with the Lennard-Jones type functions that are commonly used in molecular dynamics simulations to compute the van der Waals interactions between atoms; at large intermolecular separations molecules attract each other, while at small separations the interactions are repulsive and molecules are impenetrable. The Lennard-Jones functional form and parameters for BD simulations with coarse-grained spherical models [35, 37] can be, for instance, taken from the ones derived for large colloidal spheres [46, 47].
However, the well depth ϵ LJ needs to be treated as a free parameter and adjusted to reproduce experimental data .
Nonpolar (hydrophobic) interactions can be included in atomistic BD simulations by assuming their proportionality to the amount of the solvent accessible surface area (SASA) that gets buried when molecules come close to each other. However, computing of SASA can be slow and to make SASA-based approaches efficient, approximations are needed. Examples of how to incorporate the hydrophobic effects based on SASA in BD simulations can be found in the work of Elcock and McCammon  or Gabdoulline and Wade .
While SASA-based models are commonly used to evaluate nonpolar interactions, it has been shown that including the solvent accessible volume (SAV) terms may be essential to accurately model the nonpolar solvation, especially at atomic-length scales . Accurate treatment of nonpolar energies and forces is crucial for understanding any biomolecular process that involves changes in solvent accessibility. The described inability of SASA-only models  to reproduce and predict solvation energies and forces indicates that there is a need to include a more complete theory of nonpolar solvation in atomistic BD simulations.
It was also proposed that hydrophobic interactions can be accounted for by modifying the van der Waals interactions and assigning more favorable interaction energies to nonpolar surface atoms than to polar ones . With such an approach hydrophobic interactions can be evaluated more rapidly, however its correctness is questioned (see the discussion in ).
In BD simulations water molecules are not treated explicitly but the influence of solvent on the diffusion of suspended molecules can be included by proper treatment of hydrodynamic interactions. However, even though their importance is widely recognized [32, 37, 51], accounting for HI is no simple matter due to their long-range and many-body character  resulting in high computational cost. Therefore, one is often tempted to neglect HI and, consequently, the correlations of motions of molecules.
In the Ermak-McCammon algortihm  (Equation 5), hydrodynamic interactions between spherical particles can be included by using the Oseen , Rotne-Prager  or Rotne-Prager-Yamakawa  forms of diffusion tensors (without HI the diffusion tensor is represented by a diagonal matrix). At each step of a BD simulation, the diffusion tensor of the whole system is factorized using the Cholesky decomposition  to obtain hydrodynamically correlated displacements of particles. The Cholesky decomposition scales as N3 (with N being the number of particles), while evaluating direct intermolecular interactions, that are assumed pairwise, scales as N2. It is possible to lower the cost of evaluating HI using the Chebyshev approximation proposed by Fixman , but the HI-related computational overhead is still considerable (N2.5). Recently, Geyer and Winter proposed a faster algorithm that scales as N2 and is based on a truncated expansion of the hydrodynamic multi-particle correlations as two-body contributions .
An additional computational cost arises when HI are evaluated in finite simulation boxes containing molecules. This cost is due to the fact that for evaluating diffusion tensors Ewald summation has to be used [56, 57]. While the direct intermolecular interactions in such periodic systems can be computed using less expensive cutoff -based schemes, applying these schemes to evaluate HI leads to diffusion tensors that are not positive definite (a condition that has to be met to use diffusion tensors as covariance matrices .)
The treatment of hydrodynamics within the framework of the Oseen, Rotne-Prager or Rotne-Prager-Yamakwa tensors is appropriate for very dilute systems capturing only the far-field, two-body hydrodynamic effects. However, in crowded biological systems, the many-body HI and lubrication forces (resulting from the thin layer of solvent that separates the nearly touching particles) may significantly influence diffusion. One crude way to include the many-body HI in BD simulations is to use mean-field approaches [58–60] that are based on scaling the diffusion coefficients of particles according to local volume fractions, as for example in the work of Sun and Weinstein .
The most advanced approach to treating HI, used recently in BD simulations of a biological system , was developed by Brady and Bossis . Their Stokesian dynamics (SD) method includes both far-field and near-field many-body contributions to Brownian forces acting on particles. Moreover, particular implementations of SD enable one to efficiently simulate long-time dynamics of large multicomponent systems. For example, the accelerated SD method, described in  scales as N1.25log(N).
Unfortunately, modeling of HI in BD simulations is currently limited to systems composed of coarse-grained particles (either spherical molecules or composed of a small number of spherical units - pseudo-atoms). While such models perform well for colloids, a question yet unanswered is whether they are sufficient for crowded, heterogeneous biological systems.
Typically, the setup of BD simulations aimed at investigating diffusion in crowded solutions is the following. A number of molecules (usually 102 - 103), described either at the atomistic level or with coarse-grained models, is contained inside a primary simulation box. Periodic boundaries are applied to reduce the edge effects and the influence of the finite size of the system. BD trajectories of all molecules are generated with the algorithm of choice. The BD simulation time scale should be long enough so the relative displacements in the system be at least comparable with the dimensions of molecules; the simulation time scales range from microseconds to milliseconds, while typical time steps are about a few picoseconds. When one considers the number of simulated molecules, the need to describe them in as detailed way as possible, and the complexity of interactions calculated in a single simulation step, it becomes clear that BD simulations of crowded systems require a lot of computational resources.
Below we present a few recent (except for one) examples of BD simulations of macromolecular diffusion in crowded media. We describe the studies that best illustrate the potential and shortcomings of the BD technique applied to in vivo-like systems and the ones that introduce a significant element of novelty. Some of these computational studies relate to experiments - experimental data were either used for parameterization or for direct comparison with simulations.
One of the first studies that used the BD technique to investigate macromolecular diffusion in congested media was that of Dwyer and Bloomfield  in 1993. They performed BD simulations of probe and self-diffusion of concentrated solutions containing short DNA polymers (modeled as wormlike chains, i.e. strings of spherical beads) and the bovine serum albumin (BSA) protein modeled as a sphere. Even though their DNA model lacked the detailed all-atom description, it reproduced the overall translational and rotational behavior of linear DNA molecules . The authors used an algorithm similar to that of Ermak-McCammon , however, they neglected HI, because the pairwise approaches are inadequate for concentrated solutions. They also used the Debye-Hückel approach to model long-range electrostatics and the repulsive part of the Lennard-Jones potential to ensure volume exclusion. This model, while very simple, turned out to be realistic enough to simulate diffusion in the DNA/BSA system and accurately reproduce experimental diffusion coefficients at 0.1 M ionic strength for the self-diffusion of BSA and the probe diffusion of BSA in DNA (both over a wide range of BSA/DNA concentrations).
Simulations of the diffusion of BSA in the DNA solution conducted at low ionic strength of 0.01 M gave worse agreement with experiments which the authors attributed to the lack of HI and limitations of the applied electrostatic model. The work of Dwyer and Bloomfield is significant because it shows that even very simple modeling approaches can accurately predict the behavior of rather complex systems.
We note that the models similar to the one of Dwyer and Bloomfield are still used, for example, in computational studies of facilitated diffusion of proteins on DNA. The description of facilitated diffusion phenomenon and models of non-specific protein-DNA interactions can be found in references [65, 66].
Most biological reactions require some sort of diffusion-mediated encounter. Here we describe BD applications that directly investigate the influence of crowded environments on diffusional encounter and association of molecules.
Sun and Weinstein  simulated systems of spherical particles, mimicking globular proteins. Particles were modeled as hard spheres; the authors used the elastic collision method to effectively deal with particle overlaps . Electrostatic interactions were introduced using the Debye-Hückel approach. A notable fact is that the authors attempted to model also many-body HI. They used the translational part of the Ermak-McCammon algorithm with the diagonal form of the diffusion tensor. However, they applied an average hydrodynamic correction and scaled the diagonal terms of the diffusion tensor (i.e., the single particle diffusion coefficients) with functions that depend on the instantaneous local volume fraction of each particle [58–61]. The authors investigated how the size of background molecules, volume fraction, electrostatic interactions, and HI influence the diffusion of a tracer molecule. Their study also quantitatively described the dependence of diffusion-limited reaction rates on crowding in model systems. Using a very approximate approach, they showed a significant effect of HI on diffusion and association rates. While this review focuses on simulations of multimolecular systems, the importance of intermolecular HI in modeling the association kinetics in a two-molecular case was also described by Frembgen-Kesner and Elcock . A number of simplifications was used in the described work of Sun and Weinstein . Simulated systems were rather homogeneous in comparison with in vivo environments. Molecules were treated as spheres thus their association was non-specific (the model of particles with uniformly reactive surfaces was used). Anisotropy of particles, which may play a role upon association, and specific intermolecular interactions were neglected. However, their work is important as it showed the necessity to include HI in BD simulations of crowded systems.
One of the conclusions of the work of Sun and Weinstein  was that crowding slows down the non-specific association of molecules. The two more recent BD studies focused on the influence of crowders on specific interactions in protein-protein association and considered the fact that partner proteins need to be properly oriented to associate. Wieczorek and Zielenkiewicz  used BD to simulate the diffusional encounter of HEL and HyHEL-5 proteins in solutions of hard-sphere crowders. Both proteins were modeled with atomistic details, however, only the excluded volume interactions were evaluated; simulation steps leading to overlaps were rejected and repeated with different random displacements, so no systematic forces were explicitly applied. The crowders and proteins were modeled as rigid bodies with isotropic diffusion tensors (both HEL and HyHEL-5 were represented as spheres with equivalent hydrodynamic radii) and propagated with the algorithm described in . Hydrodynamic interactions were neglected. The authors concluded that if the definition of an encounter complex accounts for the orientation of associating molecules, crowding can accelerate the association rate.
Similar observations were made by Kim and Yethiraj . Using spherical protein models with non-uniformly reactive surfaces, they showed that association of model molecules can be either slowed down or accelerated, depending on the applied reaction criteria. These authors recently used BD to investigate the effects of crowding on the association at model membranes . Both studies of Kim and Yethiraj [69, 70] focus on how the excluded volume effects influence the reaction rates.
In general, the above simulation results are in agreement with purely theoretical predictions . However, they lack a close connection to experiments and the ability to reproduce heterogeneous in vivo conditions.
In the recent work of Elcock's  and Skolnick's  groups BD was applied to investigate the diffusion in cytosol-like systems with the Escherichia coli cytoplasm as a model system. The cytoplasm model of McGuffee and Elcock  contains 50 different types of the most abundant proteins and nucleic acids present in the E. coli cytoplasm. Additionally, a few molecules of the Green Fluorescent Protein (GFP) were included, for the sake of comparison with experimental data. All molecules were represented with atomistic details. The algorithm of rigid bodies was applied in BD but all molecules were treated as isotropic, so the average rotational and translational diffusion coefficients were used instead of 6×6 diffusion tensors .
Intermolecular interactions included electrostatics, modeled using the effective charges approach of Gabadoulline and Wade , and combined van der Waals and hydrophobic interactions described with the Lennard-Jones potential , with the well depth treated as an adjustable parameter. Hydrodynamic interactions were neglected. The authors investigated how different energy models affect the translational diffusion coefficient of GFP. When only steric interactions (modeled with the repulsive part of the 6/12 Lennard-Jones potential) were considered, the BD simulated translational diffusion coefficient of GFP was 3 to 6 times higher than the coefficient estimated experimentally. Even including electrostatics did not significantly reduce the discrepancy between the experimental and computed values. Nevertheless, the authors obtained a close agreement with experiments after adding the short-range Lennard-Jones attraction and tuning its well-depth. The observation that purely steric interactions are not able to account for the reduction in proteins' mobilities in vivo concur with previous computational studies .
While adjusting the short-range attractive interactions led to reproducing the in vivo GFP diffusion coefficient, it has been later pointed out  that the magnitude of attractive interactions can be tuned to set the diffusion coefficient to any chosen value. McGuffee and Elcock investigated also the character of translational diffusion in cytoplasm with different interaction energy models. When only the steric interactions were considered, for the three chosen proteins they observed normal diffusion at short times (ps), subdiffusion at moderate time scales (ns) and normal diffusion at long times (μs). When electrostatic and short-range attractions were included, sub-diffusive behavior was noticed even for short observation times.
The approach to model diffusion in the E. coli cytoplasm taken by Ando and Skolnick  is different than the one described above. These authors investigated the effect of macromolecular shapes on the diffusion in crowded, heterogeneous environment, the role of HI, and the interplay between hydrodynamic and non-specific attractive interactions. Their cytoplasm model consisted of 15 different types of molecules, each described with a set of spherical subunits positioned at the Cα protein atoms and the P, C4', N1, and N9 atoms of nucleic acids. The authors dealt with the anisotropy of macromolecular shapes using the rigid body BD algorithm  that describes molecules with precomputed, 6×6 diffusion tensors. Notably, this is the first BD study of a multimolecular system that takes into account the anisotropy of diffusion. To investigate the role of HI they also simulated the system of equivalent spheres corresponding to the cytoplasm model (all molecules in the system were represented as spherical particles and assigned hydrodynamic radii based on translational diffusion coefficients computed from their atomic structures ) with the far-field many-body HI interactions and lubrication forces treated with the method of Banchio and Brady . Soft-sphere harmonic potential was used to model repulsive interactions between particles. Attractive interactions were described by the Lennard-Jones potential function with a correction factor accounting for the roughness of molecular surfaces. Electrostatic interactions were modeled using the DLVO model.
Based on simulations of molecular-shaped and equivalent sphere systems, the authors showed that the effects of macromolecular shapes on long-time translational diffusion are rather small in the studied concentration range of 250 - 350 . This is an important result that justifies the use of isotropic models to study translational diffusion in the crowded systems. The influence of macromolecular shapes on rotational dynamics, which is an interesting issue as well, was not investigated. Their BD simulations with only steric repulsion between particles also resulted in overestimated translational diffusion coefficients. However, when HI were switched on, the computed long-time translational diffusion coefficient of GFP was in good agreement with experimental value; this is another important result because this model did not include any tunable parameters. The authors concluded that the macromolecular motions in vivo are likely to be dominated by steric interactions and hydrodynamics. Additionally, they observed that when non-specific attractive interactions between particles (either molecular-shaped or spherical) are included, the reduction of translational diffusion coefficients depends much more on molecular sizes in comparison with the simulations where HI were enabled.
Both studies described above [37, 43] are state of the art applications of BD in the field of macromolecular diffusion in biological systems, as they both take into account the heterogeneity of the environment and intermolecular interactions beyond simple excluded volume models. Both are also closely connected to experimental data.
This short review focuses on the current status of the particle-based Brownian dynamics technique and its applicability to model diffusive phenomena in biological environments. We described a few BD studies of three-dimensional macromolecular diffusion in aqueous cellular compartments under molecular crowding conditions. However, there are various other BD applications, not presented here, that consider diffusive phenomena in cellular architectures, such as the cytoskeleton  or membranes .
We presented different modeling approaches, from simple excluded-volume models to most sophisticated atomic level descriptions of the in vivo-like systems. The complexity of intra and extracellular environments, their dynamics, and the number of different interactions that influence the diffusive behavior of macromolecules, make computational studies of crowding extremely difficult. While it is desirable to construct models that treat biological systems at atomic resolutions and include sophisticated interactions [42, 43], the recent model of the E. coli cytoplasm  shows that less elaborate, coarse-grained approaches should not be disregarded, at least as far as diffusive transport is concerned.
Obviously, current BD algorithms need development. These include accurate and computationally efficient models of direct, both specific and non-specific, intermolecular interactions as well as HI models that up to date have been largely neglected in the simulations of biological systems.
BD simulations should be (at least ideally) able to reproduce and predict the in vivo dynamics. However, the computational cost of sophisticated BD simulations may be too high, even with the growth of computing power and available modern technologies. Thus, the long-term goal should be to systematically develop efficient algorithms interfacing different modeling approaches - from atomistic to coarse-grained models of macromolecules, mesoscopic models of biological environments, and appropriate boundary conditions for different cell compartments. New approaches and models should be verified experimentally to establish their quality and predictive power.
The authors acknowledge support from University of Warsaw (BST 2010, G31-4), Polish Ministry of Science and Higher Education (N N301 245236, N N519 646640) and Foundation for Polish Science (Focus program and TEAM/2009-3/8 project co-financed by European Regional Development Fund operated within Innovative Economy Operational Programme).
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.