Abstract
Background
Brownian dynamics (BD) simulations can be used to study very large molecular systems, such as models of the intracellular environment, using atomicdetail structures. Such simulations require strategies to contain the computational costs, especially for the computation of interaction forces and energies. A common approach is to compute interaction forces between macromolecules by precomputing their interaction potentials on threedimensional discretized grids. For longrange interactions, such as electrostatics, gridbased methods are subject to finite size errors. We describe here the implementation of a DebyeHückel correction to the gridbased electrostatic potential used in the SDA BD simulation software that was applied to simulate solutions of bovine serum albumin and of hen egg white lysozyme.
Results
We found that the inclusion of the longrange electrostatic correction increased the accuracy of both the proteinprotein interaction profiles and the protein diffusion coefficients at low ionic strength.
Conclusions
An advantage of this method is the low additional computational cost required to treat longrange electrostatic interactions in large biomacromolecular systems. Moreover, the implementation described here for BD simulations of protein solutions can also be applied in implicit solvent molecular dynamics simulations that make use of gridded interaction potentials.
Keywords:
Continuum solvent electrostatics; Ionic strength; DebyeHückel; PoissonBoltzmann equation; Brownian dynamics simulation; Protein diffusion; Discretization grid; Finite difference; Second virial coefficient; Small angle scattering intensityBackground
Simulations of concentrated solutions of macromolecules such as those designed to mimic the intracellular environment, are becoming feasible due to improvements in computational power and simulation methods [15]. Given that even for simulating a small volume of a protein solution, several hundreds of proteins have to be taken into account, coarse grained methods, which neglect atomic details, e.g. by treating each protein as a sphere, are often applied [6].
However, to understand the effects of differences in protein sequence or point mutations from simulations requires a more detailed level of modeling. Explicit inclusion of atomic detail can be computationally demanding and therefore, approximations and calculation strategies are required to make the simulations feasible. A commonly employed approach is to retain atomic detail for the macromolecules while treating them as rigid bodies in continuum solvent. Apart from restricting the number of degrees of freedom considered in the simulations, this treatment permits interaction forces between macromolecules to be computed efficiently by precomputation of their interaction potentials on threedimensional discretized grids. Thus, during the simulations, forces can be computed by considering the interactions of each atom of each macromolecule with the interaction potential grids of the other macromolecules. Grid formalisms for intermolecular interactions are extensively used for macromolecular docking methodologies [7,8], binding site determination [9], as well as in structure determination from electron microscopy maps [10,11]. A major drawback of gridded potentials is, however, the occurrence of finite size problems [3]. To minimize truncation errors in computing energies or forces, the interaction potential must be small at the edges of a grid. For molecular electrostatic potentials, the longrange nature of the Coulombic interaction, especially at low salt concentration or for highly charged macromolecules, means that very large grids are often required. For example, at 5 mM ionic strength, the Debye length of the solution is 43 Å. For a small globular protein with a radius of 20 Å and a net charge of + 10e, the electrostatic grid dimensions should be at least 200 × 200 × 200 Å to obtain an electrostatic potential of ≈ 0.1 kcal/mol/e at the grid edges. Assuming a grid spacing of 1 Å, the grid must have at least 201 × 201 × 201 points. This grid size is not a problem when a single small protein is considered but becomes an issue when simulating a periodic box containing several hundreds or thousands of proteins in solution. The grid size can also be a problem for memory usage in calculations for one or a few large macromolecules.
One solution to this problem is to use multiple focused grids with different grid spacings centered on each macromolecule: a detailed potential grid with a small grid spacing for representing the electrostatic potential at shortrange and a coarse grid with a larger grid spacing for the longrange part [1]. Another solution, which will be described in this paper, is to exploit the fact that beyond a certain distance from the surface of the macromolecule, the electrostatic potential becomes centrosymmetric. Thus, a cubic gridded potential is used for the shortrange part of the electrostatic potential up to a defined distance threshold and a continuous screened Coulomb potential is used beyond this distance. The distance threshold corresponds to the radius of the largest sphere enclosed by the grid.
We have recently developed a Brownian dynamics (BD) method for simulating many macromolecules (10^{2}10^{3}) described as atomically detailed rigid bodies in a continuum solvent in a periodic box [3]. The model used is based on that originally developed for the simulation of the diffusional association of two proteins and implemented in the SDA (Simulation of Diffusional Association) software [8]. For the simulation of many proteins, this method gives results in good agreement with experimental translational and rotational diffusion coefficients and small angle scattering structure factors for dilute [3] as well as concentrated protein solutions [12]. In this approach, intermolecular forces are computed as the sum of electrostatic interaction, electrostatic desolvation, nonpolar desolvation and softcore repulsion terms [3,8]. For computational efficiency, all these terms are precomputed on grids for each macromolecular solute before carrying out the BD simulations. To overcome errors due to the finite size of the electrostatic grids, we here describe the implementation of a longrange electrostatic correction to the model for interaction forces used in our BD simulations. The purpose of this correction is to improve the accuracy of the computed interprotein forces and extend the applicability of the approach to highly charged proteins and low ionic strength conditions. For validation, we performed BD simulations of bovine serum albumin (BSA) and hen egg white lysozyme (HEWL) with and without the longrange electrostatic correction and compared the results to experimentally determined small angle scattering structure factors and selfdiffusion coefficients. The same methodology described here for the implementation of the longrange DebyeHückel correction, should also be applicable in implicit solvent molecular dynamics simulations that make use of gridded interaction potentials [1316].
Methods
Brownian dynamics (BD) is a simulation method that employs a mesoscopic model in which the solvent is treated as a continuum and the solutes are modelled as discrete entities at a level of detail appropriate for the problem being studied. BD thus takes advantage of the large separation in timescale between the fast solvent motion and the slower motion of solute particles (polymers or colloids) which make it possible to treat the solvent implicitly. Furthermore, internal solute degrees of freedom are often neglected and macromolecules are treated as rigid bodies interacting by direct interactions (electrostatic, van der Waals, nonpolar) and solventmediated (hydrodynamic) interactions. Due to these simplifications, BD can be used to study larger biomacromolecular systems on longer timescales than is possible with classical atomicdetail molecular dynamics simulations.
Translational motion is propagated according to the following Equation [17]:
where r_{i }is the position of the center of geometry of the solute i and Δt = (t_{1 } t_{0}) is the timestep.
The effect of the solvent is described by a random displacement, R_{i}, which mimics the collision of solute i with the solvent molecules and is defined by a Gaussian distribution with mean 〈R_{i}〉= 0 and covariance . From the latter, it follows that the stochastic displacement is proportional to the square root of the translational diffusion tensor, . The second term on the r.h.s. of Equation 1, the divergence of the diffusion tensor, describes the hydrodynamic drift of the solute towards regions of high mobility. The force acting on solute i results from the sum of the forces acting on solutes j at time t_{0}, F_{j}(t_{0}), coupled with the diffusion tensor.
We employ a simplified treatment of hydrodynamic interactions to avoid the computationally expensive Cholesky factorization required to calculate the square root of the diffusion matrix. A mean field approach is used where is replaced by a volume fractiondependent diffusion coefficient, , and Equation 1 simplifies to [12]
We define the local volume, V_{i}, as the volume of the sphere of radius R^{cut }centered on solute i. The local volume fraction ϕ_{i }for the solute i is obtained by dividing the sum of the volumes of the solutes within R^{cut }by the local volume V_{i}[18]. The volume of a protein, v, is computed by approximating the protein as a sphere having a radius equal to the hydrodynamic radius (σ^{stokes}) estimated using HYDROPRO [19]. The cutoff for the local volume, R^{cut}, is set to four times the side of the largest interaction grid of the central solute. For a small simulation box, this cutoff was rescaled to a value equal to half of the simulation box size. A solute j is totally included in the local volume when the centertocenter distance d_{ij }between the central solute i and solute j is less than . When a solute k is only partially included within R^{cut}, that is, when , we account for that portion of solute volume derived by the spheresphere intersection. The volume fraction dependent shorttime translational diffusion coefficient () is then obtained using the Tokuyama model [2022], derived for a concentrated hardsphere suspension of particles interacting with both direct and hydrodynamic interactions. An equation analogous to Equation 2 is used for the rotational motion [12], with the volume fraction dependent shorttime rotational diffusion coefficient obtained using the model derived by Cichocki et al. which includes lubrication forces as well as two and threebody expansions of the mobility functions [23].
The forces, F_{i}, are computed as finitedifference derivatives of the pairwise free energies of interaction between the solutes as described in the next section.
Interaction energies and forces
For each pair of macromolecules, the interaction free energy, ΔG^{12}, is defined as:
A detailed description and parameterization of Equation 3 can be found in Refs. [3,24]. Briefly, the first two terms in Equation 3 are the interaction energies of the charges of one macromolecule ( or ) with the electrostatic potential of the other macromolecule ( or ). Charges were assigned using the effective charge approximation [25]. The third and fourth terms of Equation 3 represent the electrostatic desolvation energy arising from the introduction of the low dielectric cavity of one macromolecule in the presence of the charges of the other [25,26]. The desolvation energy is computed as the interaction of the charges of one macromolecule ( or ) with the electrostatic desolvation potential of the other macromolecule ( or ) [26], with parameterization as in Ref. [24]. The fifth and sixth terms in Equation 3 correspond to the nonpolar interactions due to the burial of the solvent accessible surface areas (SASAs) of the surface atoms The last two terms of Equation 3 describe the softcore repulsive potential introduced to avoid overlaps. The softcore potential is modelled using an inverse power function. The smoothness of the softcore potential allows abrupt changes in the forces at close contact to be avoided. In Equation 3, r specifies the atomic coordinates. For computational efficiency, all interaction potentials, Φ, are mapped onto grids centered on each of the macromolecules.
This formalism implies a truncation of the electrostatic potential in the gridcharge formalism due to the finite extent of the grids. To alleviate this problem, we here introduce an analytical longrange correction to the electrostatic interaction term that makes use of the assumption that beyond the electrostatic grid boundaries, a macromolecule can be treated as a DebyeHückel sphere.
According to the DebyeHückel theory of dilute electrolyte solutions, all ions in the solvent are treated as point charges while each pair of solutes are treated as spheres with radii a_{i}, a_{j }and net charges z_{i}e_{l}, z_{j}e_{l}, where e_{l }is the elementary charge. Then, the potential of mean force between a pair of solute molecules is
where ε_{0} is the vacuum permittivity, ε_{r }is the relative permittivity of the solvent, a = a_{i }+ a_{j}, and κ is the inverse of the Debye length, and is proportional to the ionic strength .
As shown in Equation 3, to compute the electrostatic interaction between a pair of macromolecules, the electrostatic potential of macromolecule 1 is multiplied by the effective charges of the second macromolecule. Due to the finite size of the grid, when the second macromolecule is on the border of the electrostatic potential grid of macromolecule 1, only a fraction of the effective charges on macromolecule 2 are taken into account for computing the electrostatic interaction. An isotropic distance cutoff from the center of macromolecule 1 is used in computing this interaction, so that if the effective charge is beyond this distance cutoff, its electrostatic interaction is not computed. The spherical cutoff is assigned on the assumption that the electrostatic potential becomes centrosymmetric at the grid edges and therefore a switch to the analytical DebyeHückel potential can be made beyond the cutoff. The application of the DebyeHückel potential reduces the discontinuity in the energy and forces at the grid cutoff distance.
Second osmotic virial coefficients
Osmotic virial coefficients are coefficients in the virial expansion of the state equation and they reflect deviations from ideal behaviour due to the presence of interactions. For simple cases, they can be obtained analytically. For this reason, they are commonly used to assess force field accuracy [1,3,27,28].
From classical statistical mechanics, the second osmotic virial coefficient can be obtained from [29]
Where r is the centertocenter distance and w (r) is the potential of mean force. For an isotropic potential, the corresponding equation is
Small angle scattering intensity
To assess the correctness of the interaction potentials, we compared experimental and computed small angle scattering intensities. Scattering intensities were computed from the simulations using [30]
where γ is a factor related to instrument effects, n_{p }= N/V is the protein concentration expressed as number density (N is the number of particles and V the total volume of the solution), Δ ρ is the electron density contrast between the scattering particle and the solvent, and v is the particle volume. P (q) is the form factor normalized such that P(0) = 1, S (q) is the structure factor and q is the scattering vector. The prefactor γ (Δρ)^{2}v^{2} can be obtained in experiments and then the normalized scattering intensity is expressed as
We computed the form factor for BSA using the analytical expression for the orientationally averaged form factor of an oblate ellipsoid with radii a and b where a is the semiaxis of revolution [31,32]. Following ref. [32], we set a = 17.5 Å and b = 47.4 Å.
The structure factor, S (q), was computed by Fourier transformation of the radial distribution function, g (r) [33] as follows
where n_{p} is the number density, r is the centertocenter distance, q is the magnitude of the scattering vector given by q=4πλ^{1}sin (θ/2) (where θ is the total scattering angle) and h (r) is the total correlation function which is given by h (r) = g(r)  1. The radial distribution function was computed from BD simulations using the centertocenter protein distances. We estimated the convergence of the g (r) by checking that it was not varying with increasing simulation time. This was done by computing the g (r) over the full trajectory and comparing this g (r) with an average g (r) computed from 20 segments selected sequentially from the trajectory.
Test systems of two spherical particles
For a system composed of two charged softsphere particles interacting via a DebyeHückel potential, the longrange contribution to the second virial coefficient can be computed by integrating Equation 6. This equation can be solved analytically by expanding the exponential up to the second order and substituting the DebyeHückel expression for the potential of mean force [29,34].
Only the longrange contribution to the second virial coefficient is taken into account in the analysis. Hence, the lower bound of the integration (lb) is not 0 but it is set to the sum of the protein radii (a_{i }+ a_{j}) plus one or two Debye lengths (1/κ). For example, solving Equation 5 setting the lower bound to lb = (a_{i }+ a_{j}) + 1/κ gives
where e is the base of the natural logarithm, e_{l }is the elementary charge and ρ is the concentration of the ions (equivalent to the ionic strength for monovalent ions).
The reason for considering only the longrange contribution is twofold. Firstly, our purpose is to assess the accuracy of the longrange DebyeHückel potential included in the BD simulation model. Secondly, for the expansion of the exponential e^{w/kT} up to the second order to be reasonably accurate, w/kT ≪ 1 is required. This means that the shortrange contribution of B_{22} at low ionic strength or for highly charged systems cannot be obtained using Equation 5.
In the numerical integration, the two particles were represented by spherical fullerenelike particles of radius 6 Å composed of 180 atoms. A partial point charge was placed on each atom. The total charge of each sphere was uniformly distributed over all the atoms. Different systems were simulated by varying the net charge and the ionic strength (see Table 1 and Table 2 in Results and discussion). The interaction energy between the two particles is given by
Table 1. Long range contribution to theB_{22} value at 5 mM ionic strength for the two softsphere systems
Table 2. Long range contribution to the B_{22 }values at 300 mM ionic strength for the two softsphere systems
To compute the second virial coefficient, one particle was kept fixed at the center of the simulation box and the other was moved on a regular lattice within the simulation box, avoiding overlaps with the central particle. The size of the box was set to 400×400×400 Å ^{3} and the dimension of the lattice was set to 100×100×100 vertices. The interaction energy (Equation 11) was computed for each position assumed by the second particle and the second virial coefficient was computed by integrating Equation 6 numerically with the potential of mean force, , where r is the centertocenter distance. As for the analytical computation of B_{22}, the integration was performed setting half, one, or two Debye lengths as the lower bound of the integral.
We considered two spherical particles i and j with corresponding radii a_{i} and a_{j} and net charges z_{i} and z_{j}, each resulting from 180 partial point charges uniformly distributed near the surface of each particle at a distance r from the particle’s center. Six different combinations of net charges on the particles were tested, namely: + 1/ + 1, + 5/ + 5, + 10/ + 10 and + 1/ 1, + 5/ 5, + 10/ 10 (in units of elementary charge). For each pair of particles the integration was performed at different ionic strengths, 5 mM and 300 mM. These two ionic strengths were chosen to assess the importance of the DebyeHückel term at low and high salt conditions (compared to the 150 mM physiological ionic strength). The computed values were obtained by with and without inclusion of the DebyeHückel potential.
From the set of approximately 10^{6} interaction energies computed at the lattice vertices (avoiding the overlapping region), we extracted 100 random subsets of 10^{5} values. For each subset, the second virial coefficient was computed. Then, an average B_{22} and a standard deviation over the subset was calculated.
BD Simulations of protein solutions
BD simulations were performed with SDAMM [3], a parallelized program based on the SDA software [8] capable of handling many proteins (10^{3} 10^{4}) treated as rigid bodies in atomic detail. For further details, see [3].
BD simulations were carried out for 250 protein molecules that were initially randomly positioned (avoiding overlaps) in a cubic box with periodic boundary conditions. The dimensions of the simulation box were varied according to the concentration of the protein solution.
The DebyeHückel interaction between a pair of proteins was computed up to a distance cutoff of 4 times the side of the electrostatic grid. If the simulation box was small, to avoid selfimage interactions, this cutoff was rescaled to a value equal to half of the simulation box size.
Each system was subjected to 5 or 10 μs of simulation at 300 K. Equilibration was assessed by monitoring the convergence of the radial distribution function and the stabilization of the energies. In all cases, 1 μs was sufficient to obtain an equilibrated system according to these criteria and the remaining 4 or 9 μs were used for the analysis. The integration timestep was 0.5 ps. The positions and orientations of the proteins were recorded along with energy values every 0.5 ns.
Simulations of HEWL were performed at 14, 28, 57 and 85 g/L for comparison with experimental longtime translational selfdiffusion coefficients [35]. Four sets of simulations were performed varying the ionic strength (1 mM and 5 mM) and including or omitting the analytical DebyeHückel potential. Simulations were performed for 5 μs.
Simulations of BSA were performed at 0.9, 4.5, 9, 18, 45, 90 g/L for comparison with the experimental small angle Xray scattering (SAXS) intensities described in ref. [32]. Two sets of simulations were performed. In one set, the DebyeHückel potential was included, whereas in the other set, the DebyeHückel potential was omitted. Because of the faster convergence of the higher concentration simulations, simulations at 0.9, 4.5, 9 and 18 g/L were performed for 10 μs whereas the simulations at 45 and 90 g/L were performed for 5 μs.
Protein preparation
The crystal structure of hen egg white lysozyme (HEWL) was taken from the Protein Data Bank (ref): 1hel. The structure of BSA used for the simulations was a model taken from Modbase [36]. It was obtained by homology modelling based on the crystal structure of human serum albumin (HSA) [37].
Polar hydrogen atoms were added to the structures according to the specified pH and ionic strength (IS) using the H++ software [38]. The simulations of HEWL were performed at pH 5 ; the computed net charge of HEWL was +10e. The simulations of BSA were performed at pH 7. BSA had a computed net charge of 16e.
Atomic partial charges and radii were assigned to all the atoms from the OPLS united atom force field [39]. Electrostatic potential grids Φ were computed by solving the linearized PoissonBoltzmann equation using the program UHBD [40]. The grid size was set to 100×100×100 Å ^{3} for HEWL and 200×200×200 Å ^{3} for BSA with a grid spacing of 1.0 Å. Nonpolar desolvation, electrostatic desolvation and softcore repulsion grids were set to 100×100×100 Å ^{3} for HEWL and 130×130×130 Å ^{3} for BSA, with a grid spacing of 1.0 Å.
Results and discussion
Comparison of simulations and analytical results for systems of two spherical particles
The two spheres system (see Computational Details section) was simulated with different combinations of net solute charge at two ionic strengths with and without inclusion of the DebyeHückel potential. For each system, the analytical value of the long range contribution to the B_{22} was compared to the computed one. All values are given in Table 1 for 5 mM and Table 2 for 300 mM ionic strength. For a better comprehension of the length scale of the contribution of the electrostatic potential to the second virial coefficient, the analytical B_{22} values from the analytical calculations and from the simulations were obtained using different lower bounds for integrating Equation 6. We first consider the systems at low ionic strength (5 mM).
5 mM ionic strength
Let us first consider the integration done with a lower bound of one Debye length which at 5 mM ionic strength corresponds to 43 Å. From Table 1, it is clear that when using a grid of 100×100×100 Å ^{3} without the DebyeHückel potential, the longrange decay of the electrostatic potential is not captured. This result is expected since the electrostatic potential grid size is of the same order as the Debye length. Doubling the length of the side of the grid results in a B_{22} value that is approximately 50% of the analytical value. The longrange tail (beyond 100 Å) of the electrostatic potential is missing and it is apparent that it represents an important contribution to the second virial coefficient.
By turning on the DebyeHückel potential and keeping the smaller electrostatic potential grid (side length: 100 Å), more than the 90% of the analytical B_{22} value is recovered. For systems with the highest net charge at one Debye length, the potential is too high and the integral expression in Equation 6 diverges.
For a perfectly isotropic case, such as this one, the DebyeHückel potential smoothly recovers the truncation of the electrostatic potential due to the finite grid. This can be seen from the electrostatic potential energy computed by varying the interparticle separation (see Additional file 1).
Additional file 1. Charged spheres electrostatic potential energy. Electrostatic potential energy of two uniformly charged spheres for different net charge (e) combinations. Panels A,B: +1/+1, +1/1; panels C,D: +5/+5, +5/5; panels E,F: +10/+10, +10/10. Red and green colors show interactions between charges of the same and opposite sign, respectively. The DebyeHückel analytical approximation, Equation 4 (continuous line), Brownian dynamics without DebyeHückel term (crosses) and Brownian dynamics with DebyeHückel term (circles) are shown. Interactions are computed at 5 mM (left panels) and 300 mM (right panels). The abscissa is the particle centertocenter distance divided by the particle diameter. The inclusion of the DebyeHückel potential recovers the dependence of the energy on particle separation computed with the analytical model.
Format: PDF Size: 2.1MB Download file
This file can be viewed with: Adobe Acrobat Reader
At two Debye lengths (2/κ), the B_{22} value of the systems with the smaller grid (100 Å) without the DebyeHückel potential is zero, since the grid is smaller than the Debye length. By doubling the grid dimension, the side of the grid becomes of the same order as the Debye length and the B_{22} is still not computed correctly. With the DebyeHückel potential and the smaller grid, however, the analytical second virial coefficient can be reproduced well.
300 mM ionic strength
Increasing the ionic strength to 300 mM, at lower bounds of one or two Debye lengths (5.5 Å), the B_{22} values computed using only the smaller electrostatic potential grid agree rather well with the analytical values, see Table 2. Doubling the grid dimensions or adding the DebyeHückel potential is not required because more than 90% of the interactions are captured within one Debye length. It is clear that at 300 mM ionic strength, the gridbased formalism is sufficient to properly describe the longrange electrostatic interaction, even using the smaller grid.
Protein systems modeled in atomic detail
We now turn to more complex and realistic systems composed of solutions of proteins represented in atomic detail subjected to BD simulation as described in the Computational Details section.
Scattering intensities Several BSA solutions at different concentrations were simulated for 10 μs to 20 μs using BD. To assess the effect of the DebyeHückel approximation on the BSA selfinteractions, two sets of simulations were performed. In one set, the DebyeHückel potential was included whereas in the other set, it was omitted.
Normalized small angle scattering intensities were computed using Equation 8 and compared to experimental SAXS intensities. The experiments were performed without added salt which corresponds to an ionic strength up to 5 mM [31,32]. This nonzero ionic strength arises from several factors such as dissolved CO _{2}, a residual amount of salt present in the protein solution, and the dissociation of surface groups upon solvation [31,32]. Simulations were performed at 5 mM ionic strength with a corresponding Debye length of 43.1 Å.
As shown in Figure 1, the scattering intensities obtained from the simulations with the DebyeHückel approximation reproduce experimental SAXS intensities better then the intensities calculated from simulations which do not include the DebyeHückel interaction. In particular, the greatest improvement is seen at low q values, i.e. long range interactions are accurately captured. At high concentrations, the DebyeHückel approximation tends to overestimate the height of the correlation peak seen in the normalized experimental intensities. This phenomenon can be explained considering that simulations have been performed at 5 mM ionic strength, but that at high protein concentrations, the effective ionic strength may be higher due to the presence of highly charged proteins. Indeed, the correlation peak is lower in the simulations without the DebyeHückel approximation (see also Figure 2 and Figure 3). This suggests that at low ionic strength and high protein concentration, the ionic strength of the simulation should be slightly increased to better reproduce experimentally observed scattering intensities.
Figure 1. BSA SAS intensities. Experimental [32] (dashed lines) and computed (continuous lines) normalized small angle scattering intensities at different concentrations (indicated on the plots) of BSA. Computed curves from simulations without (A) and with (B) the DebyeHückel approximation. Curves are shifted by 0.2 on the vertical axis for better visibility.
Figure 2. BSA structure factors. Experimental [32] (dashed lines) and computed (continuous lines) structure factors at various concentrations (indicated on the plot) of BSA obtained from simulations without (dark green) and with (dark red) the DebyeHückel approximation. Curves are shifted by 0.2 on the vertical axis for better visibility.
Figure 3. BSA radial distribution functions. Computed radial distribution functions at various concentrations (indicated on the plot) of BSA obtained from simulations without (dark green) and with (dark red) the DebyeHückel approximation. Curves are shifted by 0.2 on the vertical axis for better visibility. Averages and standard deviations of the g(r) are shown by the dark line and light color, respectively.
The computed static structure factors obtained from the two sets of simulations are compared in Figure 2. Focusing on the low q region (q < 0.1 nm^{1}), for a given concentration, the value of S(q) is lower when the DebyeHückel potential is used. The long wavelength limit of S (q) is proportional to the normalized isothermal osmotic compressibility, vis.:
where χ_{T }is the isothermal osmotic compressibility. (In the canonical ensemble, ), n_{p} is the protein number density, and k_{B} is the Boltzmann constant [32,41,42]. The decrease of S (q) at low q values can be explained by the decrease of the osmotic compressibility due to the longrange electrostatic repulsion introduced with the DebyeHückel potential [43].
The first peak in the S(q) represents the correlation between a pair of proteins. We observe that the simulations which include the DebyeHückel potential show a shift of the first peak to lower q values (at high concentrations) or the appearance of a peak (at low concentrations), indicating the presence of a longrange correlation between the proteins. With increasing concentration, the peak shifts to higher q values, suggesting a decrease of the correlation distance. The same effect can be seen better in real space from the radial distribution functions plotted in Figure 3 where it can be seen that the introduction of a longrange repulsion pushes the proteins away from each other. It also leads to a more structured solution, with the appearance of a second peak in the simulations at 90 g/L protein concentration.
Longtime selfdiffusion coefficients Besides the effect on proteinprotein interactions, the addition of the DebyeHückel potential also has consequences for the dynamics of the proteins. Simulations of HEWL were performed at low ionic strength (1 and 5 mM) at different lysozyme concentrations and compared to experimental diffusion coefficients obtained from pulsed gradient spin echo NMR for HEWL solutions without added salt at pH 4.9. As shown in Figure 4, the presence of the DebyeHückel potential systematically lowers the longtime selfdiffusion coefficients. This effect can be explained considering that, for a given concentration, simulations which include the DebyeHückel potential correspond to a larger effective concentration due to the longrange repulsive interaction [43,44]. In general, the magnitude of the effect on the diffusion coefficient due to the DebyeHückel potential is related to the ionic strength of the solution, the size of the protein, and the protein concentration. For proteins whose size is comparable to the Debye length, κ^{1}, as in our case, this effect can be significant. For very large proteins, the Debye length can be much smaller than the size of the protein, and hence, adding the longrange DebyeHückel interaction may lead only to small effects on the diffusion coefficient.Simulations performed at 1 mM ionic strength underestimate the diffusion coefficients compared to the experimental values (see Figure 4). As described above for the BSA case, the ionic strength of the solution is affected by several factors. Thus, it is possible that the value of 1 mM used in the simulations does not correctly describe the effective ionic strength of the experimental solutions. We therefore also performed simulations at higher ionic strength (5 mM), obtaining better agreement with the experimental data, see Figure 4.
Figure 4. HEWL translational diffusion coefficients. Normalized longtime translational selfdiffusion coefficients of HEWL at low ionic strength. Simulations were performed at 1 mM (A) and 5 mM (B) ionic strength. Experimental values from ref. [35] (black diamonds), and computed values from BD simulations with (red squares) and without (green squares) DebyeHückel potential are shown. The Tokuyama [22] analytical model is shown by the black dotted line. Insets are loglog plots of the same data.
Methodological considerations
The DebyeHückel potential has been implemented together with cubic grids for the proteins. The transition from the gridded potential to the DebyeHückel potential with increasing distance from a solute center occurs at the shortest distance to the grid boundary. Thus, cubic grids permit the most efficient implementation of the DebyeHückel correction. Their use is usually appropriate for globular proteins, however, it may become an issue when modeling large elongated molecules. For the latter, a large number of grid points on a cubic grid will have very low (negligible) values of the mapped interaction potentials, leading to an unnecessarily high memory requirement.
On the other hand, an advantage of the DebyeHückel implementation is that it removes the requirement for the electrostatic potential to have very small values at the grid edges; the electrostatic potential is only required to be centrosymmetric. This means that smaller grids can be used with the longrange interactions being captured by the DebyeHückel with only a small computational cost (see Additional file 2).
Additional file 2. Runtime versus protein concentration for the simulations of BSA. The DebyeHückel correction requires very little additional computational effort. Indeed, at low concentrations, the inclusion of the DebyeHückel correction keeps the likecharged molecules apart and therefore reduces the number of pairs of molecules for which gridtype potentials are used to compute intermolecular forces, thus leading to reduced runtimes.
Format: PDF Size: 604KB Download file
This file can be viewed with: Adobe Acrobat Reader
Using the DebyeHückel correction may be an issue for some highly or nonuniformly charged systems as it can lead to force discontinuities at the grid boundaries. A possible solution to this problem, currently not implemented, is to apply an interpolation function between the electrostatic potential grid and the DebyeHückel potential for computing the forces at the grid boundary.
Conclusions
We have here described the implementation of a DebyeHückel correction for the computation of gridbased electrostatic interaction energies and forces for use in atomically detailed manyprotein Brownian dynamics simulations. The ability of this manyprotein BD method to correctly reproduce small angle scattering data and diffusion coefficients, was previously shown for several proteins [3,12]. Due to computational limitations on the size of the electrostatic interaction grids, the method could not be applied to highly charged systems or low ionic strength conditions without impairing the accuracy of the resulting simulations. The introduction of the simple DebyeHückel correction described in this paper with its very low associated computational costs allowed us to extend the range applicability of this BD method to highly charged systems at low ionic strength. In particular, comparison of the model with the DebyeHückel correction to analytical results for spherical solutes, as well as to experimental SAXS intensities for BSA protein solutions, and to longtime selfdiffusion coefficients of HEWL protein solutions, showed good agreement. Some other potential applications of the methodology are the simulation of protein crystallization, of proteinsurface adsorption, and of heterogeneous crowded protein solutions. Furthermore, the DebyeHückel correction described here should be of value in implicit solvent molecular dynamics simulations which make use of gridded interaction potentials [1316].
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
PM and RCW devised and planned the project. PM and MM implemented the method. PM carried out the simulations. PM, MM and RCW analyzed the data. PM and RCW wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
Financial support from the Center for Modelling and Simulation in the Biosciences (BIOMS) (P.M.), the German Ministry for Education and Research (BMBF) Virtual Liver Network (Grant No. 0315749), and the Klaus Tschira Foundation, and a grant for supercomputing time at the Environmental Molecular Science Laboratory (Grant No. 30994) are gratefully acknowledged. We would like to gratefully acknowledge Susanne Krömker for helpful discussion and Stefan Richter for his help in code development.
References

