INTRODUCTION

TOP ABSTRACT INTRODUCTION RATE CONSTANT CALCULATION WITH... BIASED BROWNIAN DYNAMICS TEST CASES DISCUSSION REFERENCES

Many biochemical reactions are classified as diffusion-limited, which means the encounter of the reactants is the most time-limiting part of the reaction. A number of reactions, such as enzyme-substrate and antibody-antigen association, belong to this category. Understanding the nature of the encounter in molecular detail will provide insights into underlying mechanisms of association. Brownian dynamics (BD) simulations have been widely used to calculate the rates of the diffusion-limited reactions. BD theory was developed by Einstein in 1905, Smoluchowski in 1907, 1913, 1915, Fokker in 1913, 1914, and Planck in 1917. The early work in BD is reviewed in an excellent paper by Chandrasekhar (1943). Kramers in 1940 first used BD to study the diffusion-limited reactions. There is a good review of reaction-rate theory by Hänggi et al. (1990). Computer simulations of diffusional motion in biomolecules were first developed by Ermak and McCammon (1978) and Northrup et al. (1984). A state-of-the-art BD simulation code for biomolecules, UHBD, which calculates biomolecular rates of association, is provided by McCammon and co-workers (Davis et al., 1991).

BD simulation is a CPU-intensive computation. The accuracy of the result is typically proportional to the inverse square root of the CPU time. To reduce the simulation time, especially for those reactions with high energy or entropy barriers, an enhanced sampling methodweighted ensemble Brownian dynamics (WEBD)was developed by Huber and Kim (1996). WEBD maintains an ensemble of particles (in configuration space), and associates a probability weight with each particle. All particles carry out Brownian motion without influencing each other. Periodically, WEBD splits and merges particles according to their position and weight. This procedure enables better sampling in regions inaccessible to standard BD; thus particles are sampled more near to reaction regions and the accuracy of the result is improved.

Here we propose a novel sampling methodbiased Brownian dynamics (bias-BD)for efficient computation of reaction rates. Bias-BD also associates a weight with each particle and changes the weight as the particle moves. Bias-BD introduces a biasing force in addition to the electrostatic force between the reactants. A particle loses weight when moving along the biasing force and gains weight when moving against the biasing force. As a consequence, particles are sampled more near to the reaction site, albeit with less weight. Bias-BD samples more of the important trajectories and reduces the variance of the result. To obtain the same computational precision in our test case, bias-BD needs only one-seventh the CPU time of the unbiased BD. By selecting a suitable biasing force one can control the sampling density in different regions. Bias-BD simulates only one particle at a time. There is no coupling between particles, as in WEBD, and hence bias-BD is embarrassingly parallel.

Variance-reduced methods for stochastic simulations are discussed in detail in Kloeden and Platen (1992, pages 511-528), which are based on Milstein (1988, page 225) and Wagner (1987, 1988). Milstein (1988) introduces an additional drift term in the stochastic differential equation, which is the same as doing a Girsanov transformation to the underlying probability measure. Wagner (1987, 1988) derives unbiased estimators for functional integrals of stochastic processes based on the general principles of Monte Carlo integration. These methods are used by researchers in other areas, such as Melchior and Öttinger (1995, 1996), and Zuckerman and Woolf (1999). Bias-BD is similar to these variance-reduced methods, but differs from them in three aspects. First, bias-BD is deduced from the discretized stochastic differential equation. The weight calculation is based on the random number distributions at each step and bias-BD is valid for each discretized step. Second, bias-BD does not require a unified termination time. In rate constant calculation, each trajectory has a different lifetime. The use of variance-reduced methods in Milstein (1988), Wagner (1987, 1988), and Kloeden and Platen (1992) is not straightforward here. Third, our choice of bias force is different.

In the following sections, we first give an overview of the Northrup, Allison, and McCammon (NAM) method (Northrup et al., 1984); however, our overview is based on a paper by Zhou (1990). Then we describe bias-BD and its weight calculation. Finally, we discuss two test cases and show our choice of bias force. The first is a 3D pure diffusion model, which we can compare with the analytical result. The other is the reaction rate between O₂ ions and super-oxide-dismutase (SOD). In the test cases we implement bias-BD in UHBD (Davis et al., 1991) and show a sevenfold improvement in CPU time.

	RATE CONSTANT CALCULATION WITH BROWNIAN DYNAMICS

