INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The kinetics of bimolecular associations play an important role in biological processes. An improved understanding of how molecules use electrostatic and hydrodynamic steering to speed up their association rate could lead to the design of genetically engineered enzymes with a high turnover rate. Enzyme-substrate bimolecular reaction simulations using the Northrup-Allison-McCammon method have been applied to numerous bimolecular encounters (Antosiewicz et al., 1995; Antosiewicz and McCammon, 1995; Davis et al., 1991b; Gabdouline and Wade, 1997; Kozack et al., 1995; Kozack and Subramaniam, 1993; Luty et al., 1993a,b; Wade et al., 1994). For reaction systems with large free energy barriers, enzyme-substrate encounter events occur infrequently, thus warranting extensive sampling of the phase space. Consequently, standard University of Houston Brownian Dynamics (UHBD) simulations (Davis et al., 1991a) require many Brownian dynamics trajectories for reaction rates to be predicted with high confidence. In some cases, the number of trajectories required may become so prohibitively large that the simulation becomes impractical.

Huber and Kim presented the idea of Weighted Ensemble Brownian (WEB) dynamics as an alternative method for simulations of enzyme-substrate reactions in which reaction events occur infrequently (Huber and Kim, 1996). This algorithm is based on the flux overpopulation method, in which reaction rates are determined from the reactive flux (molecule reacted per unit time) in steady-state simulations. This refers to a type of simulation where multiple Brownian dynamics trajectories are simulated at once and the trajectories are immediately reinitialized once their destination is reached. The main difficulty in determining the reactive rates by either simulation method (traditional BD or the one based on the flux overpopulation method) is that most interesting systems require particles to overcome large free energy barriers before reaching their destination. As a result, the particles spend most of their time wandering within local free energy wells and only rarely surmount the barriers. To obtain good estimates of the reactive rates with traditional BD, prohibitive amounts of computer time are necessary. However, as we have shown before (Rojnuckarin et al., 1998), the flux overpopulation method can be modified to obtain meaningful results in a much more efficient manner.

Each particle (diffusing substrate) in a WEB simulation is regarded as a collection of probability packets; i.e., a single particle is split and combined, depending on its grid location, into weighted probability packets. The configuration space available to the particles is divided into regions (or bins) along the reaction coordinate. Regions that represent high free energy barriers may never be sampled in a standard BD simulation. Whereas, in the WEB method, because a single particle is split into pseudoparticles (probability packets) with smaller weights, there is a likelihood of sampling more regions of configuration space, even if only with smaller weights. To prevent the number of pseudoparticles from exploding, we also combine these pseudoparticles that arrive in the same bin in configuration space during the simulation. The net probability is conserved and the weights associated with the pseudoparticles satisfying the reaction condition reflect the probability of reaction. A schematic diagram of the WEB method is presented in Fig. 1. The main advantage of this method is the opportunity to use a diffusing substrate to sample regions of configuration space that represent high free energy barriers.

View larger version (43K): [in this window] [in a new window]   FIGURE 1   (a) Schematic diagram of the NAM method. (b) Each configuration in the simulation is assigned a weight that a particular realization of the system is in that configuration. The pseudoparticles are split and combined to get better sampling in the interesting parts of configuration space.

The probability distribution at the top of a large free energy barrier is poorly sampled in a distribution of equal weights; however, if the particles have unequal weights, the probabilities of the particles can be adjusted such that the region is adequately sampled. The bias results in more particles crossing free energy barriers and reaction events occurring more frequently. The probabilistic weight adjustment through splitting and combining rigorously corrects for the bias such that the weighted-average result from WEB simulation should be identical to the standard BD result; although the reaction events occur more frequently in WEB simulations, each event carries only a minuscule probabilistic weight. In addition, because reaction events occur more frequently in WEB simulations, the reaction rate can be determined with higher confidence in less time than it can be determined with standard BD. We can also trade the WEB algorithmic speed-ups for more realistic interaction models that may otherwise lengthen the simulation time. Complete details of the WEB dynamics algorithm have been published elsewhere (Huber and Kim, 1996; Rojnuckarin et al., 1998) (the WEB/UHBD package; the WEB part of the package is available from the NCSA Computational Biology Group upon request); only the basic notions will be presented here.