McGuffee SR, Elcock AH: Atomically detailed simulations of concentrated protein solutions: The effects of salt, pH, point mutations, and protein concentration in simulations of 1000molecule systems.

Geyer T, Winter U: An O(N^{2}) approximation for hydrodynamic interactions in Brownian dynamics simulations.

Mereghetti P, Gabdoulline RR, Wade RC: Brownian dynamics simulation of protein solutions: structural and dynamical properties.

Ando T, Skolnick J: Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion.

McGuffee SR, Elcock AH: Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm.

Wieczorek G, Zielenkiewicz P: Influence of macromolecular crowding on proteinprotein association rates  a Brownian dynamics study.

Goodsell DS, Olson AJ: Automated docking of substrates to proteins by simulated annealing.

Gabdoulline RR, Wade RC: Simulation of the diffusional association of barnase and barstar.

Goodford PJ: A computational procedure for determining energetically favorable binding sites on biologically important macromolecules.

Wu X, Milne J, Borgnia M, Rostapshov A, Subramaniam S, BR B: A coreweighted fitting method for docking atomic structures into lowresolution maps: application to cryoelectron microscopy.

Sherwood P, Brooks BR, Sansom MSP: Multiscale methods for macromolecular simulations.

Mereghetti P, Wade RC: Atomic detail Brownian dynamics simulations of concentrated protein solutions with a mean field treatment of hydrodynamic interactions.