TOP ABSTRACT INTRODUCTION RATE CONSTANT CALCULATION WITH... BIASED BROWNIAN DYNAMICS TEST CASES DISCUSSION REFERENCES

Rate constant

We consider the reaction between two types of molecules E (enzyme) and S (substrate). The molecules diffuse in solvent and have certain concentrations. A reaction happens when an S molecule is close enough to an E molecule's binding site. We will compute the reaction rate in a unit volume.

We assume the concentrations do not decrease because of reactions. We also assume concentrations are dilute, i.e., an S molecule's movement is affected by at most one E molecule, and is not affected by others of type S. Under these assumptions, we may consider the relative movement between one pair of molecules. Let k be the reaction rate for a simplified system, which has one E molecule and S moving relative to the E molecule with a unit relative concentration at infinity. If the actual concentration of S at infinity is p_S, the reaction rate will be kp_S. If we assume E's concentration to be p_E, or there are p_E E molecules in a unit volume, the reaction rate in a unit volume is kp_Sp_E, because each E molecule provides the same reaction rate independently. In the following discussion we will regard the pair E, S as a "particle" in configuration space.

The movement of S is described by a stochastic differential equationa diffusion limit Langevin equation, or Langevin equation under strong frictionwhich means the particle has no inertia, or the particle loses its momentum in a very short time (Ermak and McCammon, 1978):

(1)

In Eq. 1, R(t) represents the translational, rotational, and internal coordinates of particle S relative to E. D(r) is the relative diffusion coefficient between E and S as a function of coordinates in the configuration space. F(r) = U(r) is the electrostatic force that E exerts on S. The force vanishes at infinity, and we may let U(r) = 0 at infinity. W(t) is a vector of independent standard Wiener processes (Kloeden and Platen, 1992, page 70). In the following discussion we treat r as a 3D vector; however, the results extend to higher dimensions.

The system is equivalently described by the relative concentration p(r). In steady state, p(r) satisfies the steady-state Smoluchowski equation (Chandrasekhar, 1943; Northrup et al., 1984; Zhou, 1990):

(2)

where

(3)

As mentioned before, we assume p() = 1. Let the boundary of molecule E be _E. _E is partitioned into two parts: the binding site _E1 and the remainder _E2. The surface _E1 is an absorption boundary with p(r) = 0. The surface _E2 is a reflection boundary with n · j(r) = 0, where n is the outward normal direction on _E2.

The rate constant k is given by

(4)

where  is any surface that encloses _E (Zhou, 1990).

Domain truncation

The infinite domain rate constant k is hard to compute directly. If p(r) and U(r) are spherically symmetric for r  q with negligible error where r is the distance between E and S, the Smoluchowski Eq. 2 reduces to an ordinary differential equation for r  q that can be solved by 1D quadrature. The infinite domain problem is then truncated to a finite domain problem. Let

(5)

We call p_f(r) the finite domain concentration. It satisfies the same equation as p(r) and similar boundary conditions:

(6)

Let the rate constant for p_f(r) be k_f. From Eqs. 3-5 we have

(7)

In Eq. 4, let  be a spherical surface with radius r  q. Due to spherical symmetry we have

(8)

So,

(9)

Integrating Eq. 9 from q to +, we get

(10)

Dividing Eq. 10 by k on both sides and substituting in Eq. 7, we have

(11)

In Eq. 11, we see that to compute the infinite domain rate constant k, we can compute the finite domain rate constant k_f and a 1D quadrature. The finite domain problem has boundary condition p_f(q) = e^U(q)/k_BT. Eq. 11 can be explained in the following manner: the term 1/k is similar to the resistance in an electrical system. Two resistances are serially connected. One resistance is from the q surface to infinity whose value equals the 1D quadrature, and the other resistance is from the q surface to the reaction site whose value is 1/k_f. The term similar to the voltage in an electrical system is e^U(r)/k_BT p(r). To measure the resistance, we apply the voltage 1 at one end and voltage 0 at the other end, and measure the current. Eq. 11 tells us that the resistance of two serially connected resistances is the sum of the two.

FOP and NAM methods