The WEB algorithm is incorporated into UHBD to utilize the sophisticated electrostatic routines built into UHBD (Davis et al., 1991a; Holst et al., 1994) to calculate forces between enzymes and substrates. Furthermore, incorporation of WEB into UHBD allows the UHBD-user community to take advantage of this efficient simulation algorithm. To outline the capabilities and limitations of the WEB algorithm, we describe here the application of the WEB method to three model systems. First we attempt to reproduce the Smoluchowski analytical reaction rate for the association of two spheres. We also test the algorithm with two protein-substrate systems: the association reaction of superoxide (O₂) with the enzyme copper-zinc superoxide dismutase (CuZn SOD) and the binding of the monoclonal antibody NC6.8 fragment with the anti-sweet-taste ligand N-(p-cyanophenyl)-N'-(diphenylmethyl)-guanidinium acetic acid.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

UHBD

The UHBD suite of programs (Davis et al., 1991a; Madura et al., 1995) can be used to simulate diffusion-limited protein-substrate reactions. In the UHBD framework, the protein is modeled at the atomic level of detail, while the substrate is modeled as a collection of spherical subunits, each of which may represent an atom or a group of atoms in the substrate. The enzyme and substrate interact via electrostatic interaction and volume exclusion. The electrostatic force calculation is based on the so-called point charge approximation, where point charges, q, placed on the subunits interact with the electric field, E, around the enzyme by a simple relation, F = qE. The electric field around the enzyme corresponds to the solution of the Poisson-Boltzmann equation (PBE), which is solved only once in the absence of the substrate. Charge distributions on the enzyme and protein low-dielectric interior are taken into account in the numerical solution of the PBE. The enzyme-substrate relative diffusion constant is calculated with the Rotne-Prager-Yamagawa correction for the hydrodynamics interaction between subunits (García de la Torre and Bloomfield, 1981). The Ermak-McCammon Brownian dynamics algorithm (Ermak and McCammon, 1978) and SHAKE-HI algorithm (Allison and McCammon, 1984) are used, respectively, to generate the configuration at the next time step and to preserve geometric constraints between subunits.

UHBD calculates the diffusion-limited bimolecular reaction rate via the Northrup-Allison-McCammon (NAM) method (Northrup et al., 1984). This method involves constructing a sphere of radius b large enough such that the interaction between substrate and protein depends only on the distance from the protein, and a larger sphere of radius q around the enzyme (Fig. 2). The probability, , that a Brownian trajectory starting on the b surface arrives at the active site before leaving the termination sphere of radius q is estimated by the UHBD Brownian dynamics simulation described previously. The reaction rate, k, is related to the probability  by the following relation:

(1)

In the above expression, D denotes the relative diffusion constant of enzyme and substrate, and U denotes the interaction potential between the protein and the substrate.

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Schematic diagram of the NAM method. Particle A reacts at the active site, while particle B proceeds beyond the termination radius.

Because WEB dynamics simulation is a steady-state simulation based on the flux-overpopulation method, we can estimate  from the ratio of steady-state reactive flux, f_SS^React, and the steady-state flux across the q surface, f_SS^QSurf, provided that we always reinitialize the reaction system on the b surface:

(2)

Substituting Eq. 2 into Eq. 1 gives the expression for the reaction rate used in this work.

Reaction coordinates adaptation

For the WEB dynamics algorithm to bias the enzyme and substrate toward a reacting configuration, it divides the configuration space into regions or bins along the reaction coordinate parameter, a measure that indicates the progress toward the reaction. The WEB algorithm splits and combines the probabilistic weight, which is assigned to each Brownian dynamics particle in such a way that the configuration space in different bins is sampled equally. This bias improves the sampling in the regions near the active site that are not easily accessible with standard BD because of the free energy barriers. Consequently, defining the reaction coordinate is the key step of the WEB dynamics algorithm.

