Long range Debye-Hückel correction for computation of grid-based electrostatic forces between biomacromolecules
- Paolo Mereghetti^{1, 2},
- Michael Martinez^{1} and
- Rebecca C Wade^{1, 3}Email author
DOI: 10.1186/2046-1682-7-4
© Mereghetti et al.; licensee BioMed Central Ltd. 2014
Received: 11 February 2014
Accepted: 4 June 2014
Published: 17 June 2014
Abstract
Background
Brownian dynamics (BD) simulations can be used to study very large molecular systems, such as models of the intracellular environment, using atomic-detail 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 three-dimensional discretized grids. For long-range interactions, such as electrostatics, grid-based methods are subject to finite size errors. We describe here the implementation of a Debye-Hückel correction to the grid-based 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 long-range electrostatic correction increased the accuracy of both the protein-protein 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 long-range 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 Debye-Hückel Poisson-Boltzmann 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 [1–5]. 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 three-dimensional 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 long-range 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 short-range and a coarse grid with a larger grid spacing for the long-range 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 short-range 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, non-polar desolvation and soft-core 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 long-range 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 inter-protein 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 long-range electrostatic correction and compared the results to experimentally determined small angle scattering structure factors and self-diffusion coefficients. The same methodology described here for the implementation of the long-range Debye-Hückel correction, should also be applicable in implicit solvent molecular dynamics simulations that make use of gridded interaction potentials [13–16].
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, non-polar) and solvent-mediated (hydrodynamic) interactions. Due to these simplifications, BD can be used to study larger biomacromolecular systems on longer time-scales than is possible with classical atomic-detail molecular dynamics simulations.
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 $\u3008{\mathbf{R}}_{i}{\mathbf{R}}_{j}\u3009=2{D}_{\mathit{\text{ij}}}^{t}\mathrm{\Delta t}$. From the latter, it follows that the stochastic displacement is proportional to the square root of the translational diffusion tensor, ${D}_{\mathit{\text{ij}}}^{t}$. 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 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 center-to-center distance d_{ ij }between the central solute i and solute j is less than ${R}_{\mathit{\text{cut}}}-{\sigma}_{j}^{\mathit{\text{stokes}}}$. When a solute k is only partially included within R^{ cut }, that is, when ${R}^{\mathit{\text{cut}}}-{\sigma}_{k}^{\mathit{\text{stokes}}}<{d}_{\mathit{\text{ik}}}<{R}^{\mathit{\text{cut}}}+{\sigma}_{k}^{\mathit{\text{stokes}}}$, we account for that portion of solute volume derived by the sphere-sphere intersection. The volume fraction dependent short-time translational diffusion coefficient (${D}^{{t}_{\mathit{\text{short}}}}\left({\varphi}_{i}\right)$) is then obtained using the Tokuyama model [20–22], derived for a concentrated hard-sphere 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 short-time rotational diffusion coefficient obtained using the model derived by Cichocki et al. which includes lubrication forces as well as two- and three-body expansions of the mobility functions [23].
The forces, F_{ i }, are computed as finite-difference derivatives of the pairwise free energies of interaction between the solutes as described in the next section.
Interaction energies and forces
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 (${q}_{{i}_{2}}$ or ${q}_{{j}_{1}}$) with the electrostatic potential of the other macromolecule (${\Phi}_{e{l}_{1}}$ or ${\Phi}_{e{l}_{2}}$). 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 (${q}_{{i}_{2}}$ or ${q}_{{j}_{1}}$) with the electrostatic desolvation potential of the other macromolecule (${\Phi}_{\mathit{\text{edesol}}{v}_{1}}$ or ${\Phi}_{\mathit{\text{edesol}}{v}_{2}}$) [26], with parameterization as in Ref. [24]. The fifth and sixth terms in Equation 3 correspond to the non-polar 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 soft-core repulsive potential introduced to avoid overlaps. The soft-core potential is modelled using an inverse power function. The smoothness of the soft-core 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 grid-charge formalism due to the finite extent of the grids. To alleviate this problem, we here introduce an analytical long-range 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 Debye-Hückel sphere.
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 $\left({\kappa}^{2}=\frac{{e}_{l}^{2}\beta}{{\epsilon}_{0}{\epsilon}_{r}}\sum _{i}{\rho}_{i}{z}_{i}^{2}\right)$.
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 cut-off 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 cut-off is assigned on the assumption that the electrostatic potential becomes centrosymmetric at the grid edges and therefore a switch to the analytical Debye-Hückel potential can be made beyond the cutoff. The application of the Debye-Hückel potential reduces the discontinuity in the energy and forces at the grid cut-off 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].
Small angle scattering intensity
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 semi-axis of revolution [31, 32]. Following ref. [32], we set a = 17.5 Å and b = 47.4 Å.
where n_{ p } is the number density, r is the center-to-center 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 center-to-center 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 soft-sphere particles interacting via a Debye-Hückel potential, the long-range contribution to the second virial coefficient can be computed by integrating Equation 6. This equation can be solved analytically by expanding the exponential ${e}^{-w\left(r\right)/{k}_{B}T}$ up to the second order and substituting the Debye-Hückel expression for the potential of mean force [29, 34].
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 long-range contribution is two-fold. Firstly, our purpose is to assess the accuracy of the long-range Debye-Hückel potential included in the BD simulation model. Secondly, for the expansion of the exponential e^{-w/k T} up to the second order to be reasonably accurate, |w/k T| ≪ 1 is required. This means that the short-range contribution of B_{22} at low ionic strength or for highly charged systems cannot be obtained using Equation 5.
Long range contribution to the B _{ 22 } value at 5 mM ionic strength for the two soft-sphere systems
z _{ i } | z _{ j } | ${B}_{22}^{\left(E{P}_{100A}\right)}$ | ${B}_{22}^{\left(E{P}_{200A}\right)}$ | ${B}_{22}^{(E{P}_{100A}+\mathit{\text{DHP}})}$ | ${B}_{22}^{\left(A\right)}$ | ||||
---|---|---|---|---|---|---|---|---|---|
(e_{ l }) | (e_{ l }) | 1/κ | 2/κ | 1/κ | 2/κ | 1/κ | 2/κ | 1/κ | 2/κ |
+1 | +1 | 0.0 | 0.0 | 13.0^{(0.0)} | 0.0 | 29.0^{(0.0)} | 14.1^{(0.0)} | 30.4 | 16.1 |
+1 | -1 | 0.0 | 0.0 | -13.3^{(0.0)} | 0.0 | -29.4^{(0.0)} | -14.2^{(0.0)} | -30.7 | -16.2 |
+5 | +5 | 0.0 | 0.0 | 261.5^{(0.1)} | 0.0 | 636.9^{(0.1)} | 339.3^{(0.2)} | 675.8 | 397.4 |
+5 | -5 | 0.0 | 0.0 | -432.7^{(0.3)} | 0.0 | -853.8^{(0.3)} | -367.8^{(0.3)} | -878.0 | -429.8 |
+10 | +10 | 0.0 | 0.0 | 603.6^{(0.2)} | 0.0 | 1916.9^{(0.2)} | 1216.5^{(0.7)} | 1490.0 | 1428.1 |
+10 | -10 | 0.0 | 0.0 | -5037.6^{(9.7)} | 0.0 | -7338.9^{(10.5)} | -1682.2^{(1.7)} | -4738.0 | -1867.6 |
Long range contribution to the B _{ 22 } values at 300 mM ionic strength for the two soft-sphere systems
z _{ i } | z _{ j } | ${B}_{22}^{\left(E{P}_{100A}\right)}$ | ${B}_{22}^{\left(E{P}_{200A}\right)}$ | ${B}_{22}^{(E{P}_{100A}+\mathit{\text{DHP}})}$ | ${B}_{22}^{\left(A\right)}$ | ||||
---|---|---|---|---|---|---|---|---|---|
(e_{ l }) | (e _{ l } ) | 1/κ | 2/κ | 1/κ | 2/κ | 1/κ | 2/κ | 1/κ | 2/κ |
+1 | +1 | 0.3^{(0.0)} | 0.1^{(0.0)} | 0.3^{(0.0)} | 0.1^{(0.0)} | 0.3^{(0.0)} | 0.1^{(0.0)} | 0.3 | 0.1 |
+1 | -1 | -0.3^{(0.0)} | -0.1^{(0.0)} | -0.3^{(0.0)} | 0.1^{(0.0)} | 0.3^{(0.0)} | 0.1^{(0.0)} | -0.3 | -0.1 |
+5 | +5 | 6.4^{(0.0)} | 3.3^{(0.0)} | 6.6^{(0.1)} | 3.5^{(0.0)} | 6.7^{(0.0)} | 3.6^{(0.0)} | 7.4 | 4.0 |
+5 | -5 | -9.4^{(0.0)} | -3.7^{(0.0)} | -9.6^{(0.2)} | -3.9^{(0.1)} | -9.7^{(0.1)} | 4.0^{(0.1)} | -11.6 | -4.6 |
+10 | +10 | 18.1^{(0.0)} | 11.4^{(0.0)} | 18.8^{(0.1)} | 12.0^{(0.1)} | 19.4^{(0.1)} | 12.6^{(0.1)} | 4.2 | 12.8 |
+10 | -10 | -97.6^{(0.7)} | -18.4^{(0.1)} | -98.1^{(0.4)} | -19.0^{(0.1)} | -98.9^{(0.7)} | -19.6^{(0.1)} | -72.4 | -22.0 |
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, $w\left(r\right)=\Delta {G}_{\mathit{\text{Debye}}}^{1-2}$, where r is the center-to-center 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 Debye-Hü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 Debye-Hü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 Debye-Hü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 self-image 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 long-time translational self-diffusion coefficients [35]. Four sets of simulations were performed varying the ionic strength (1 mM and 5 mM) and including or omitting the analytical Debye-Hü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 X-ray scattering (SAXS) intensities described in ref. [32]. Two sets of simulations were performed. In one set, the Debye-Hückel potential was included, whereas in the other set, the Debye-Hü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 Poisson-Boltzmann 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 Å. Non-polar desolvation, electrostatic desolvation and soft-core 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 Debye-Hü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 Debye-Hückel potential, the long-range 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 long-range 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 Debye-Hü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 Debye-Hü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 inter-particle separation (see Additional file 1).
At two Debye lengths (2/κ), the B_{22} value of the systems with the smaller grid (100 Å) without the Debye-Hü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 Debye-Hü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 Debye-Hü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 grid-based formalism is sufficient to properly describe the long-range 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 Debye-Hückel approximation on the BSA self-interactions, two sets of simulations were performed. In one set, the Debye-Hü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 non-zero 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 Å.
where χ_{ T }is the isothermal osmotic compressibility. (In the canonical ensemble, ${\chi}_{T}=-V{\left(\frac{\partial V}{\partial \Pi}\right)}_{T}={\left\{{n}_{p}{\left(\frac{\partial \Pi}{\partial {n}_{p}}\right)}_{T}\right\}}^{-1}$), 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 long-range electrostatic repulsion introduced with the Debye-Hü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 Debye-Hü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 long-range 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 long-range 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.
Long-time self-diffusion coefficients
Methodological considerations
The Debye-Hückel potential has been implemented together with cubic grids for the proteins. The transition from the gridded potential to the Debye-Hü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 Debye-Hü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 Debye-Hü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 long-range interactions being captured by the Debye-Hückel with only a small computational cost (see Additional file 2).
Using the Debye-Hückel correction may be an issue for some highly or non-uniformly 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 Debye-Hückel potential for computing the forces at the grid boundary.
Conclusions
We have here described the implementation of a Debye-Hückel correction for the computation of grid-based electrostatic interaction energies and forces for use in atomically detailed many-protein Brownian dynamics simulations. The ability of this many-protein 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 Debye-Hü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 Debye-Hückel correction to analytical results for spherical solutes, as well as to experimental SAXS intensities for BSA protein solutions, and to long-time self-diffusion coefficients of HEWL protein solutions, showed good agreement. Some other potential applications of the methodology are the simulation of protein crystallization, of protein-surface adsorption, and of heterogeneous crowded protein solutions. Furthermore, the Debye-Hückel correction described here should be of value in implicit solvent molecular dynamics simulations which make use of gridded interaction potentials [13–16].
Declarations
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.
Authors’ Affiliations
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 1000-molecule systems. J Am Chem Soc. 2006, 128: 12098-12110.View ArticleGoogle Scholar
- Geyer T, Winter U: An O(N^{2}) approximation for hydrodynamic interactions in Brownian dynamics simulations. J Chem Phys. 2009, 130: 114905-View ArticleGoogle Scholar
- Mereghetti P, Gabdoulline RR, Wade RC: Brownian dynamics simulation of protein solutions: structural and dynamical properties. Biophys J. 2010, 99 (11): 3782-91.View ArticleGoogle Scholar
- Ando T, Skolnick J: Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion. Proc Natl Acad Sci. 2010, 107 (43): 18457-18462.View ArticleGoogle Scholar
- McGuffee SR, Elcock AH: Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm. PLoS Comput Biol. 2010, 6 (3): e1000694-View ArticleGoogle Scholar
- Wieczorek G, Zielenkiewicz P: Influence of macromolecular crowding on protein-protein association rates - a Brownian dynamics study. Biophys J. 2008, 95: 5030-5036.View ArticleGoogle Scholar
- Goodsell DS, Olson AJ: Automated docking of substrates to proteins by simulated annealing. Proteins: Struct Function Genet. 1990, 8 (8): 195-202.View ArticleGoogle Scholar
- Gabdoulline RR, Wade RC: Simulation of the diffusional association of barnase and barstar. Biophys J. 1997, 72 (5): 1917-1929.View ArticleGoogle Scholar
- Goodford PJ: A computational procedure for determining energetically favorable binding sites on biologically important macromolecules. J Med Chem. 1985, 28 (7): 849-857.View ArticleGoogle Scholar
- Wu X, Milne J, Borgnia M, Rostapshov A, Subramaniam S, BR B: A core-weighted fitting method for docking atomic structures into low-resolution maps: application to cryo-electron microscopy. J Struct Biol. 2003, 141: 63-76.View ArticleGoogle Scholar
- Sherwood P, Brooks BR, Sansom MSP: Multiscale methods for macromolecular simulations. Curr Opin Struct Biol. 2008, 18 (5): 630-640.View ArticleGoogle Scholar
- Mereghetti P, Wade RC: Atomic detail Brownian dynamics simulations of concentrated protein solutions with a mean field treatment of hydrodynamic interactions. J Phys Chem B. 2012, 116: 8523-8533.View ArticleGoogle Scholar
- Sharp K: Incorporating solvent and ion screening into molecular dynamics using the finite-differnece Poisson-Boltzmann method. J Comp Chem. 1991, 12 (4): 454-468.View ArticleGoogle Scholar
- Lu BZ, Chen WZ, Wang CX, Xu Xj: Protein molecular dynamics with electrostatic force entirely determined by a single Poisson-Boltzmann calculation. Proteins: Struct Function Bioinform. 2002, 48 (3): 497-504.View ArticleGoogle Scholar
- Prabhu NV, Zhu P, Sharp KA: Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smooth-permittivity finite difference Poisson-Boltzmann method. J Comp Chem. 2004, 25 (16): 2049-2064.View ArticleGoogle Scholar
- Wang J, Tan C, Lu Q, Luo R, Tan Yh: Poisson-Boltzmann Solvents in Molecular Dynamics Simulations. Commun Comput Phys. 2008, 3 (5): 1010-1031.Google Scholar
- Ermak DL, McCammon JA: Brownian dynamics with hydrodynamic interactions. J Chem Phys. 1978, 69: 1352-1360.View ArticleGoogle Scholar
- Urbina-Villalba G, Garcia-Sucre M, Toro-Mendoza J: Average hydrodynamic correction for the Brownian dynamics calculation of flocculation rates in concentrated dispersions?. Phys Rev E. 2003, 68: 061408-View ArticleGoogle Scholar
- de la Torre JG, Huertas ML, Carrasco B: Calculation of hydrodynamic properties of globular proteins from their atomic-level structure. Biophys J. 2000, 78: 719-730.View ArticleGoogle Scholar
- Tokuyama M, Oppenheim I: Dynamics of hard-sphere suspensions. Phys Rev E. 1994, 50: R16-R19.View ArticleGoogle Scholar
- Tokuyama M, Oppenheim I: On the theory of concentrated hard-sphere suspensions. Physica A. 1995, 216: 85-119.View ArticleGoogle Scholar
- Tokuyama M, Moriki T, Kimura Y: Self-diffusion of biomolecules in solution. Physical Rev E. 2011, 83 (5): 1-8.View ArticleGoogle Scholar
- Cichocki B, Ekiel-Jezewska ML, Wajnryb E: Lubrication corrections for three-particle contribution to short-time self-diffusion coefficients in colloidal dispersions. J Chem Phys. 1999, 111 (7): 3265-3273.View ArticleGoogle Scholar
- 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. J Am Chem Soc. 2009, 131 (26): 9230-9238.View ArticleGoogle Scholar
- Gabdoulline RR, Wade RC: Effective charges for macromolecules in solvent. J Phys Chem. 1996, 100: 3868-3878.View ArticleGoogle Scholar
- Elcock A, Gabdoulline R, Wade R, McCammon J: Computer simulation of protein-protein association kinetics: Acetylcholinesterase- Fasciculin. J Mol Biol. 1999, 291: 149-162.View ArticleGoogle Scholar
- Neal BL, Asthagiri D, Lenhoff AM: Molecular origins of osmotic second virial coefficients of proteins. Biophys J. 1998, 75: 2469-2477.View ArticleGoogle Scholar
- Elcock AH, McCammon JA: Calculation of weak protein-protein interactions: the pH dependence of the second Virial coefficient. Biophys J. 2001, 80: 613-625.View ArticleGoogle Scholar
- Hill TL: An introduction to statistical thermodynamics. 1986, New York: Dover Publications IncGoogle Scholar
- da Silva MA, Itri R, Arêas EPG: Lysozyme viscoelastic matrices in tetramethylurea/water media: a small angle X-ray scattering study. Biophys Chem. 2002, 99 (2): 169-179.View ArticleGoogle Scholar
- 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. J Phys Chem B. 2007, 111: 251-9.View ArticleGoogle Scholar
- Heinen M, Zanini F, Roosen-Runge 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. Soft Matter. 2012, 8 (5): 1404-View ArticleGoogle Scholar
- Allen MP, Tildesley DJ: Computer Simulations of Liquids. 1989, USA: Oxford University PressGoogle Scholar
- Asthagiri D, Paliwal A, Abras D, Lenhoff AM, Paulaitis ME: A consistent experimental and modeling approach to light-scattering studies of protein-protein interactions in solution. Biophysical J. 2005, 88: 3300-3309.View ArticleGoogle Scholar
- Price WS, Tsuchiya F, Arata Y: Lysozyme aggregation and solution properties studied using PGSE NMR diffusion measurements. J Am Chem Soc. 1999, 121: 11503-11512.View ArticleGoogle Scholar
- Pieper U, Eswar N, Webb BM, Eramian D, Kelly L, Barkan DT, Carter H, Mankoo P, Karchin R, Marti-Renom MA, Davis FP, Sali A: modbase, a database of annotated comparative protein structure models and associated resources. Nucleic Acids Res. 2009, 37 (suppl 1): D347-D354.View ArticleGoogle Scholar
- Balbo J, Mereghetti P, Herten DP, Wade R: The shape of protein crowders is a major determinant of protein diffusion. Biophys J. 2013, 104 (7): 1576-1584.View ArticleGoogle Scholar
- Gordon JC, Myers JB, Folta T, Shoja V, Heath LS, Onufriev A: H++: a server for estimating pKas and adding missing hydrogens to macromolecules. Nucleic Acids Res. 2005, 33: W368-W371.View ArticleGoogle Scholar
- Jorgensen WL, Tirado-Rives J: The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J Am Chem Soc. 1988, 110 (6): 1657-1666.View ArticleGoogle Scholar
- 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. Comput Phys Commun. 1995, 91 (1-3): 57-95.View ArticleGoogle Scholar
- Balescu RC: Equilibrium and Non-Equilibrium Statistical Mechanics. 1975, New York: John Wiley & SonsGoogle Scholar
- Nägele G: On the dynamics and structure of charge-stabilized suspensions. Phys Rep. 1996, 272: 215-375.View ArticleGoogle Scholar
- Fukasawa T, Sato T: Versatile application of indirect Fourier transformation to structure factor analysis: from X-ray diffraction of molecular liquids to small angle scattering of protein solutions. Phys Chem Chem Phys. 2011, 13 (8): 3187-3196.View ArticleGoogle Scholar
- 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. Biophys J. 2010, 99 (5): 1342-1349.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. 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 credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.