To get the finite domain rate constant k_f, one approach is to solve the Smoluchowski equation directly. This can be extremely difficult, since it can be a high-dimensional partial differential equation (PDE) if we consider the rotational and internal degrees of freedom of the molecule. Another way is to use random walk, or BD, simulations. Generally speaking, the random walk method provides a simple way to solve high-dimensional PDEs.

One approach to BD simulations is the flux-over-population (FOP) method. In FOP, we set up a time period  and a thin shell with radii from q to q + q, and compute the expected number of particles in the shell n = 4q²qe^U(q)/k_BT. At the start of time period , generate n particles in the shell randomly. During , let all particles, those in the shell and those inside the q surface, do Brownian dynamics. At the end of , count the reacted particles, divide by  to get the average reaction rate for this  time period, and destroy all the particles beyond the q surface. The system reaches steady state after a number of time periods. After that, we calculate a rate for each time period and their mean value.

Alternatively, we can calculate a reaction probability P(q, ) by noting that particles in FOP method do not affect each other and each particle contributes to the reaction rate independently. A particle is generated randomly in the shell. It is allowed to follow a Brownian trajectory where its position is checked every  time. P(q, ) is the probability that the particle terminates by reacting instead of escaping out of the q surface. Then,

(12)

There is a restriction on q and  in FOP. Note that a particle outside the q + q surface has a finite probability of entering the q surface in  time. We require the probability to be small so that it can be neglected. For example, if q > 4 · (2D)^1/2, the probability is very small.

There is another method to compute k_fthe NAM methodproposed by Northrup, Allison, and McCammon in 1984 and further explained by Zhou in 1990.

The NAM method further assumes that U(r) is spherically symmetric for r  b for some b < q. Let (r) be the probability that a particle starting at r reacts instead of escaping out of the q surface, and let (r) = 1  (r). (r) satisfies the steady-state backward Smoluchowski equation and the given boundary conditions (Zhou, 1990):

(13)

Comparing Eq. 6 and Eq. 13, we see that (r) and e^U(r)/k_BT p_f(r) satisfy the same PDE and same boundary condition; so

(14)

We have

(15)

Taking  to be a sphere with radius r  b and noting that U(r) is spherically symmetric, we have

(16)

Let

(17)

then

(18)

Integrating Eq. 18 from b to q,

(19)

or

(20)

where

(21)

Thus, the NAM method connects the finite domain rate constant with a probability that could be calculated by BD simulations. Eq. 20 can be explained in the following manner: two resistances are connected in parallel. One resistance is from the b surface to the q surface whose value equals the 1D quadrature, and the other resistance is from the b surface to the reaction site whose value is unknown. Now we add some voltage at the b surface, and voltage 0 at both the q surface and the reaction site. With the value (b) that is the percentage of the current going through the unknown resistance, we can calculate the resistance when two resistances are serially connected. In the following discussion we will use  to represent (b): the average reaction probability starting from a randomly chosen point on the b surface.

Error estimation

If the 1D integrals in Eq. 11 and Eq. 20 are evaluated with negligible error, the only error source is from the computation of P(q, ) or the  value.

To calculate , let us start a particle S from b surface, simulating its movement until either it reacts or escapes. We then assign X_i = 1 or 0 for this trajectory, depending on whether it reacts or not. We simulate many trajectories and obtain many X_i values. Then we compute their mean as an approximation of ,

(22)

and estimate X's variance:

(23)

Let be the 90% confidence interval of the estimated . We have

(24)

The relative errors for k_f and k are

(25)

(26)

Typically,  is very small due to the energy and entropy barriers between the molecules. The relative error / is large. In order to get a suitable relative error, one has to sample many trajectories.

	BIASED BROWNIAN DYNAMICS

TOP ABSTRACT INTRODUCTION RATE CONSTANT CALCULATION WITH... BIASED BROWNIAN DYNAMICS TEST CASES DISCUSSION REFERENCES

As mentioned, bias-BD associates a weight with a particle, and changes the weight as the particle moves. To calculate , we assign X_i = w_i or 0 for each trajectory, depending on whether it reacts or not, where w_i is the final weight of the particle. We then calculate  and with Eqs. 22-24, as usual. The effect of bias-BD is that we have many reacted trajectories with small weights instead of a few reacted trajectories with unit weight. Most of the particles that escape have a weight larger than 1. The average weight of all trajectories is ~1. Bias-BD samples more near the binding site and less away from binding site. Thus it reduces the variance of the estimate.

