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.
Brownian dynamics of rigid, arbitrarily shaped objects
The propagation scheme for a number of rigid bodies, described either at atomic level or with coarse-grained models, can be written as [23, 26]:
where i runs over molecules, Δt is the time step, and is the vector describing the position of the center of diffusion () and orientation () of the i-th molecule, . is a generalized force acting on a given particle having two components: the force acting on the diffusion center () and the torque . (Δt) is a random displacement vector resulting from the Brownian noise, with zero mean and the variance-covariance given with:
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 ).
Brownian dynamics with hydrodynamic interactions
When hydrodynamic interactions (HI) are considered in BD simulations, the diffusion tensor of the entire system (6N×6N, with N being the number of spherical molecules or the number of spherical pseudo-atoms in case of more elaborate coarse-grained molecular models ) is used to propagate particles. This tensor depends on the system's configuration and varies with time. The propagation scheme can be written as [23, 24]:
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.
When diffusing molecules are described in a reduced way with spherical subunits (e. g., as in case of model spherical proteins [35–37] or bead-models of flexible polymers ), the DLVO (Derjaguin - Landau - Verwey - Overbeek) approach (or the Debye-Hückel approximation), commonly applied in the studies of colloid suspensions, can be used . Spherical molecules (of radii R
) with central net charges (Q
) interact via pairwise additive potentials:
is the vacuum permittivity, ϵ is the dielectric constant of the immersing medium, κ is the inverse of the Debye screening length, and r
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].
In BD simulations employing all-atom models, the forces acting between the surface atoms of molecules can be evaluated using standard (6/12) Lennard-Jones potentials:
However, the well depth ϵ
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.
Applications to diffusion in crowded environments
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].
Protein-protein association in crowded environments
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.
Modeling bacterial cytoplasm
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.