For a smooth transition from standard UHBD to WEB/UHBD, we automatically define the reaction coordinate based on the definition of the reaction criteria from the standard BD simulation. In UHBD, the reaction criterion is defined for each reaction site by a set of distances between subunits of the substrates and atoms of the enzymes. If a substrate configuration satisfies all of the criteria for a reaction site, it is considered to have formed the association complex. If we define 1) the violation (in Å) of a reaction criterion as the distance between the subunit/atom pair specified by the criterion less the required distance, and 2) the maximum violation (in Å) of a reaction site as the maximum violation of all reaction criterion for that particular reaction site, we can use the minimum of the maximum violation of all reaction sites as our reaction coordinate. For example, the monoclonal antibody NC6.8 has a single binding site, and its reaction criteria are that the positive subunit of the ligand must be within 5 Å of Glu¹⁶⁹ OE2, and its negative subunit within 5 Å of Arg¹⁶² NH1 (See Simulation Setup below). If we have a ligand configuration whose positive subunit is 4 Å from Glu¹⁶⁹ OE2 and whose negative subunit is 7 Å from Arg¹⁶² NH1, the violations for the two reaction criterion would be 1 Å and 2 Å, respectively. As a result, the reaction coordinate in this case would be 2 Å.

It has been pointed out that any measure that corresponds to the progress of the reaction can be used as the reaction coordinate in WEB dynamics simulations, and the definition of the reaction coordinate in any WEB dynamics system is often left to the implementer (Rojnuckarin et al., 1998). The reaction coordinate defined above can certainly be used as a measure of progress of a reaction. A large violation implies that the substrate is far from the active site, a small violation implies that the substrate is close to reacting, and a negative maximum violation means that the substrate has reacted. The advantage of the above definition is that because the reaction criteria are already defined in the standard UHBD simulation, users do not need to supply any extra parameters in the definition of the WEB dynamics reaction coordinate. In addition, this definition is computationally efficient, as the routine that decides whether the reaction occurs and the routine that calculates reaction coordinates can be easily combined.

We build the bins along reaction coordinates by the procedure described by Huber and Kim (1996). In this procedure, 100 substrate molecules are first initialized on the q surface. Next substrate molecules are moved by one Brownian/SHAKE step and the reaction coordinates for these molecules are calculated and ranked; the reaction coordinate that corresponds to the first quintile is marked as a bin boundary. Substrates with reaction coordinates larger than the bin boundary (large violation) are removed and randomly inserted at the positions of molecules with reaction coordinates smaller than the bin boundary. The process is repeated until the bin boundary becomes negative. These bin boundaries are then rescaled to the desired number of bins specified by the users.

Program design consideration

While UHBD is designed to simulate only one BD trajectory at a time, the WEB dynamics algorithm requires that many BD trajectories be simulated concurrently so that the associated probabilistic weight can be split and combined. To be able to use as much of the original UHBD Brownian dynamics code as possible, we implement the WEB dynamics algorithm by storing the current enzyme-substrate configurations in a separate array. For each trajectory, we copy the current configuration into the appropriate FORTRAN common block variables that the original code uses and call the appropriate functions to generate the configuration at the next time step. Once we have stepped all trajectories forward by one time step, the routine that splits and combines the associated probabilistic weight is called, and the process starts over for the next time step. In addition, to save memory, the WEB-related routines reuse the memory allocated in the electrostatic calculations.

Simulation setup

Smoluchowski reaction rate

(3)

(4)

Superoxide dismutase

View larger version (68K): [in this window] [in a new window]   FIGURE 3   Structure of bovine Cu, Zn superoxide dismutase dimer (entry 2sod). The Cu actives sites are displayed as solid spheres.

Monoclonal antibody NC6.8

View larger version (77K): [in this window] [in a new window]   FIGURE 4   Structure of the Fv fragment of monoclonal antibody NC6.8 (entry 2cgr) complexed with its substrate.