Before examining bias-BD, we examine how standard BD does sampling. By discretizing the diffusion limit Langevin Eq. 1 with the Euler-Maruyama method, introduced into biophysics by Ermak and McCammon (1978), we get

(27)

where each Zⁿ is a vector of independent standard Gaussian distributed random numbers that comes from the discretization of the Wiener term W(t). Zⁿ has the probability density function (PDF)

(28)

Bias-BD selects random vector Z'ⁿ with a different PDF p'(z) and associates with each vector a weight function w(z), which maintains the following equality:

(29)

Thus, the weight distribution product preserves Gaussian distribution. The w(z) here is the weight for one step of integration. When a particle starts from b surface, its weight is 1. When a particle moves a step, the new weight equals the old weight times w(Z'ⁿ). The final weight is a product w(Z'¹)w(Z'²) ··· w(Z'^N).

The new PDF p'(z) should make reaction more likely. To get a systematic way to derive p'(z), let us introduce a biasing force F_b(r). F_b(r) acts on the particle like the electrostatic force F(r). The one step integral changes to

(30)

or,

(31)

where Z'ⁿ = Zⁿ + z₀, z₀ = F_b(Rⁿ)(Dt/2)^1/2/(k_BT). We see that Eq. 31 is the same as Eq. 27, except it uses random vector Z'ⁿ instead of Zⁿ. We know Zⁿ is Gaussian-distributed with PDF p(z). Z'ⁿ has PDF p'(z) = p(z  z₀). So, the weight function is w(z) = p(z)/p(z  z₀).

The final part is the determination of the biasing force. The biasing force affects the sampling density significantly. It acts like an additional potential. For instance, we can select this potential artificially to cancel the electrostatic force completely, so we will have equal sampling in all the configuration space. We can add a centripetal biasing force, so we get more sampling near the binding site. There are a multitude of choices.

	TEST CASES

TOP ABSTRACT INTRODUCTION RATE CONSTANT CALCULATION WITH... BIASED BROWNIAN DYNAMICS TEST CASES DISCUSSION REFERENCES

3D pure diffusion

We calculate the finite domain rate constant k_f of a 3D pure diffusion problem as described by the SDE Eq. 1 or by the PDE Eq. 6 with U(r) = 0 everywhere, k_BT = 1, D = 1, q = 15, and reaction occurring at r = a = 1. The analytical value of the rate constant is k_f = 4D/(1/a  1/q) = 13.464.

View larger version (24K): [in this window] [in a new window]   FIGURE 1   Measurement in million integration steps needed to achieve kf/kf < 1% in 3D pure diffusion test.

The 3D pure diffusion can be transformed to a 1D diffusion with an external force. Due to spherical symmetry of p_f(r), Eq. 6 can be transformed to

(32)

Letting p_1D(r) = 4r²p_f(r) and U_1D(r) = 2k_BT ln r, we have

(33)

which is the 1D Smoluchowski equation with external force (d/dr)U_1D(r) = 2k_BT/r.

If the bias force cancels the external force in 1D diffusion, we expect to obtain better results. These experiments use the centripetal bias force F_b(r) = bias × 2k_BTr/r², where bias ranges from 0 to 2. When bias = 0, it is standard BD; when bias = 1, the biased system is equivalent to a 1D pure diffusion. We try both FOP and NAM methods with the different strength of bias force. The NAM method is used with various choices for b. Numerical integration uses variable time steps. In the NAM method, time steps become smaller as r approaches a or q. In the FOP method, time steps become smaller as r approaches a and are fixed to the time period  as r approaches q. We calculate k_f with Eq. 12 for the FOP method and Eq. 20 for the NAM method. To measure computational cost, an estimate was calculated of the number of trials required to achieve a relative error of 1% with 90% confidence in the estimated value of k_f. This is multiplied by the average number of time steps per trial to yield one cost measure, as in Fig. 1, and by seconds of CPU time per trial to give a second cost measure, as in Fig. 2. Table 1 shows the best choice of bias in each configuration and the speedup relative to the standard unbiased BD when measured in number of time steps. Fig. 3 plots the number of reacted trajectories for different choices of bias. Fig. 4 compares the computed k_f and k_f values with the analytical k_f value.