Sharp K: Incorporating solvent and ion screening into molecular dynamics using the finitediffernece PoissonBoltzmann method.

Lu BZ, Chen WZ, Wang CX, Xu Xj: Protein molecular dynamics with electrostatic force entirely determined by a single PoissonBoltzmann calculation.

Prabhu NV, Zhu P, Sharp KA: Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smoothpermittivity finite difference PoissonBoltzmann method.

Wang J, Tan C, Lu Q, Luo R, Tan Yh: PoissonBoltzmann Solvents in Molecular Dynamics Simulations.

Ermak DL, McCammon JA: Brownian dynamics with hydrodynamic interactions.

UrbinaVillalba G, GarciaSucre M, ToroMendoza J: Average hydrodynamic correction for the Brownian dynamics calculation of flocculation rates in concentrated dispersions?

de la Torre JG, Huertas ML, Carrasco B: Calculation of hydrodynamic properties of globular proteins from their atomiclevel structure.

Tokuyama M, Oppenheim I: Dynamics of hardsphere suspensions.

Tokuyama M, Oppenheim I: On the theory of concentrated hardsphere suspensions.

Tokuyama M, Moriki T, Kimura Y: Selfdiffusion of biomolecules in solution.

Cichocki B, EkielJezewska ML, Wajnryb E: Lubrication corrections for threeparticle contribution to shorttime selfdiffusion coefficients in colloidal dispersions.