View larger version (24K): [in this window] [in a new window]   FIGURE 5   The dumbbell (solid spheres) is overlaid on the ligand structure (ball and stick). The positive lobe is rendered in black, and the negative lobe is rendered in gray. Side chains of Glu H:50 and Arg H:57 are also shown.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Smoluchowski reaction rate

In the Smoluchowski setup we compare the analytical value of  (Eq. 4) to the standard BD and WEB results, using time steps of 0.2 ps. Solving Eq. 4 with the parameters given in the Methods section yields a value of 0.111 for . Simulation results of  values equal 0.111 ± 0.006 and 0.112 ± 0.004 for the standard BD and WEB runs, respectively. Using a variable time step with WEB improves the accuracy of the result. In this run, the time step used, t, is a function of separation distance r (in Å):

(5)

This simulation yields a  value of 0.111 ± 0.006, which agrees very well with the analytical result.

Superoxide dismutase

CuZn SOD catalyzes the reaction that converts the toxic superoxide free radical (O₂) to oxygen and hydrogen peroxide. The bovine enzyme is a homodimer consisting of 151 amino acid residues, with one copper and one zinc ion per subunit. Each monomer is composed of a slightly flattened Greek key -barrel core. Dismutation of superoxide occurs by the alternate reduction and oxidation of the active site copper ion, yielding one molecule each of oxygen and hydrogen peroxide. The enzyme is extremely efficient, with a rate constant close to that obtained in the diffusion limit. However, the active site constitutes a very small portion of the enzyme surface, resulting in a failure of a "uniform collisions" mechanism to accurately predict the enzyme's efficiency. The reaction rate has been shown experimentally (Lepock et al., 1985) to be dependent on the ionic strength of the solvent, suggesting that the association is steered by electrostatic guidance. Because of the electrostatic guidance mechanism, CuZn SOD has been extensively studied through the use of standard Brownian dynamics (Allison et al., 1988; Sines et al., 1990, 1992; Holst et al., 1994), making it an ideal test case for the WEB/UHBD program.

The predicted rate constants from the SOD simulations are reported in Table 1. The calculated rates of 5.92 × 10⁹ M¹ s¹ (standard BD) and 6.08 × 10⁹ M¹ s¹ (WEB) agree within the computed confidence intervals and are comparable to the reported experimental value of 3.92 × 10⁹ M¹ s¹ (Polticelli et al., 1994). From simulation, the bimolecular reaction rates for various CuZn SODs have been consistently calculated to be in the 10⁹ M¹ s¹ range. Sines et al. (1990) calculated the ionic strength dependence of bovine CuZn SOD and found the rates to occur in this range. Polticelli et al. (1994) showed by both experiment and BD simulations that the reaction rate of CuZn SODs from Bos taurs and Shark Prionace glaca both occur in this range. The reaction rate of human CuZn SOD has been calculated to be slightly less than our bovine results (2.2 × 10⁹ M¹ s¹) but still falls in the × 10⁹ M¹ s¹ range (Fisher et al., 1997). It should be noted that the standard error of the WEB result (for the CuZn SOD case) is significantly larger than that of the standard BD result. Because of strong electrostatic steering forces, CuZn SOD is known to be a hyperefficient enzyme without appreciable free energy barriers. The lack of free energy barriers, coupled with an efficient steering mechanism, thus makes the WEB dynamics approach redundant in this case. However, the proximity of values computed by the standard NAM and WEB methods to each other as well as to the experimental value validates the WEB method.

Monoclonal antibody NC6.8

