Many-particle Brownian and Langevin Dynamics Simulations with the Brownmove package
© Geyer; licensee BioMed Central Ltd. 2011
Received: 11 January 2011
Accepted: 13 April 2011
Published: 13 April 2011
Skip to main content
© Geyer; licensee BioMed Central Ltd. 2011
Received: 11 January 2011
Accepted: 13 April 2011
Published: 13 April 2011
Brownian Dynamics (BD) is a coarse-grained implicit-solvent simulation method that is routinely used to investigate binary protein association dynamics, but due to its efficiency in handling large simulation volumes and particle numbers it is well suited to also describe many-protein scenarios as they often occur in biological cells.
Here we introduce our "brownmove" simulation package which was designed to handle many-particle problems with varying particle numbers and allows for a very flexible definition of rigid and flexible protein and polymer models. Both a Brownian and a Langevin dynamics (LD) propagation scheme can be used and hydrodynamic interactions are treated efficiently with our recently introduced TEA-HI ansatz [Geyer, Winter, JCP 130 (2009) 114905]. With simulations of constrained polymers and flexible models of spherical proteins we demonstrate that it is crucial to include hydrodynamics when multi-bead models are used in BD or LD simulations. Only then both the translational and the rotational diffusion coefficients and the timescales of the internal dynamics can be reproduced correctly. In the third example project we show how constant density boundary conditions [Geyer et al, JCP 120 (2004) 4573] can be used to set up a non-equilibrium simulation of diffusional transport across an array of fixed obstacles. Finally, we demonstrate how the agglomeration dynamics of multiple particles with attractive patches can be analysed conveniently with the help of a dynamic interaction network.
Combining BD and LD propagation, fast hydrodynamics, a flexible protein model, and interfaces for "open" simulation settings, our freely available "brownmove" simulation package constitutes a new platform for coarse-grained many-particle simulations of biologically relevant diffusion and transport processes.
Before any reaction can occur in a biological cell between, for example, an enzyme and its substrate, or before two or more proteins can form a functional complex, the respective partners have to find each other in the crowded interior of the cell. For a full understanding of these association and dissociation processes, be it for a general picture or aimed at designing a drug that enhances or suppresses a certain reaction specifically, a lot of details need to be put together into a consistent picture, rate constants need to be determined, and effects of mutations need to be understood. Many of the details can be determined experimentally. Some of these information are microscopic detailed spatial pictures like crystal structures while others come from macroscopically measured data about turnovers or global reaction kinetics. All these parts of the puzzle can be assembled and studied conveniently by combining them into a computer model and performing simulations. In these in silico experiments all parameters of the system can then be varied to investigate their effect on the association rates and pathways. However, to deal with the large volumes and particle numbers required to describe a biological environment and the often slowly proceeding kinetics, the simulation model has to be both detailed and efficient enough. One often used workhorse technique is Brownian Dynamics, which is named after the botanist Robert Brown, who first described the microscopic random motion of pollen grains in water in a letter to his friends. Nearly eighty years later Einstein could explain his observations . Brown could not see the individual water molecules in his microscope, to him it was a continuous solvent in which the visible pollen grains were floating. Einstein realised that it was the thermal motion of the water molecules which made the pollen grains move. He found that one does not have to know their individual trajectories, but that a heat bath and a Stokesian friction term are enough to describe how they push the large pollen particles around. Einstein's idea was then later cast into the Brownian Dynamics simulation technique . With its continuous solvent only the trajectories of the larger particles of interest are evaluated, which allows for large simulation volumes with many particles and long simulation times. This method has been applied successfully to study, for example, bimolecular association reactions [3–6], the dynamics of colloidal suspensions and polymers [7, 8], or, recently, the dynamics inside the crowded cytoplasm of a cell . Einstein's ansatz to replace the solvent molecules by an effective heat bath and a friction term explains the diffusive behaviour of a single large "Brownian" particle, but it neglects all other effects of the solvent on the interactions between multiple particles. Electrostatic interactions, e.g., are shielded by the ions contained in a physiological solvent. These ions have a size which is comparable to the water molecules and thus they are normally not modelled explicitly. Their effect on the interaction between charges is usually included via Debye-Hückel continuum electrostatics. Other solvent mediated effects are the short ranged hydrophobic and hydrophilic interactions between proteins and the so-called hydrodynamic interactions which stem from the displacement of the solvent by the moving proteins, giving rise to many-particle velocity correlations.
A Brownian dynamics simulation package therefore has to deal with many different kinds of interactions, some of which are direct physical interactions while others are a consequence of the continuum solvent ansatz. For some applications like the association of two proteins, details of their spatial forms and their charge distributions are important, whereas for other applications a model of simple spheres may already capture the important features. A general purpose BD package therefore has to be very flexible and should be easily extensible.
While most of the interactions implemented in "brownmove" follow the usual tried strategies and approximations used in Brownian dynamics simulations of biological systems, there is only limited experience with hydrodynamic interactions (HI) in this context. The basic principles can be found in textbooks since decades, but due to the until now rather expensive numerical evaluation, simulations with HI have essentially only been performed as dedicated tests of the analytical theories while in biologically inspired projects one rather tried to get away without them.
The basic ideas of hydrodynamic interactions are simple. When a particle is moved through the solvent it drags a part of the solvent with it and thus creates a flow field that moves in the same direction as the particle. The second ingredient is Faxen's theorem which states that it takes the same force to move a particle through a solvent at rest as to keep the particle at rest in a solvent flowing around the particle with the opposite velocity. For a sphere moving with a constant velocity through a continuous solvent the flow field was calculated already by Stokes (hence the term "Stokes" friction) and can be found in most textbooks that treat fluid mechanics. A corresponding expression exists for the flow field generated by the rotation of a sphere. For translation the fluid moves in the same direction as the particle with a velocity that--to first order--decays proportional to the inverse distance. Consequently, HI are a long-range phenomenon, coupling all particles in the simulation. For accelerating particles there is, however, no closed form for the resulting flow field, because the information about any changes of the particle velocity cannot travel faster than the speed of sound, which leads to retardation effects.
In simulations with an explicit solvent the solvent molecules take care that the displacements are propagated between the particles, but in any implicit solvent method this hydrodynamic coupling has to be added back explicitly. As flow fields are associated with motion one sees that hydrodynamic interactions are not conservative forces derived from a potential energy but that they describe how the velocities of the particles influence each other. The evaluation of the actual velocity coupling is then an iterative procedure. It starts from the unperturbed, zeroth order velocities of the particles that they would have if there were no HI. From theses the resulting flow fields of the particles are evaluated and then, at the locations of the other particles, converted via Faxen's theorem into an effective force acting on them, which modifies their initial velocities. This process is iterated until convergence, which is slow due to the 1/r-dependence. In addition, the higher order corrections couple three, four, and more particles. For practical applications this series has to be truncated. For spheres this can be calculated analytically and results in a so-called diffusion matrix which converts the external forces acting on each of the particles into the resulting hydrodynamically correlated velocities (the mathematical details can be found, e.g., in ref. ). When only the first iteration is considered, this matrix is called the Oseen tensor . It describes the long-range interactions, where "long-range" means that the particle diameters are much smaller than the separation. This is often expressed as that the Oseen tensor describes HI between point particles. Consequently, as points cannot rotate, there is no rotational coupling on the Oseen level. When additionally the back-coupling from the second particle back to the first is included but no three-particle terms, it is called the Rotne-Prager-Yamakawa (RPY) tensor [21, 22]. This approximation is more accurate than the Oseen tensor but it still underestimates the coupling as the particles come closer. The main reasons why HI is often treated on the RPY level is that it gives reasonably accurate results for most practical applications and that for setting up this hydrodynamic tensor only pairs of particles need to be considered, which results in a runtime that scales quadratically with the number of particles. For any further orders three or more particles have to be considered. Also, rotational coupling can be included on the RPY level . Higher order corrections up to 1/r -7 can be found for example in reference . Forms of the RPY tensor for spheres of different sizes can be found in [25, 26].
When two spherical particles come very close the above explained iteration needs to consider impractically many terms. A more efficient approach is then to expand the hydrodynamic coupling between two particles in powers of their separation, leading to the so-called lubrication corrections (see, e.g.,  for an application to a three-body problem).
From the above explanations one sees that HI add back the mechanical coupling between the particles that was lost in the implicit solvent approximation, albeit on a coarse and approximative level. Also, many proteins are not truly spherical and the correct hydrodynamics is different from the interaction of perfect spheres. The diffusional properties of non-spherical particles can be determined from a multi-bead-model [28, 29], but there is no explicit form of the diffusion matrix for the interaction between such non-spherical objects. Here, one has to resort to a number of spheres per particle. Consequently, there has been a lot of (off the record) debate about whether the computationally expensive HI is worth the effort in BD simulations of biological scenarios where already the protein models and the interactions are modelled by crude approximations only. Due to the high numerical costs of the conventional algorithms this questions--whether and how much HI affect the dynamics in many-protein simulations--has not been addressed on a broad range. Also it is not clear yet, how much higher order corrections affect the results. With our recently presented fast TEA-HI algorithm , which approximates the expensive matrix factorization in the evaluation of the HI but uses the same hydrodynamic matrix, it is now possible to compare simulations with and without HI for many different scenarios. However, we were not the first ones who tried to speed up the evaluation of many-particle HI. Most prominent is Fixman's Chebyshev approximation to the expensive matrix factorization . Other approaches are the Accelerated Stokesian dynamics by Banchio and Brady  or the mean-field-hydrodynamics of Heyes . As should be clear from the above descriptions, HI on this RPY level with its multiple approximations will not be really accurate. However, a comparison between simulations with and without HI can show whether HI makes a difference. If the mere inclusion of RPY-HI proves critical for a certain process and one is interested in accurate quantitative results then a more elaborate method has to be used.
There is a number of methods that range in accuracy and effort between a fully atomistic model, which incorporates all details, and the simplified implicit solvent BD. One approach is to numerically solve the Navier-Stokes equation [32, 33]. Other Methods like Dissipative particle dynamics [34, 35], Multi-Particle Collision Dynamics [36, 37], or the so-called Lowe-Anderson thermostat  use virtual particles to represent momentum "units" of the solvent. Yet another approach is the grid-based Lattice-Boltzmann method where a linearized Boltzmann equation is solved .
Another implication which has to be considered when using the simple and efficient RPY hydrodynamics in Brownian dynamics simulations are the short timescales that are required to describe the fast protein dynamics. Then, as already mentioned above, the flow fields are not the stationary Stokes solutions in an incompressible fluid anymore. On these short time and length scales an explicitly time-dependent method should be used (see for example the discussion in ). However, at the current stage it seems more important to map out for which types of scenarios HI does make a difference and for which it can be neglected.
This publication, which we also use to present our "brownmove" simulation package, is organised as follows. After the above introduction, the following Results and Discussions section starts with two examples on how important hydrodynamic interactions are in coarse-grained simulations of flexible proteins. The first example investigates how the stiffness of a bead-spring polymer affects its diffusional properties, while in the second example a flexible model of a compact protein is built from a number of small beads. Both cases show that hydrodynamic interactions have to be included when both rotational and translational diffusion shall be modelled correctly. In the third example non-equilibrium diffusional transport of simplified proteins through an array of fixed obstacles is simulated. The last example then demonstrates that many-particle simulations can be conveniently analysed and understood quantitatively with the help of a dynamic interaction network. The technical details of the implementation, i.e., of the propagation algorithms, the efficient hydrodynamics, the various interactions, and the available boundary conditions are given in the Methods section after the Conclusions.
In this section we present some example scenarios to illustrate the kind of applications for which brownmove was developed. From the vast number of possible settings we chose two examples that highlight the importance of hydrodynamic interactions in coarse-grained simulations of flexible protein models and two many-particle scenarios. One describes non-equilibrium transport and the other deals with the analysis of many-particle agglomeration.
So far, two projects have been published in which the brownmove simulation package was used. These were many-particle simulations of the association of cytochrome c at a charged membrane  and a coarse-grained model of a small peptide . Some other projects dealing with many-particle agglomeration and transport are currently performed. In the first project , the cytochrome c were set up with the dipolar sphere model which consists of a van-der-Waals sphere of 1.66 nm radius and three charges, a central charge of +7.5 e for the total charge of the horse-heart cytochrome c and two small charges of ± 1.75 e placed on opposite sides of the spherical protein to mimic the dipole . The 30 × 30 nm2 patch of membrane was modelled with a planar van-der-Waals surface, the charges of the lipids were implemented via a Guy-Chapman electrostatic potential plus, for some tests, nine to 25 point charges. The rectangular simulation box was only 20 nm high, so that at the highest concentrations some 300 cytochrome c were enough to fill up the simulation volume. The large bulk of the in-vitro experiments was described via our constant density interfacing algorithm, which simulates particle exchange with an infinitely large reservoir of a given density . With the particle insertion interface providing a constant bulk density we determined binding isotherms which showed that even a few localised point charges on an otherwise continuously charged membrane greatly enhance the adsorption.
In the other project  we built a coarse-grained model of a small peptide in which each residue was set up from one to three van-der-Waals spheres and some point charges. The first van-der-Waals sphere was placed around the Cα atom. For most of the residues a second van-der-Waals sphere was enough to "cover up" the remaining part of the sidechain. Only for the largest residues, a third sphere was needed. The point charges were taken from the crystal structure and placed at the positions of the respective atoms (which generally did not coincide with the centers of the van-der-Waals spheres). Then these rigid residues were connected to their neighbours by springs between the Cα atom positions and, as required, by additional springs between residues spaced further apart. The spring constants were optimised manually against an all-atom molecular dynamics simulation. With this hand-parameterized model peptide we could show how to combine the finitely damped Langevin dynamics and the conventional Rotne-Prager hydrodynamics to yield a fast and stable propagation scheme for these small scales where BD is not applicable anymore . Currently, we are working on a parameterization that allows to convert arbitrary protein structures into flexible coarse-grained models for such LD simulations.
A second way to implement a polymer with constrained internal dynamics in "brownmove" is shown in Figure 2B. Now the beads may rotate and the springs, which only connect direct neighbours, are attached off-center. Together with the effective short-range repulsion between the beads this restricts how far the next bead may move to the side. Note that in this implementation the motion of neighbouring (and potentially non-spherical) beads is essentially free up to the point where they touch while in the first variant the bending angle is confined harmonically. In this second implementation also rotational hydrodynamics have to be considered. A brownmove example for such a bead-train polymer is given as additional file 3. However, for results which are directly comparable to the usual setups we used variant A here. Simulations with polymers of N = 3, 5, 10, and 15 beads, respectively, were run both with and without Rotne-Prager hydrodynamics using our TEA-HI algorithm. The radius and the translational diffusion coefficient D bead of the beads was set to 1, which required a simulation timestep of Δt = 5 × 10-4. The springs between the direct neighbours had a spring constant of k N = 1000 units and the stiffness of the angle confining next neighbour springs of length 5 was varied from k 2N = 0.1 to 1000. The springs between the direct neighbours had a length of 2.5 such that neighbouring beads did not overlap. Overlap between all other pairs of beads was prevented by a short range repulsive potential.
From these observations the following picture emerges. Translational diffusion of our polymers was only weakly affected when the more globular structure of a flexible polymer was "unfolded" by increasing the chain stiffness. However, it was crucial to include HI for the correct scaling of D CM with the chain length. This is even more important when proteins are simulated which can fold into a much more compact structure than these bead-spring polymers with only repulsive interactions between their beads. On the other hand, the rotational motion was very sensitive to the "folding state" of the polymer but essentially unaffected by HI. Thus, if one wants to use a bead-spring polymer as a model for a protein which can change its conformation from a folded globular to a denatured unfolded state, then HI has to be included. Otherwise, only the translational or the rotational diffusional behaviour can be modelled correctly, but not both.
In these bead-spring polymers the direct interactions between the beads, i.e., the springs and the van-der-Waals interactions, were very simple and thus the relative costs of including the hydrodynamic coupling were relatively high. On average, the simulations with HI took about three times as long as the corresponding simulations without HI. However, due to the O(N 2) scaling of our TEA-HI algorithm this ratio was roughly constant for all polymer lengths and thus one can state that when a simulation can be afforded without HI, then normally the corresponding scenario with HI can now be done, too.
For comparison, the translational diffusion coefficient of a sphere scales proportional to its inverse radius . Assuming a constant density, we find that D tr decreases with the third root of the volume. We can thus estimate that in our case where the volume is proportional to the number of beads, N, D tr should decrease as N -0.33. Our results were slightly slower because both the Rotne-Prager tensor and our TEA-HI algorithm underestimate the hydrodynamic coupling at close distances [18, 24]. Also, we did not investigate how the radius of the individual beads affects the diffusion coefficients of the assembled proteins (see, e.g., the discussion in ).
When a protein is modelled from individual beads one can adjust the diffusion coefficient D bead of the individual beads until a best match between an experimentally determined diffusion coefficient and the D tr observed in the simulation is achieved. Here we could not compare to experimental values but we ran another set of simulations without HI where D bead was adjusted such that the center-of-mass diffusion reproduced the values obtained with HI. The required single bead diffusion coefficients were 2.4 × 10-4, 5 × 10-4, 9.7 × 10-4, and 1.2 × 10-3 nm2 ps-1 for the particles with N = 4, 13, 39, and 57 beads, respectively. Consequently, in Figure 6 the obtained D tr values with the faster beads coincide with the results with HI.
However, as shown by Elcock et al. , a multi-bead model "automatically" yields a much more reliable description of the hydrodynamic interactions between different particles. Additionally, the Rotne-Prager tensor is a far-field approximation which works reasonably well when two spherical beads are separated by more than a radius inbetween them. With the smaller beads of a multi-bead model HI thus becomes more accurate for smaller separations between the particles.
Another reason to use multi-bead models is that they allow to describe the internal flexibility of a protein, where often conformational changes occur upon binding or during the catalytic activity or, most prominently, when a protein folds. In these scenarios it would be a bad idea to rescale the single bead diffusion coefficients because even though the overall diffusive behaviour might be described better, the internal dynamics would take place on the too fast timescale of the individual beads. Our "spherical" particles do not undergo conformational changes but it is nevertheless interesting to investigate the dynamics of the individual beads. At very short timescales the elastic springs do not affect the small thermal fluctuations of the beads whereas for very long observation intervals the beads are effectively frozen to their positions in the particle. We thus expect that D bead extracted from a simulation decreases with increasing observation interval from the single bead value D bead(Δt) to the center-of-mass D tr.
Comparing the runtimes for the various protein sizes we again find that the effort per timestep scales quadratically with the number of beads, N, both when HI are included and when not. With HI, the simulations took about three times as long as without. From these runtimes we can estimate that up to a million timesteps of a many-particle simulation with a hundred simple rigid particles or individual shapes can be computed in about one hour on one core of a 2 GHz Core2 Duo CPU. With a typical timestep of 10 ps for a many-protein scenario this amounts to a total simulated time of 10 μs per hour.
Here, however, a word of caution is required. From a theoretical point of view it is not correct to use RPY hydrodynamics which are based on stationary flow fields together with the LD propagation algorithm with its acceleration phases. Here an explicitly time dependent ansatz for HI should be used . Thus, using BD with HI seems to be the better choice regarding theoretical consistency--but using BD for fast processes is questionable, too. To add to this uncertainty, it is also not clear how one should interpret hydrodynamic interactions inside a compact protein where there are usually only a few scattered water molecules? Yet another view to this problem would be to observe that for very short Δt HI effectively become irrelevant for both BD and LD and to conclude from that that for practical applications HI should just be used. They are required for the long time dynamics and do not hurt the details.
Probably, one should use a different interpretation for inter- and for intra-molecular HI in simulations as presented above. The inter-molecular HI describe the solvent-mediated velocity coupling via the resulting flow fields and is thus the "true" HI whereas the intra-molecular HI recover a part of the internal viscosity of the protein: when a multi-bead-model is used to coarse-grain a protein then each bead is a rigid representation of a part of the originally continuously elastic protein. Instead of the continuous deformations of the original protein, now these rigid blocks move relative to each other and HI may be a way to recover at least a part of the viscosity of the protein. Consequently, with all these uncertainties in mind, we find that some more research is required in this area of elastic coarse-grained protein models and fast algorithms for time-dependent HI to finally arrive at both a sufficiently accurate mathematical formulation and the correct interpretation.
Diffusive particle transport through an array of fixed obstacles.
no obstacles (control)
single layer, repulsive, A o = 5 nm
double layer, repulsive, A 0 = 5 nm
double layer, repulsive, A 0 = 6 nm
single layer, attractive, A o = 5 nm
double layer, attractive, A 0 = 5 nm
double layer, attractive, A 0 = 6 nm
double layer, attractive, A 0 = 7.3 nm
The results in table 1 show that with purely repulsive interactions the number of particles in the simulation volume, N p, and the diffusion rate R D decreased as expected with each additional layer of obstacles and also with the obstacle radius A 0. It also took the particles slightly longer to travel through the simulation volume when the second layer was added because now the shorter direct path was blocked and the particles had to diffuse around the obstacles. When the obstacle radius was increased beyond A 0 ≈ 6.5 nm the diffusion current was blocked within the scope of our simulation. Then also T D become very large and no particles reached the ρ = 0 side within the given simulation duration of 100 μs. Intrestingly when a short range attractive interaction was added between the proteins and the obstacles, diffusion was much less hindered. Even two layers of obstacles with A 0 = 5 nm did not reduce the diffusion rate R D. Due to the attraction the proteins stayed closer to the obstacles which reduced the effective "bulk" density. Then, more particles were inserted and N p increased. When the proteins are attracted to the obstacles they temporarily slide along their surfaces and are thus effectively funnelled through the gaps between the obstacles. Consequently, it took them less time to pass the obstacle "barrier" and T D decreased considerably. To shut off the particle transport the now attractive obstacles had to be made so large that the proteins did not fit through the pores anymore. This occurred at A o ≥ 7.4 nm. The efficiency of the surface-induced funneling was so high that even at A 0 = 7.3 nm, where the pores between the obstacles were only slightly larger than the proteins, the diffusion rate was about as high as with a single layer of much smaller repulsive obstacles.
In this simplified scenario both the mobile "proteins" and the obstacles were perfect spheres without any surface roughness and without hydrodynamic interactions which would slow down the protein diffusion close to the large obstacles . In a more realistic setup one would therefore expect that an unspecific attraction makes diffusion between fixed obstacles faster while a corrugate surface would introduce sticking and thus break the surface induced funnelling observed here. Then also HI should be included by implementing the obstacles as mobile proteins that are fixed to their positions by harmonic external potentials. However, such a detailed project with realistic shapes and interactions is beyond the scope of this publication but it can be implemented straightforwardly in brownmove based on the this example.
The simulation was performed with 27 particles in a cubic simulation box of 30 nm length with 3D periodic boundary conditions for 20 μs. Every 10 ns the positions and orientations of the particles were saved to disk. The particle definition and the brownmove simulation setup are given as additional files 8 and 9. After the simulation, which took about 35 minutes on one core of a 2 GHz Core2 Duo CPU, the positions of the red vdW spheres were extracted from the trajectory dump and for each timestep an interaction network was constructed with a distance criterion as indicated in Figure 10B and 10C. For this, each red vdW sphere corresponds to a node of the network and a link was added between all nodes that had a center-to-center distance of less than 4 nm. A spatial configuration as shown in Figure 10B where some of the distance checks are indicated by the green circles would then result in a network as shown in panel C. For illustration purpose the nodes of the network are placed at similar positions as the corresponding vdW spheres, but for our actual network analysis the spatial coordinates were not used anymore once the network was set up. From the dynamic network then typical network measures like the total number of links or the size distribution of the clusters were extracted at each output step.
Panel B of Figure 12 gives the average degree <k>, i.e., the average number of links per node, and the maximal degree during the same simulation run, max(k). To a first approximation each link contributes the same amount to the total binding energy of the system. This panel therefore can be compared to a plot of the total energy of the simulation in a conventional spatial analysis. We see that except for the beginning of the simulation and during the "reorganisation phase" T = 15...17 μs both the averaged <k> and the maximal degree only fluctuated around their typical values and especially the formation of the large clusters cannot be identified from this plot. This can be understood because with the many independent fast binding and unbinding events between the particles each particle will on average have a similar number of binding partners and the maximal degree is limited by the size of the neighbouring particles. Thus only a certain number of particles can come close enough to a given particle to be counted as bound. However, from the comparison of panels A and B we find that the large clusters identified in the cluster size distribution were not very compact. Otherwise, <k> and max(k) would peak at the same time points as the cluster size.
A third view onto the simulation is presented in panel C which gives the number of clusters (connected components) of sizes larger than a single particle, #CC(N>1), and the average size of the clusters of at least two particles, <N>(N>1). These two graph measures again indicate that the simulation was highly dynamic. The number of clusters of at least two connected particles fluctuated between one and about five, indicating that quite often during the simulation there was only one large cluster plus a number of unconnected particles. These events coincide with a large average cluster size. In fact, when there is only one large cluster then the average size is the same as its size.
In general we find that the more local measures such as the degree of a node give less information about the overall state of the simulation than the more global ones like the cluster size distribution. Other helpful measures which were not presented here include the clustering coefficient which quantifies how well or how regular the neighbourhood of a given particle is connected, or the distribution of shortest paths which gives a measure about how densely packed the clusters are.
In this publication we gave four simple examples for simulation projects that can be performed "out of the box" with our brownmove simulation package. The main features presented here are the fast hydrodynamics algorithm, the easy and flexible setup of mobile proteins and of fixed structures in the simulation volume, and how a dynamic network can be used to conveniently analyse a many particle simulation quantitatively and to visualise the fast changes of the time dependent spatial properties.
In the first example we investigated how the stiffness of a bead-spring polymer which was controlled by additional springs between next neighbours, affects the translational and rotational diffusion coefficients. This model system can be compared to a denatured protein in unfolded states ranging from a molten globule like structure of the completely flexible polymer up to a rather rod-like structure when the 1-3-springs are made stiffer. In experiments one finds that a more stretched configuration generally diffuses slower. When hydrodynamic interactions (HI) are omitted--as was often done in BD simulations of protein--the long-time translational diffusion coefficient was completely insensitive to the "folding state" of the polymer. Only when HI were included we observed that the more stretched stiffer polymers diffused slower than the more compact flexible versions. Interestingly, the rotational diffusion, characterised by the relaxation time τrot of the end-to-end vector, was only very weakly affected by HI. While the (constrained) bead-spring polymer can be interpreted as an unfolded protein, the flexible particles of our second example represent folded proteins. They were built from different numbers of small beads placed on a hexagonal close packing lattice and connected by springs to their direct neighbours. Again we examined translational and rotational diffusion and compared our results to the predictions for a sphere with an equivalent radius. In the simulations with HI the predicted scaling both of the center-of-mass diffusion coefficient D tr and of the orientational relaxation time τrot with the number of beads was observed, whereas without HI we could only reproduce either D tr or τrot by rescaling the diffusion coefficients of the individual beads, but not both at the same time. Generally, without HI either rotation was too fast or translation too slow. From these two examples one finds that HI must be included in simulations of flexible protein models when translation, rotation, and internal dynamics are investigated. With HI, the translational diffusion coefficient of the individual beads is then the only parameter that needs to be adjusted. The rotational properties and also the relative timing of internal motions then come "for free". This was already briefly discussed recently by Frembgen-Kesner and Elcock . Whereas there still the usual numerically expensive Cholesky factorisation was used to evaluate the hydrodynamic correlations, we used the recently introduced truncated expansion approximation with its much faster O(N 2) scaling. For these two examples with their very simple beads the inclusion of HI only slowed down the simulations by a factor of about three. When the individual sub units are more complex, i.e., when additionally point charges are included or more than a single van-der-Waals sphere are used to model non-spherical blocks, then the relative costs of HI can even decrease to less than ten percent of the total simulation time. This had been the case for our recently published simulation of a small peptide .
The third example demonstrates how constant density boundary conditions  can be used to model non-equilibrium transport scenarios with the brownmove package. Here, a diffusion current of the small mobile particles across an array of fixed spherical obstacles developed. This can be seen as a greatly simplified model of diffusive transport in a cell where various fixed structures like vesicles or the cytosceleton obstruct free diffusion. As expected, the resulting diffusion current decreased with increasing size and number of the obstacles. Interestingly, when an additional attractive interaction was added between the mobile particles and the obstacles, the diffusional transport was much less affected, because now the mobile "proteins" were guided along the surfaces of the obstacles through the narrow holes between them. One can speculate that a certain "stickiness" between the proteins and the structural elements of a cell might actually help with the efficient transport and thus at least partly compensate for blocking the direct path. In our example the mobile "proteins" and the obstacles were extremely simplified. However, it is straightforward to build much more realistic models of both the proteins and the cellular structures in brownmove by using multiple van-der-Waals spheres per bead or by adding any number of point charges to both the proteins and the obstacles. A more realistic fibrillar structure of the obstacles could for example be implemented by many small spheres along straight or curves lines within the simulation volume. Also crowding could be included by adding a fixed number of mobile particles of another species which then are not exchanged at the reservoirs.
In the last example we simulated a scenario where particles with a "sticky" patch formed temporary agglomerates. Usually such simulations with their fast and frequent association and dissociation events are tedious to analyse. Following a recent project , we mapped the spatial positions onto a dynamic interaction network where each particle is represented by a node and a link is added when the attractive patches of two particles are closer to each other than a specified minimal distance. From this dynamic network we then extracted the distribution of cluster sizes, the average and maximum number of links per node, and the number and average size of the clusters of at least two particles. Together with images from only two snapshots these time dependent network measures allowed to obtain a quantitative picture of the dynamics during the simulation.
The examples presented here are rather templates than actual projects. However, they not only emphasise the importance of HI but they also demonstrate that our freely available "brownmove" simulation package is very flexible and allows to easily investigate a number of scenarios for which a specialised software had to be written before. This especially applies to many-particle simulations with different kinds of particles in "open" systems where, e.g., non-equilibrium transport or reactions, or the association of particles from a large bulk phase onto a surface are considered.
As mentioned above, the particles in a "brownmove" simulation are set up hierarchically from various "objects" (see Figure 1). The "brownmove" package is written in C++ and each of these "objects" is implemented as a class. This design allows to easily extend the current protein model by new interaction types or new functionalities. In addition to the moving particles also fixed objects can be defined which interact with the mobile particles. These fixed objects are used to implement the walls of the simulation box or any rigid structures within this volume.
To model a particle, first a "Protein" object is defined. When this Protein object describes a rigid particle like a folded protein or a colloidal particle then it contains a single "Gestalt" object, which in turn contains the "Shape" objects that are responsible for the interactions. To implement a bead-spring polymer or a flexible protein, the "Protein" object contains multiple "Gestalt" objects and the definitions of the springs between these Gestalten. With this local definition of the connecting springs a single template molecule can be set up from which then multiple copies are drawn and inserted into the simulation during runtime. A "Protein" object thus defines an entity which is inserted or removed from a simulation as a whole.
The next level in the model hierarchy are the independently moving "Gestalt" objects. They keep track of their position and orientation and contain the "Shape" objects. The "GeomShape" is always present. It handles the generation of the random forces and then converts the total force accumulated during the current timestep into the corresponding displacement. Here the Langevin and the Brownian Dynamics propagation schemes are implemented. The other Shape objects model the physical interactions with the other particles, that is, during the force evaluation every pair of "Gestalt" objects compares whether they have Shapes of the same type. Then these two Shapes calculate their contribution to the mutual interaction.
Currently, five types of interactions are implemented. These are shielded Coulombic interactions between sets of point charges in the "EstatShape", effective short range van-der-Waals type interactions in the "vdWShape", bonds from harmonic and quartic terms in the "BondShape", external forces via the "ExternalShape", and hydrodynamic interactions via the "HIShape" objects.
Here, ε is the relative dielectric constant of water and ε0 the vacuum dielectric constant. The shielding from ions in the solvent is captured in the inverse Debye length κ= 1/l D. In most cases, the point charges are embedded in the proteins where there are no counter ions to shield the interaction. This is accounted for by B ik = b i + b k , the sum of the effective burial depths b i and b k . The point charges are defined in an internal coordinate system of the "EstatShape", which in turn can be placed at an arbitrary position and orientation into the coordinate system of the enclosing "Gestalt".
For each van-der-Waals sphere a position within the "vdWShape", a radius, and a "colour" index are specified. For each pair of colours different interaction parameters can be specified to model, e.g., short-range hydrophobic attractions or purely repulsive hydrophilic interaction patches. For numerical stability, the diverging Lennard-Jones interaction can be linearized when the spheres overlap more than a specified distance.
Bonds between "Gestalt" objects of the same "Protein" can be hooked up at arbitrary positions on the "BondShapes". This means that bonds are not restricted to the centers of a "Gestalt" but that they can be attached eccentrically. Currently, harmonic and quartic terms are implemented for the bond potentials.
Similar "hooks" for externally specified forces are provided by the "ExternalShape" objects. Currently implemented are a harmonic potential which allows to confine the position of a "Gestalt" to a certain position, a constant vectorial force which can be used to model, e.g., gravitational forces or a constant fluid velocity, and a shear field.
where the velocity follows the force instantaneously and the displacement due to the external forces increases linearly with the timestep Δt.
While in practical applications both for LD and BD there is the usual upper limit for the integration timestep where the numerical accuracy deteriorates, there is also a conceptual lower limit for Δt in the BD approximation. As a rough estimate, Δt should not be smaller than about ten times the velocity relaxation time τrel. For BD simulations with a single particle type one usually finds an integration timestep which is large enough to be conceptually usable and also short enough for a numerically stable propagation, but this may not work any more when proteins of different sizes are considered. Then, a timestep which is stable enough for the faster smaller particles may be unphysically short for the slower larger ones. For more details on this problem see reference .
In brownmove both algorithms are implemented and can even be used within the same simulation. When a mass is defined for a given particle then the LD propagation scheme is used by the "GeomShape" object, otherwise the BD algorithm. Even though the LD equations of motion (8) and (9) look more complicated than the simple BD equations (10), the additional numerical costs are negligible, because most of the terms are constants and the effort for the actual propagation step scales linearly with the number of particles, while the evaluation of the pairwise forces scales quadratically. We therefore advocate to always use the LD propagation scheme for which only the easily controllable numerical accuracy puts a constraint on the timestep.
Hydrodynamic interactions, which describe the coupling of the particle velocities via the displaced solvent, are modelled in "brownmove" via the Rotne-Prager-Yamakawa (RPY) tensor extended to handle rotation and particles of different radii [21–23, 25, 26]. There is no theoretically rigid formulation for hydrodynamic interactions of overlapping spheres of different sizes. However, a working ansatz was proposed by Durchschlag and Zipper . For practical applications the extended RPY HI-tensor can handle a slight overlap of the beads before the results become numerically unstable. Consequently, in brownmove projects in which the hydrodynamic interaction between different physical particles is defined by their outer surface, any HIShape should be accompanied by a VdwShape with a repulsive short ranged potential to prevent the particles to come closer than their actual, hydrodynamically relevant radii.
where ε = <D ik /D ii> is the average of the normalised off-diagonal entries of the diffusion matrix. Apart from the faster evaluation, this approximate form of the HI has the advantage that the sums in equations (12), (13), and (14) can be evaluated simultaneously with very low memory requirements from the temporarily set up two-body submatrices of the diffusion tensor. With this even large many-particle simulations fit into the fast level-3 cache of current CPUs.
In "brownmove" the above TEA-HI algorithm is implemented in the "HiShape" objects in which currently a single hydrodynamic sphere can be defined. How the hydrodynamic coupling, which is based on an instantaneous flow field, can be combined with the LD algorithm is explained in detail in reference . The basic idea is to apply the HI correlations to average forces, which would lead to the same displacements during the timestep in the BD picture, i.e., when the acceleration is ignored.
Coming back to the problem of bead overlap, in the original algorithm of Ermak and McCammon the results degrade with increasing overlap because the RPY tensor does not cover this regime whereas with Fixman's Chebychev approximation additionally the convergence becomes slower due to the diverging range of the eigenvectors when the spheres start to overlap . This then requires more terms for the approximation to converge to the specified accuracy. On the other hand, our TEA-HI approximation resembles a Taylor expansion that is always truncated after the first correction term regardless of the achieved accuracy. Consequently, the runtime is not affected by the particle separation. As detailed in , the truncation errors increase up to 5% in a dimer of touching spheres which is less than the errors involved in the long-range expansion RPY tensor. Consequently, with a diffusion tensor that correctly describes near-field HI, our method would lead to a relative error of up to 5-10% for the pairs of close particles and less for all the further separated pairs of particles in the simulation.
In "brownmove" various boundary condition can be used. The most simple scenario is an infinite simulation volume which would be used, e.g., to verify the long time diffusional behaviour of a flexible protein assembled from multiple sub-units. To confine the particles, "brownmove" allows to specify combinations of simple reflecting walls, one, two, or three dimensional periodic setups, and walls with a "Gestalt". By defining a "Gestalt" for a wall, not only planar van-der-Waals surfaces can be defined but also static structures like membrane proteins built from van-der-Waals spheres and point charges. A "Wall" object can thus also be used to model a rigid network of microtubili or non-mobile vesicles around which the proteins have to find their way. When periodic boundary conditions are specified, "brownmove" uses image particles during the evaluation of the interactions.
A special type of boundary condition is implemented with the use of particle acceptors and injectors, which allow to define constant density reservoirs and implement reactions at a membrane. The idea of a constant density interface is the following . When a large simulation volume is divided by a virtual wall then particles will cross this boundary with an average rate that depends on their density and their diffusion coefficient. Now the volume on the other side of the virtual wall can be omitted when all particles that cross from the cis to the trans side are removed from the simulation and at the same time new particles are inserted randomly close to the interface with a rate and distribution that corresponds to the assumed density on the trans side. When the density on the trans side is higher more particles will be inserted into the simulation than leave until the density inside the simulation equals the density specified at the interface. Such an interface is especially useful for simulations where the adsorption of particles to surfaces is studied . Here, the constant density interface behaves like an infinitely large bulk and the density above the surface will remain constant no matter how many particles bind. The same algorithm can also be used to model finite reservoirs or non-equilibrium conditions that lead to a diffusion current through the simulation volume. By taking out one type of particles and inserting a different type at the same position reactions like charge transfer can be modeled. More details are given in reference .
For maximal flexibility the output of a brownmove simulation is not directly saved to disk, which would often result in unnecessarily large output files. The particle positions are rather piped into a command that is specified in the setup file. In the most simple case this command dumps the particle positions to a file. If, e.g., from a simulation of a bead-spring polymer only the center of mass and the vector from the first to the last bead is required at each output interval, the output command would extract that information on the fly and only save the processed output to disk. The output commands can be simple scripts, full-fledged analysis tools, or even multiple analysis programs chained together via pipes.
The brownmove simulation package which was presented here is freely available for academic use. The latest version can be downloaded together with documentation and some examples at http://service.bioinformatik.uni-saarland.de/brownmove.
The author thanks Jörg Niggemann and Christian Gorba for fruitful discussions on the design of the brownmove package and for help in the initial steps of the C++ implementation. Uwe Winter implemented the first versions of the Langevin propagation and of the fast HI approximation.
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.