Gabdoulline RR, Wade RC: On the contributions of diffusion and thermal activation to electron transfer between Phormidium laminosum plastocyanin and cytochrome f: Brownian dynamics simulations with explicit modeling of nonpolar desolvation interactions and electron transfer events.

Gabdoulline RR, Wade RC: Effective charges for macromolecules in solvent.

Elcock A, Gabdoulline R, Wade R, McCammon J: Computer simulation of proteinprotein association kinetics: Acetylcholinesterase Fasciculin.

Neal BL, Asthagiri D, Lenhoff AM: Molecular origins of osmotic second virial coefficients of proteins.

Elcock AH, McCammon JA: Calculation of weak proteinprotein interactions: the pH dependence of the second Virial coefficient.

Hill TL: An introduction to statistical thermodynamics. New York: Dover Publications Inc; 1986.

da Silva MA, Itri R, Arêas EPG: Lysozyme viscoelastic matrices in tetramethylurea/water media: a small angle Xray scattering study.

Zhang F, Jacobs RMJ, Martin CM, Schreiber F, Skoda MWa: Protein interactions studied by SAXS: effect of ionic strength and protein concentration for BSA in aqueous solutions.

Heinen M, Zanini F, RoosenRunge F, Zhang F, Hennig M, Seydel T, Schweins R, Sztucki M, Antalík M, Schreiber F, Nägele G, Fedunová D: Viscosity and diffusion: crowding and salt effects in protein solutions.