View larger version (26K): [in this window] [in a new window]   FIGURE 2   Measurement in CPU seconds needed to achieve kf/kf < 1% in 3D pure diffusion test.

View larger version (24K): [in this window] [in a new window]   FIGURE 3   Number of reacted trajectories growing with bias in 3D diffusion test.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Comparison of simulated kf value and analytical value in 3D diffusion test.

We see that the FOP method has the same efficiency as the NAM method when the b surface is very close to the q surface. The closer the b surface is to the reaction site, the more efficient is the NAM method. But for a general problem, another error becomes significant for the small b surface: the asymmetry of the potential U(r). Hence, there is a tradeoff in the selection of b. The best bias value is between 1.2 and 1.5. This implies that the biasing force should be somewhat more centripetal than what is needed to get pure 1D diffusion.

Rates between O₂ ions and SOD

In this experiment we calculate the reaction rate between O₂ ions and SOD. Data of SOD are taken from the Protein Data Bank and the partial charge table is taken from the standard CHARMM force field (Brooks et al., 1983) incorporated into UHBD. SOD is represented by a set of point charges in ionized solvent and fixed in space. The electrostatic potential surrounding SOD is computed by solving the Poisson-Boltzmann equation on a 3D grid. O₂ is treated as a point charge and does BD movements relative to SOD under the influence of the electrostatic field. Reaction is said to happen when O₂ is within 7 Å of the copper atom of SOD. O₂ and each atom of SOD has a hard sphere radius, which prevents O₂ from penetrating into the SOD. When O₂ has a hard sphere contact with the nonreacting region of SOD, one simply retries the integration step with another Gaussian-distributed random vector.

We use F_b(r) = 2k_BT(r  r₀)/|r  r₀|² as the biasing force, where r₀ is the position of the copper atom of SOD. The effect of this force is to have unbiased sampling in the radial direction if there is no electrostatic force. The biasing force and the weight calculation are implemented in UHBD (Davis et al., 1991). Numerical integration uses variable time steps that become smaller as reactants are closer.

Table 2 shows the results. Though the average cost per trajectory of the bias-BD is a little higher than that of the standard BD, the 90% confidence interval of the bias-BD is much less. The last line in the table shows that the overall cost of the bias-BD to get 1% relative error is about 1/7 of that of the standard BD.

Fig. 5 compares fluctuations of the  value of bias-BD and standard BD. The solid line is that of bias-BD, and the dashed line is that of standard BD. Each point in the figure represents an estimated  value from 100 trajectories. The horizontal axis shows there are 100 independent runs. We note that bias-BD yields much less fluctuation of the  value, i.e., reduces variance.

View larger version (31K): [in this window] [in a new window]   FIGURE 5   Fluctuation of  in bias-BD and standard BD in SOD test.

	DISCUSSION

TOP ABSTRACT INTRODUCTION RATE CONSTANT CALCULATION WITH... BIASED BROWNIAN DYNAMICS TEST CASES DISCUSSION REFERENCES

In previous sections we reviewed the theory for reaction rate computation and proposed bias-BD to do enhanced BD sampling. We see that the reaction system is similar to an electric system, where the conductance is the rate constant, the voltage is e^U(r)/k_BT p(r), and the current density is the flux. Bias-BD greatly reduces variance in BD simulations of reactions with high energy or entropy barriers if we select the random variable X_i to be w_i or 0, depending on whether it reacts or not, where w_i is the final weight of the particle.

The main limitation is the ignorance of the choice of bias force that would be optimal for minimal variance. Our choice that the bias force renders the 3D diffusion into a 1D diffusion problem without an external force while being a good choice is perhaps not the best. Our experiments suggest that the bias force could be more centripetal than the force needed to make the 3D diffusion into a 1D diffusion problem. Further study is in progress to find better choices of biasing forces.

The advantage of BD occurs primarily for high-dimensional configuration spaces. For such problems the entropy barrier will be higher, and the reaction probability will be less. We expect that biased BD can give better improvement. It remains a challenge to find suitable biasing forces in high dimensions that can provide steering automatically.