Antibody-antigen binding is an important process in our immune response to foreign substances. The binding process can also be thought of as a enzyme-substrate reaction in which the antigen corresponds to the substrate, the antibody corresponds to the enzyme, and the binding event corresponds to the reaction event (Kozack et al., 1995). Because the antibody-antigen binding equilibrium constant is usually very high, the binding process can be treated as an irreversible reaction. In the case of the monoclonal antibody NC6.8 and its anti-sweet-taste ligand (N-(p-cyanophenyl)-N'-(diphenylmethyl)-guanidinium acetic acid) the electrostatic steering effect generated by the antibody is not as strong as that of the CuZn SOD, and the reaction rate is expected to be a few orders of magnitude smaller. Consequently, the simulation study of NC6.8 binding would likely benefit from the use of the WEB dynamics algorithm.

To understand the chemical basis of association between mAb NC6.8 and the anti-sweet-taste ligand, we have previously used continuum electrostatics calculations (Livesay et al., manuscript submitted for publication). Complex formation is mediated by a pair of salt bridges between the ligand and the antibody. The carboxy moiety of the ligand forms a salt bridge with the side-chain guanidinium ion of Arg H:57. The guanidinium ion moiety of the ligand forms a salt bridge with the carboxy side chain of Glu H:50.  stacking interactions also play a large role in the specificity of the association. The sweetener is made up of three phenyl rings, and there are no less than 14 aromatic residues within a 15-Å sphere about the ligand. Tyr L:101 is involved in extensive  stacking interactions with the p-cyanophenyl moiety of the sweetener. Tyr L:37 is in contact with one phenyl of the diphenylmethyl moiety, whereas Tyr H:100 contacts the other. Trp H:33 forms a van der Waals contact with both carbons of the acetic acid moiety of the ligand. In addition, two H-bonds between the ligand and protein add to the stability of the complex. Thus unlike in the SOD association with superoxide ion, the association of the sweetener with the antibody fragment is weakly steered and stereochemically involves a more intricate encounter.

When large energy barriers must be crossed for associations to occur, as in this case of NC6.8:ligand binding, WEB algorithm proves to be a much better alternative. The large energy barriers associated with NC6.8:ligand binding arise from two sources. First, compared with CuZn SOD:superoxide association, there is a decreased amount of electrostatic steering in the mAb NC6.8:ligand association. Second, the multiple reaction criteria ensure a decrease in the probability of achieving a successful reaction. Defining the reaction criteria in this way means that even though the substrate is in the proximity of the binding site, a reaction does not occur until the substrate orientates itself in a very particular way in relation to the antibody. Both result in the need for many standard BD trajectories to achieve good statistics. The considerable free energy barrier in this case results in an association constant approximately two orders of magnitude lower than that of CuZn SOD (exact experimental data for the mAb NC6.8:ligand example are unavailable). Table 2 shows that the WEB result has a significantly smaller confidence interval than the standard BD result. The smaller confidence interval results from the improved efficiency of WEB compared to standard BD. Therefore, simulation times needed for low 90% confidence intervals can be drastically reduced using WEB. If it is assumed that the width of the confidence interval is inversely proportional to the square root of the computation time, standard BD would require approximately eight times more CPU cycles than WEB to reach the same confidence interval.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Quantitative assessment of macromolecular association is of paramount importance in the study of biochemical phenomena. Both experimental stopped-flow techniques and computational approaches that involve calculating encounter rates from diffusion equations have been employed in the assessment of on-rates of macromolecular association (Hibbits et al., 1994; Northrup and Erickson, 1992). The computational approach is limited by two factors: first, the ability to sample configuration space efficiently, and second, incorporation of intermolecular forces, effects of the medium, and stereochemical factors. We show here that the WEB method is very efficient in sampling the configuration space. The WEB method yields analytically accurate results in the Smoluchowski model and results comparable to both experimental and other accurate theoretical methods in the case of association of the enzyme superoxide dismutase with the superoxide ion. In the case where substantial free energy barriers exist, such as in the association of the sweetener ligand with the monoclonal antibody fragment NC6.8, the WEB dynamics algorithm is much more efficient. Standard BD approaches fail to yield acceptable confidence levels in comparable computing time scales. The WEB method has been integrated into the UHBD suite of programs and can be used to simulate bimolecular encounters (Rojnuckarin et al., 1998).