Allen MP, Tildesley DJ: Computer Simulations of Liquids. USA: Oxford University Press; 1989.

Asthagiri D, Paliwal A, Abras D, Lenhoff AM, Paulaitis ME: A consistent experimental and modeling approach to lightscattering studies of proteinprotein interactions in solution.

Price WS, Tsuchiya F, Arata Y: Lysozyme aggregation and solution properties studied using PGSE NMR diffusion measurements.

Pieper U, Eswar N, Webb BM, Eramian D, Kelly L, Barkan DT, Carter H, Mankoo P, Karchin R, MartiRenom MA, Davis FP, Sali A: modbase, a database of annotated comparative protein structure models and associated resources.

Balbo J, Mereghetti P, Herten DP, Wade R: The shape of protein crowders is a major determinant of protein diffusion.

Gordon JC, Myers JB, Folta T, Shoja V, Heath LS, Onufriev A: H++: a server for estimating pKas and adding missing hydrogens to macromolecules.

Jorgensen WL, TiradoRives J: The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin.

Madura JD, Briggs JM, Wade RC, Davis ME, Luty BA, Andrew I, Antosiewicz J, Gilsong MK, Bagheri B, Scott LR, McCammon JA: Electrostatics and diffusion of molecules in solution: simulations with the University of Houston Brownian Dynamics program.

Balescu RC: Equilibrium and NonEquilibrium Statistical Mechanics. New York: John Wiley & Sons; 1975.

Nägele G: On the dynamics and structure of chargestabilized suspensions.

Fukasawa T, Sato T: Versatile application of indirect Fourier transformation to structure factor analysis: from Xray diffraction of molecular liquids to small angle scattering of protein solutions.

Stylianopoulos T, Poh MZ, Insin N, Bawendi MG, Fukumura D, Munn LL, Jain RK: Diffusion of particles in the extracellular matrix: the effect of repulsive electrostatic interactions.