Biophysical Journal 87:1478-1497 (2004)
© 2004 The Biophysical Society
Dependence of DNA Polymerase Replication Rate on External Forces: A Model Based on Molecular Dynamics Simulations
Ioan Andricioaei *,
Anita Goel * 
,
Dudley Herschbach * and
Martin Karplus * 
* Department of Chemistry and Chemical Biology, Harvard University, Cambridge, Massachusetts 02138 USA; Department of Physics, Harvard University, and Harvard MIT Joint Division of Health Sciences and Technology, Cambridge, Massachusetts 02138 USA; and Institut Le Bel, Université Louis Pasteur, Strasbourg, France
Correspondence: Address reprint requests to Martin Karplus, E-mail: marci{at}tammy.harvard.edu.
 |
ABSTRACT
|
|---|
Molecular dynamics simulations are presented for a Thermus aquaticus (Taq) DNA polymerase I complex (consisting of the protein, the primer-template DNA strands, and the incoming nucleotide) subjected to external forces. The results obtained with a force applied to the DNA template strand provide insights into the effect of the tension on the activity of the enzyme. At forces below 30 pN a local model based on the parameters determined from the simulations, including the restricted motion of the DNA bases at the active site, yields a replication rate dependence on force in agreement with experiment. Simulations above 40 pN reveal large conformational changes in the enzyme-bound DNA that may have a role in the force-induced exonucleolysis observed experimentally.
|
INTRODUCTION
|
|---|
DNA replication, a crucial process for the propagation of the genome, is catalyzed by DNA polymerases (Kornberg and Baker, 1991
). These motor enzymes function in repeated cycles during which they move along a partly single-stranded, partly double-stranded DNA chain (Joyce and Steitz, 1994). In each catalytic cycle, a single nucleotide, complementary to the nucleotide on the template strand (if no error is introduced), is added to the primer strand of the double-stranded DNA. Thus, the double-stranded DNA chain with N basepairs changes to one with N + 1 basepairs. Because the fidelity of the base addition is fundamental for genetics and elucidating its mechanism may aid in the design of therapeutic agents against rapidly mutating pathogens, an understanding of the incorporation step performed by the polymerases is important. It is essential for a full description of DNA replication.
Based on sequence similarity, polymerases can be classified in seven different families (Patel and Loeb, 2001). The most extensively studied, both kinetically and structurally, are those in family A (which includes the prokaryotic and archaea Pol I, the eukaryotic polymerases 
and 
, as well as the viral T3, T5, and T7 polymerases) and those in family B (which include the prokariotic Pol II, several eukariotic and archaeal polymerases, as well as the viral adenovirus, HSV, RB69, T4, and T6 polymerases). Despite their diverse biological functions (replication, recombination, repair) the structural and chemical mechanisms of base incorporation by the polymerases in these two families seem to be very similar (Patel and Loeb, 2001; Steitz, 1998; Brautigam and Steitz, 1998). For a number of polymerases in family A (T7 polymerase, Doublié et al., 1998; Taq polymerase I, Li et al., 1998b) and in family B (RB69 polymerase, Franklin et al., 2001) there exist crystal structures for ternary intermediate complexes (i.e., DNA primer/template and incoming nucleotide in complex with the intermediate (active) state of the enzyme). These ternary structures are important for understanding the mechanism of the incoming nucleotide incorporation. The gross features of these DNA polymerases have been likened to a right hand (Ollis et al., 1985), with the fingers interacting with the incoming nucleotide and the template, the palm containing the catalytic site and binding to the incoming nucleotide, and the thumb binding the double-stranded DNA. The incoming nucleotide forms H-bonds with its partner (template base) and is stacked onto the last primer base (see Fig. 1).
View larger version (102K):
[in this window]
[in a new window]
|
FIGURE 1 Atomic detail of the active site of the open (left) and closed (right) conformations of Taq polymerase (in blue surface representation) complexed with the DNA, depicted in licorice representation (primer in red, template in gray, incoming dCTP in yellow). Only protein regions within 7 Å of the DNA bases shown are represented. The hydrogen atoms on the sugar-phosphate backbone are omitted for clarity.
|
|
A common mechanism for these polymerases has been proposed (Patel and Loeb, 2001; Steitz, 1999), based on the available crystal structures. In crystal structures of the apo form (Kim et al., 1995; Nayal et al., 1995) (i.e., with neither DNA nor the nucleotide bound), in those of the binary complexes with dNTP (Ollis et al., 1985; Li et al., 1998a) or with DNA (Beese et al., 1993; Eom et al., 1996; Kiefer et al., 1998; Li et al., 1998b), and in that of a ternary complex (i.e., with DNA and incoming nucleotide; Li et al., 1998b) a relaxed, "open" configuration of the fingers domain has been observed. In other ternary complex structures (Doublié et al., 1998; Li et al., 1998b; Franklin et al., 2001), this domain is in a "closed" configuration. In going from the open structure to the closed structure, the fingers domain closes by 
40–60° and the template DNA base that is to pair with the incoming nucleotide rotates inward by >90° (see Fig. 2), placing the incoming nucleotide in an optimal alignment for its subsequent chemical incorporation into the DNA primer. Although the details of the transition are not entirely understood, this major conformational change takes place before chemistry can occur, and is believed to be the rate-limiting step in a kinetic mechanism (Patel et al., 1991; Dahlberg and Benkovic, 1991) for the overall reaction. Additionally, this mechanism is consistent with an early suggestion (Bryant et al., 1983) (i.e., one proposed before the availability of crystal structures) that dNTP first binds to a polymerase/primer/template complex independently of the template information, and that only after the rate-limiting step (believed to be the closure of the fingers) does tight binding become template dependent. Among the polymerases for which ternary closed structures exist (i.e., the intermediate (active) ternary structures mentioned above with the fingers in the closed configuration) the large fragment of Taq polymerase I (Li et al., 1998b) (Klentaq1) is particularly appropriate for a computational analysis because, in addition to the closed ternary complex structure with ddCTP bound, a ternary open structure with ddCTP bound has been reported (Li et al., 1998b). Both open and closed structure crystals diffracted to similar resolution (2.3 Å) in the same space group, and were grown using identical crystallization procedures, with the exception that the open form was obtained from a selenomethionine-substituted protein and, after crystallization, was incubated with a washing solution that partially depleted the population of bound nucleotide. This procedure is believed (Li et al., 1998b) to have caused the transformation of an originally closed crystal structure to the open form. Closed ternary structures have been determined also for ddATP-, ddTTP-, and ddGTP-trapped complexes, in addition to that with the ddCTP nucleotide (Li and Waksman, 2001), but no open ternary structure with the last three nucleotides is available.
View larger version (43K):
[in this window]
[in a new window]
|
FIGURE 2 "Minimalist" representation of the open and closed conformational change of Taq-DNA complex in the immediate vicinity of the active site; this can be compared with Fig. 1. The ds DNA template extending to the 3' side comprises nucleotides labeled +1, +2, ...; template nucleotides on the 5' side are labeled 0, –1, ...; the primer strand is in red. In the open complex, the incoming nucleotide (in yellow) encounters a Tyr (blue licorice) whereas its destined partner base (gray licorice, at 0) is tilted sideways and inaccessible. In the closed complex, the Tyr is displaced by the change in the position of the O helix (blue cylinder) and the partner base is flipped into a bonding configuration in the space previously occupied by the Tyr.
|
|
Recently, the structural and kinetic data for family A polymerases have been supplemented by single-molecule experiments, in which it was shown that the rate of the replication reaction catalyzed by the DNA polymerases is altered when a force is applied to the template strand (Wuite et al., 2000; Maier et al., 2000). T7 DNA polymerase was studied with an optical trap (Wuite et al., 2000) and the large (Klenow) fragment of Escherichia coli, its 3'–5'exonuclease deficient mutant (Klenow exo–), and a 3'–5'exonuclease deficient mutant of T7 DNA polymerase (Sequenase) were studied by use of a magnetic trap (Maier et al., 2000). The two sets of experiments showed similar behavior. It was observed that the replication rate decreases at high forces and appears to increase at low forces. The experimental uncertainties are such that, although the rate decrease at high applied force is unequivocal, the rate increase at low force is within the error bars. The qualitative demonstration that the rate of polymerization depends strongly on the applied force f is of interest per se; i.e., it shows that the limiting (rate-determining) step involves work by the polymerase complex (and therefore motion) against an external force. However, a quantitative interpretation of the experimental results requires introduction of a model expressing the replication rate in terms of the force-induced work.
It is possible, in principle, that the application of the force transforms a nonlimiting step (e.g., translocation) into a limiting step. Therefore, the fact that rate depends on force does not always mean that motion is involved in the limiting step. However, the opposite is true: if there is no rate-force dependence, then there is no motion in the rate-limiting step. For RNA polymerases (Wang et al., 1998) and DNA helicases (V. Croquette, ENS, Paris, personal communication), analogous single-molecule experiments were crucial in showing that the rate-limiting step is independent of an applied force.
The models proposed previously to interpret the single-molecule experiments are phenomenological in nature and assume Arrhenius behavior (Wang et al., 1998), as in other studies of the effect of external forces on molecular reactions, such as protein unfolding (Carrion-Vazquez et al., 1999) or ligand unbinding (Moy et al., 1994). The rate constant in the absence of the force, k(0), is related to that in the presence of the force, k(f), by the equation
 | (1) |
with 
g(f) = G(f) – G(0) the force-induced barrier change relative to the zero-force activation barrier, G(0). Equation 1 implicitly assumes that k(0) also has an Arrhenius-type temperature dependence and that the magnitude of the force is such that the system is not in the so-called diffusive limit (Izrailev et al., 1997; Evans and Ritchie, 1997). Using Eq. 1, two types of models have been proposed. Both assume the closed conformation of the enzyme complex can be taken as a surrogate for the transition state. One model (Wuite et al., 2000; Maier et al., 2000), which we refer to as "global," evaluates the force-induced activation barrier from experimentally determined extension versus force curves for "bare" single-stranded (ss) and double-stranded (ds) DNA polymers, thousands of basepairs long and not interacting with enzyme. The other model, referred to as "local" (Goel et al., 2001, 2002), focuses solely on conformational changes of just two DNA segments (each involving one nucleotide) in the neighborhood of the active site of the enzyme.
From the global model it was concluded that the best fit to the replication rate data was obtained if, in the rate-determining step, more than one (n = 2 or 4, depending on the enzyme) of the single-stranded nucleotides at the end of the duplex were converted from ss to ds geometry, only to have n – 1 of them (i.e., all but the one that becomes part of the ds DNA) return to ss geometry before the next catalytic cycle. If correct, this conclusion would have important mechanistic implications. The local model also obtained fair agreement with the rate data, but assumed that only one ss base (n = 1) of the two bases considered in the ss portion of the template is converted to ds geometry in the closed state of the enzyme complex. This is consistent with the ternary crystal structures (Doublié et al., 1998; Li et al., 1998b; Franklin et al., 2001) that "catch" the polymerase in the act of incorporating the incoming nucleotide; these structures indicate that only one template base becomes part of the ds DNA. It has been argued (Wuite et al., 2000), however, that the n = 2 interpretation remains tenable because not only the first template segment, but also the second one is ordered in the closed state crystal structures and the interphosphate distance in the second is close to that in ds DNA.
Neither model is capable of calculating the effect of the force on the rate-determining step without introducing certain assumptions. To provide first-principle-based information that makes possible a more complete understanding of the experimental results, we use all-atom molecular dynamics simulations and explore the dynamical processes in the open and closed structures as a function of the applied force. Given these results, we evaluate some of the parameters of the rate-force models used to describe experimental data. Specifically, we address the issue of the geometry of the DNA segments at the active site by monitoring distance and angle time series obtained during the molecular dynamics simulations at different values of the applied force, and use these to calculate the force-dependent barrier of the reaction.
In what follows, we begin with a brief review of the global and local models. We then present the results of the molecular dynamics simulations in the presence of an external force and use them to evaluate the parameters that appear in the local model. This permits us to refine the local model and present a restricted-cone local model (RCLM), which takes into account the restrictions on the ss DNA motion due to the presence of the protein. An evaluation of the results is given in the Concluding Discussion section. It is followed by a brief Methods section, which describes the molecular dynamics simulation protocol.
|
THE GLOBAL AND LOCAL MODELS
|
|---|
Global model
The global model (GM) (Wuite et al., 2000; Maier et al., 2000) for the rate-force dependence assumes that the force-dependent part of the activation free energy, i.e., g(f) in Eq. 1, has the form
 | (2) |
The first term on the right-hand side is the force-dependent part of the activation enthalpy, which is equal to the reversible work required to convert n bases from the ss geometry to the ds geometry. The values of xss(f), xds(f), measured in different experiments, are the extensions (along the direction of the force) per base, at force f, for ss and ds DNA chains, respectively; the values are taken from data on stretching of ss or ds polymers (Smith et al., 1996; Maier et al., 2000). The second term, involving the force-dependent entropy of activation, is evaluated, in the GM version of Wuite et al. (2000), from the areas under the experimental curves by plotting the forces fss,ds versus the extensions xss,ds produced by those forces,
 | (3) |
This term is neglected in the GM version of Maier et al. (2000). Both versions of the GM implicitly assume that the ds geometry is formed in the transition state leading to the closed state.
The model is global in the sense that the functions xss(f) and xds(f), which govern both terms in Eq. 2, pertain to the contributions per base to the end-to-end extension of the entire ss or ds polymers on which the pulling experiments are performed. These polymers typically consist of over 10,000 basepairs, whereas the portion of the DNA interacting with the polymerase consists of only 10 basepairs (Doublié et al., 1998; Li et al., 1998b) of the duplex and 4 ss bases of the template (Turner et al., 2003). In effect, data on the elasticity of ss and ds DNA, obtained from experiments in the absence of an enzyme, are used to calibrate the amount of work involved in converting a single ss DNA base to ds geometry, thereby slightly shortening the template. This calibration is then assumed to apply as well in the presence of the enzyme. Accordingly, the activation barrier has no direct contribution from enzyme-DNA interactions; the only parameter having to do with the enzyme is n, the number of ss DNA bases converted to ds geometry in the transition state. Fitting this model to the observed variations of replication rates with tension indicated, on the basis of the calibrated elasticity changes, shortenings that correspond to n = 2 for T7 DNAp (Wuite et al., 2000) and Sequenase (Maier et al., 2000) and n = 4 for the Klenow fragment (Maier et al., 2000).
Molecular dynamics (MD) simulations are well suited to test both conceptual and specific assumptions of the global model. The prime factors employed in the model, xss(f) and xds(f), the projections of DNA segments on the direction of the force, are averages over thousands of residues. Consequently, there is no distinction between the orientation, relative to the direction of the applied force, of the DNA segments at the active site and the orientation of segments far away from the active site. This is inconsistent with a key structural feature, incorporated in the MD simulations. One of the two DNA segments at the active site of the closed complex (see Figs. 2 and 3) is clearly kinked, approximately perpendicular to the axis of the duplex. Because ds DNA is extremely stiff compared with ss DNA, the applied force is essentially along the axis of the duplex, and therefore, being approximately perpendicular to the ss overhang, the resulting torque (that tries to move the ss overhang from a perpendicular to a parallel direction relative to the axis) will contribute significantly to the energetics.
View larger version (53K):
[in this window]
[in a new window]
|
FIGURE 3 (Top panel) Top view (down the axis of the double helix) showing the kinked ss DNA protruding outside the helical boundaries as a consequence of the interactions with the polymerase in both the open (left) and closed (right) structures (same color code as in Fig. 2; protein not shown for clarity). Because tension is applied in the direction coming out of the plane of the article, the resulting torque, acting to align the overhang from an in-plane to an out-of-plane orientation, is expected to have a significant effect on the kinked DNA region. (Bottom panel) Side view of the open and closed structures, (same format as top panel). The direction of the applied force is illustrated by the dotted line. Again, it is apparent that the ss overhang is the "lever arm" of a significant torque contribution forcing the ss segments to orient from a vertical to horizontal position.
|
|
The MD simulations, which evaluate the range of fluctuations in the orientations of the segments at the active site, demonstrate that for both the kinked segment and others, the averaged projections on the force direction differ markedly from the "generic" results given by xss(f) and xds(f). This shows that the basic "calibration" procedure adopted in the global model is not appropriate.
The inference from the global model that n > 1 in Eq. 2 for the three enzymes studied is therefore based on assumptions that are not appropriate, particularly because they do not take into account the importance of the specific (kinked) geometry of the enzyme-DNA interactions.
Local model
The fact that the global model ignores the specific geometry and interactions of the DNA segments at the active site, other than that a number n of bases are converted from ss to ds geometry, led Goel et al. (2001, 2002) to suggest a local, "structurally guided" model for the rate-force dependence. In this model, the force-dependent activation energy depends on the behavior of two DNA segments neighboring the active site. Vectors a and b in Fig. 4 are introduced by connecting equivalent atoms in adjacent bases along the template DNA strand associated with the +1, 0 segment and with the 0, –1 segment (in the local model (LM), the vectors connecting neighboring C1' atoms were used); both the orientation and length of a and b can change in going from the open to the closed state. The additional activation enthalpy barrier in the presence of the force was assumed to equal the average mechanical work done by the enzyme against the external force f, in the process of converting the two segments a, b from their conformation in the open form of the enzyme complex (the reactant state), to their conformation a', b' in the closed complex (the transition state); i.e.,
 | (4) |
where 
indicates an average over the probability distribution functions for the scalar products of f with the base-associated vectors.
View larger version (23K):
[in this window]
[in a new window]
|
FIGURE 4 Views of open and closed Taq polymerase complexes, focused on DNA segments at the exit from the protein. Arrows (a, b in open, a', b' in closed) show the designation of the segment vectors between adjacent C5' atoms and the direction of the applied force f (parallel to axis of the ds DNA helix). Other DNA segments than these two schematized here do not change significantly in going from their open to close states. Portions of the protein (blue) within 10 Å of the segments a or a' are shown; the atoms in this part of the protein are the main contributors to the restriction of orientations available to the DNA segments. The explicit solvent and ions included in the simulation are omitted.
|
|
The local model was introduced to test whether a model that assumes only one base changes geometry from ss to ds, could fit the replication rate data. In the local model it was assumed, in accord with the structural data, that, upon going from the open to the closed complex, the a vector, connecting bases 0 and +1, goes from ss to ds geometry whereas the b vector, connecting –1 and 0, retains ss geometry, although the orientation of b does change (see Fig. 5). Thus only one nucleotide changes its conformation from ss-like to ds-like geometry to be incorporated into the DNA duplex, although the motion of two nucleotides (i.e., 0, –1) is actually involved. Consequently, the distinction between n = 1 as used to describe the LM model (Goel et al., 2001) and n = 2 is somewhat arbitrary and the essential part is that the displacements of two nucleotides are included. The lengths of the a and b vectors (Fig. 2) were constrained to agree approximately with the average contour lengths per residue, Lss = 7 Å and Lds = 2.6 Å, consistent with structural data, as well as the stretching curves used in the global model.
View larger version (38K):
[in this window]
[in a new window]
|
FIGURE 5 Views of first five bases in the open (left) and closed (right) template DNA strand at the active site: bases –1, 0, +1, +2, +3 in blue, red, gray, orange, and yellow, respectively. Clearly, base 0 changes between the open and closed configuration by stacking onto the double-stranded part of the template (i.e., onto base + 1).
|
|
The duplex bound to the polymerase is in B-DNA form except for the three pairs toward the ss junction, which are in A-DNA form. Although the B-form is more stable for free DNA, the DNA segments at the active site of family A polymerases and HIV-1 reverse transcriptase are observed to be in the A-form (Patel and Loeb, 2001). The a segment is thus in A-form, which justifies the A-form value of lds = 2.6 Å rather than the value of 3.4 Å for B-form DNA.
In Eq. 4 the work terms contributing to the force dependence of the activation barrier then become
 | (5a) |
and
 | (5b) |
Here the angles 
and ß specify the orientation, with respect to the direction of f, of the a and b vectors for the open complex; ' and ß' denote the same for the closed conformation (see Fig. 6). The external force is assumed constant and locally directed along the duplex axis, as a consequence of the long persistence length of ds DNA. The averages over the LM angular orientations in Eqs. 5a and 5b correspond to the GM projections in Eq. 2; thus, 
for = , ß, ß' and 
These angular orientations were evaluated using the freely jointed chain (FJC) approximation, which averages over a Boltzmann distribution corresponding to the force-dependent potential energy of the FJC segments (Bueche, 1962). Over most of the pertinent range of f, the FJC results are fairly similar to the experimental xss(f) and xds(f) functions. The FJC angular averages are given by the Langevin formula,
 | (6) |
for 
or ß', where 
with kB the Boltzmann constant, T temperature, and 
or dds the polymer Kuhn length (the Kuhn length is the characteristic length scale describing the flexibility of the chain, and equals twice the persistence length); dss = 14 Å and dds = 1000 Å were used, conventional values (Rouzina and Bloomfield, 2001) in accord with the experimental stretching curves. With the FJC model, simple analytic formulas can be obtained (Goel et al., 2002) for 
analogous to Eqs. 2 and 3, with each the sum of contributions from the a and b segments.
By design, this implementation of the LM facilitates comparison with the GM. In effect, the GM postulates n terms like 
wa(f)
in Eq. 5, corresponding to n segments, each of which shrinks in length 
between the open and closed complex. Instead, the LM assumes that only the leading segment a shrinks, but the neighboring segment b nonetheless can contribute, despite no shrinkage in length, if in wb(f) the averaged angular motion with respect to f differs appreciably between the open and closed complex. This point was illustrated by computing k(f)/k(0) for two limiting cases (Goel et al., 2001, 2002). In Case I, the angular fluctuation of the b segment was considered unhindered in both the open and closed complex, so cos ß = cos ß'; in Case II, the fluctuations remained unhindered in the open complex but are strongly restricted by interaction with the enzyme to be near 90° in the closed complex, resulting in cos ß' = 0. Case I (based on the free energy of activation) gave results nearly identical to the GM with n = 1, as expected, whereas those for case II resembled the GM with n = 2–4.
This LM approach appears to have served its heuristic purpose, demonstrating that the effect of tension on replication rate need not involve "extra" length shrinkage, as inferred from the GM (i.e., n > 1), but might instead be strongly influenced by angular conformational changes induced by the enzyme. The LM in this form cannot otherwise be useful, however. The FJC approximation and other expedient assumptions, although of some use in describing "generic" or global behavior, are not appropriate for local interactions. A realistic version of the LM requires the correct description of the effect of the enzyme-DNA interactions. A purpose of the MD simulations is to address this issue.
Before ending this section, we consider the relation between the LM and GM, and a linear-response model proposed by Bell (1978) and refined by Evans et al. (Evans and Ritchie, 1997). This model has been used for interpreting pulling experiments on ligand-receptor unbinding (Moy et al., 1994), on protein unfolding (Carrion-Vazquez et al., 1999) or DNA unzipping (Strick et al., 2000). In the linear-response model, the amount by which the barrier is changed is g = 
f, where is the width of the activation barrier (i.e., distance along reaction coordinate from the "reactant" minimum to the top of the barrier) and can also be looked at as a characteristic length over which the force acts. This linear form for the barrier change leads to a force-dependent rate equal to
 | (7) |
This is to be compared to the LM result in the limit of large tension
 | (8) |
obtained by adding up Eqs. 4 and 5 in Goel et al. (2001), where 
and 
are positive constants calculated from the model. Equation 8 indicates that, for relatively large forces (higher than 7 pN), the LM is in the linear-response regime, but it includes a constant offset to the barrier change (see Table IV in Goel et al., 2002). In the GM, by contrast, the rate expression, when put in the form of Eq. 7, yields a value for that depends on f as estimated from the experimental curves of free DNA.
|
MOLECULAR DYNAMICS SIMULATIONS
|
|---|
We present a molecular dynamics analysis of the DNA polymerase binding site in this section. By simulating the system at atomic detail in the presence of an external force, we have direct access to the various parameters assumed by the LM, and, as it turns out, are able to introduce a refined version called the RCLM as described in "The restricted-cone local model" section.
This study employed protocols developed for "all-atom" simulations, making full use of the data available from crystal structures for both the open and closed conformations of the Taq polymerase complexed with DNA (Li et al., 1998b) as well as techniques to include solvation and dielectric screening of ionic interactions (as described in Methods). This resulted in a simulation system explicitly treating over 12,000 atoms, of which over 5000 atoms, located within an explicit water sphere of 25-Å radius centered on the binding site, are allowed to move.
The closing of the finger domain is assumed not to be affected by the force; this "first-order" approximation is expected to be reasonable because the force is applied on the DNA and not the protein (Wuite et al., 2000; Maier et al., 2000). In other words, the force-dependent work term depends on f only through the force-dependent properties of the DNA segments (i.e., its orientations) and not through any force-dependent properties of the protein.
Also, in common with the phenomenological treatments, the observed closed structure before the incorporation of the base (with the incoming nucleotide in Watson-Crick pairing with the 0-th template base situated on segment a', see lower panel in Fig. 2), is taken to correspond to the transition state; the open structure corresponds to the reactant state. Consequently, the force-dependent change in the barrier height is calculated between these two states, i.e., between the open and the closed structures. The absence of a force dependence of the protein transition is consistent with the fact that, in the x-ray structures, the ss DNA overhang is free to sweep inside the space of a cone centered at the active site, without steric hindrance from the surrounding protein atoms (see Fig. 1); this unhindered motion has been observed in the molecular dynamics simulations.
The key change of the DNA determining the overall reaction is the 
transformation; i.e., the change in the position of the 0-th DNA base from an ss to a ds configuration, after which the chemical bond-forming reaction is a relatively fast exothermic process. The open and closed crystal structures for Klentaq1 show that this transition involves a rotation of the a segment by 20° toward the helical axis to allow the ss base at position 0 to flip in and hydrogen-bond with the incoming nucleotide. Because the force acts directly on the DNA alone, the assignment of the closed DNA geometry as the transition state geometry is expected to be a good approximation for the calculation of the force-dependent barrier. Even if the transition state for the conformational change were at an intermediate position between the open and the closed states in the absence of an external force, the force-dependent change of the activation energy from the reactant to this transition state could be well approximated by the energy change from the reactant to the closed state. This assumption comes from arguments invoking the Brønsted relationship (Brønsted, 1928), which belongs to the broader class of linear free energy relationships, in which the logarithm of a rate constant is a linear function of the logarithm of an equilibrium constant; a model explaining this observation has been put forth by Evans and Polanyi (1938).
Analysis of origin of force dependence
In Eq. 4, the key variables are the vectors a, b, a', b' (see Fig. 6). They are replaced in Eq. 5 of the LM by the force-dependent average cosines and the lengths of the vectors. In Goel et al. (2001), the average cosines were calculated using a FJC model in an external field, whereas the lengths of the segments were ascribed constant values equal to the interbase distances of ss and ds DNA measured from the crystal structures. Because molecular dynamics simulations allow us to give a time-dependent description of the key variables, we can calculate 
in Eq. 4 directly. To do this, we run constant-temperature molecular dynamics in the open state and in the closed state in the presence of a range of applied forces. This type of simulation, in which the timescale is orders of magnitude shorter than the experimental timescale, has been shown to be useful for mapping out energy profiles for protein folding and unfolding (Paci and Karplus, 2000), ligand-receptor binding (Izrailev et al., 1997; Merkel et al., 1999) or nucleic acid structural transitions (Konrad and Bolonick, 1996; MacKerell and Lee, 1999). To our knowledge, the work presented here is the first simulation of an applied force acting on a protein-DNA complex.
The local force exerted on the DNA in the enzyme complex is assumed to be adequately approximated by the force applied to the entire DNA strand and to be directed parallel to the axis of the double helix portion. Whereas the local instantaneous force acting on a string of arbitrary shape is expected to be directed along the local tangent, we seek to model an average global tension as it is measured in the experiments (and as it is modeled in the GM and LM); its direction is parallel to the double helix axis. Also, the force is assumed to remain constant during the open to closed transition. Actually, the instantaneous force felt by the leading segments at the active site will fluctuate and thereby differ from the global tension applied to the entire DNA strand. The timescale for such fluctuations can be estimated (Goel et al., 2001) from the Zimm model (Grosberg and Khokhlov, 1994) in terms of the Kuhn lengths and the solvent viscosity. This indicates that the timescale for fluctuations of individual Kuhn segments of ds DNA is vastly longer than our computation intervals of 3 ns. Accordingly, the helix axis remains practically stationary in our MD simulations. The stiffness of ds DNA likewise ensures that for it the local and global f are nearly the same (except for very weak forces). The fluctuation timescale for Kuhn segments of ss DNA are, in contrast, substantially shorter than 3 ns. These fluctuations thus are averaged over in the MD simulations. In the simulation, a point was placed along the direction of the double helix axis 40 Å away from the O5' atom of the –1 nucleotide (see Fig. 7) and a constraint force with the desired magnitude was applied to this atom, which is the outmost nonhydrogen atom of the modeled template; i.e., the O5' atom of the b or b' segments. (Details concerning the way the force was simulated are given in the Methods section.)
View larger version (20K):
[in this window]
[in a new window]
|
FIGURE 7 Force direction (shown here in the open state). The incoming nucleotide is in green, the template in gray, and the primer in red. The dotted line points to the position of the fictitious particle used as reference for the application of the force. The direction of the force is parallel to the direction of the DNA helix. For clarity, only the protein within 10 Å of segment a is shown (in blue) and the explicit solvent and ions included in the simulation are omitted.
|
|
Fig. 8 specifies the coordinate system and angles employed in the MD simulations to describe the orientation of the DNA segments. For both the open and closed complexes, the z axis is along the axis of the duplex helix; the polar angle = , ß, ', ß' then specifies the angle between the z-direction and the segment vectors. The azimuthal orientation of the segments about the duplex axis is specified by 
a, b, 'a, 'b, measured counterclockwise from the x axis, which is chosen to lie in a fixed plane that contains the z axis and the most probable direction of the a vector in the open complex. We use the distances between consecutive C5' atoms; this choice is employed in all subsequent analysis reported in this article. The C5' distances are preferable to the C1' distances originally used in the local model (Goel et al., 2001). This is because the C1' atoms are not aligned with the sugar-phosphate backbone to which the pulling force is applied in the simulations and can rotate about that backbone. However, the simulations are rather insensitive to the particular choice of defining points; choosing the distances between the P atoms or the O5' atoms or between the centroid of the sugar-phosphate backbone all give similar results.
We performed two independent runs, Run I and Run II, for forces up to 20 pN, in steps of 2 pN. (The constant-force experiments cover the range 0–20 pN for Klenow and Sequenase (Maier et al., 2000), and 0–15 pN for T7 DNAp (Wuite et al., 2000); some constant-elongation data points are available for up to 35 pN for T7 DNAp (Wuite et al., 2000).) At each force value, 3-ns simulations were run and the last configuration at a given force was used as the starting configuration for the next (higher) force. In these calculations, the length of a(t), b(t) (in the open state), and a'(t), b'(t) (in the closed state) were computed between the C5' atoms of the template nucleotides +1–0, and from 0 to –1, respectively (see Fig. 2). Fig. 9 shows the time series of the angles made by the two DNA segments with the direction of the force, over the range between 0 and 20 pN, and Fig. 10 plots the corresponding lengths of the two segments. We note that, even in the open state, a is much more restricted in its motion (by its interactions with the surrounding enzyme) than b. More importantly, we observe that the angular orientations of a, a', b', are less influenced by the external force, whereas the orientation of b is strongly affected by forces larger than 8–10 pN.
View larger version (64K):
[in this window]
[in a new window]
|
FIGURE 10 Time series for length of DNA segments (a in black, b in red) during the pulling simulations; format as in Fig. 9. Results of run I are shown; run II yields similar plot, with nearly identical mean and standard deviation, and is not shown.
|
|
Fig. 11 plots histograms displaying MD results for the azimuthal angles, defined as in Fig. 8. As with the polar angles, the range of azimuthal orientations is seen to be much more constrained for the a vector than the b vector, in both the open and closed conformations. For forces up to 20 pN (and even at higher tensions) the azimuthal orientations of the a and a' vectors remain within a few degrees of their orientations at f = 0. In contrast, for the open conformation, the most probable azimuthal angle for b, which occurs near 55° for f = 0, shifts markedly to near 90° for f = 20 pN (and remains there at higher forces). In the closed conformation, the histogram of the azimuthal angle for b' is trimodal for f = 0, with a prominent peak near 30° and other sizable ones near 80° and 100°. For f = 20 pN (and higher forces), the most probable azimuthal angle for b' is the most prominent peak, which remains in the vicinity of 30°. These results indicate that in the force range 0–20 pN, the changes in azimuthal orientations of the DNA segments are roughly comparable to those for the polar angles.
View larger version (25K):
[in this window]
[in a new window]
|
FIGURE 11 Histograms for the DNA segment azimuthal angles for small and large applied pulling forces: f = 0 (blue), 20 pN (red), 40 pN (green), and 60 pN (violet); dashed lines for segment a (or a' in the bottom panel), continuous for b (or b' in the bottom panel). The limits of the available azimuthal angles for all simulated forces up to 20 pN indicate that the ratio  is approximately unity and thus there is little contribution to the force-dependent thermodynamic functions from the azimuthal angles.
|
|
Under the assumptions of ergodicity, we can calculate for each force f that we apply in the simulation, the ensemble averaged work w, as required for Eq. 4, from the difference of time averages,
 | (9a) |
 | (9b) |
where the first sum on the right-hand side of each equation is calculated from the trajectory in the open state, and the second sum is from the trajectory in the closed state.
The results of the simulations were used to calculate wa(f)MD and wb(f)MD in Eqs. 9a and 9b. From these values, 
was calculated with
 | (10) |
Fig. 12 shows the various contributions and the predicted values of 
Although there are significant differences between the two runs, particularly for the more flexible open state, the overall trends are clear. For low forces (2–8 pN) the values are negative, corresponding to an increase of rate, but for higher forces is positive and the rate decreases.
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 12 Dependence on pulling force of enthalpy of activation for converting the open-to-closed complex, obtained from MD simulations determining work done by the enzyme on DNA segments. Data represented with lines are averages from Runs I and II; diamonds and squares represent separately h from Run I and Run II, respectively, with the error bars indicating standard deviations. For segments a, a', and b', at forces up to 20 pN, the dependence of the work on force is nearly linear, indicating the orientation of these segments is relatively insensitive to the force; however, segment b undergoes a marked reorientation near 10 pN that shifts the enthalpy of activation from negative to positive values.
|
|
This sign change is due to orientation, relative to the force direction 
of the most mobile segment, i.e., b in the open state. In the crystal structure as well as at low forces, the angle between 
and b is obtuse, so the average projection 
is negative, whereas larger forces are able to orient b to make an acute angle with f so that becomes positive (see Fig. 9). By contrast, for the more rigid closed state, 
has negative signs throughout the entire force range; the orientation of b' is significantly affected only for forces larger than 40 pN, as described in the "Higher force regime" section.
Fig. 13 shows the values of k(f)/k(0) calculated from in Eq. 1 assuming that Ts(f) is negligible, as it appears to be in the LM. The figure also shows the experimentally measured values for the Klenow fragment and the results of the RCLM model described in the next section. It can be seen that there is good agreement with the experimental decrease in the rate for forces between 8 and 20 pN. The apparent increase in the measured rate for two related enzymes is also qualitatively reproduced, although the value and position of the peak in the k–f curve are somewhat different from those of the experiments. As noted in the Introduction, the observed increase in rate at low force is within the experimental error bars; the calculations suggest that the behavior is real.
Fig. 9 leads to the following picture for how the force affects the rate through the orientation of the segments. As explained above, the increase in the rate is due to the change of the sign of 
which is negative for low forces and positive for forces larger than 8 pN. The interactions that keep b at angles ß above 90° in the open state are overcome by forces higher than 8–10 pN. By contrast, in the closed state, the corresponding angle ß' remains at the values close to the zero-force case, for this range of forces. The two distinct behaviors of the angles ß and ß' above 8 pN, combined with the relatively constant contributions of angles and ', yield a force dependence of 
for this force regime that explains the MD-calculated rate-force curve. Moreover, for comparison, our MD estimate of 8–10 pN for the force required to "break" b free from its position in the open structure, is within the range of experimental forces needed to break DNA hairpins (9 ± 3 pN for A-T and 20 ± 3 pN for G-C hairpins; Rief et al., 1999).
Comparison with GM and LM
In Fig. 14 we contrast the variation with force of the average cosines determined from the MD simulations with experimental force-extension data, x(f)/L, used in the GM and with the similar FJC curves from Eq. 6 used in the LM treatment. Such functions, in effect, restrict to acute angles and typically give values of cos greatly different in magnitude and/or sign from the MD results; in Case I of the LM, the segments are free to rotate in a solid angle of 4
, so their most favorable orientation is along f (angles of 0° are most favorable), which yields large average cosines for a, b, a', and b'. By contrast, in the force regime between 0 and 20 pN, the MD values are always very small (cos 
.15). In Case II, the average cosines for a, b, a' are large (angles close to 0°), whereas the average cosine for b' is postulated by the LM to be 0 (due to what is called in the LM the kink, which keeps the respective angle at 90°). There is a similar contrast with respect to segment lengths. In the GM and the LM, shrinkage accompanying the conversion of an ss segment to a ds segment has a major role; the values used in the LM are Lss = 7 and Lds = 2.6 Å (see above). However, whereas the mean interbase distances do show such shrinkage for free DNA strands, the differences obtained from MD are much smaller and even opposite in sense for segments near the active site when DNA is in an enzyme complex. In the MD simulations the DNA segment lengths do not differ much between the open and closed complexes; moreover, the nominally ss segments, b and b', are actually shorter than the a and a' segments. The agreement of the MD results and the experimental rate-force curve is due to the force-dependent orientation of the segments in the open relative to the closed state. In particular, the main contributor for the 0–30 pN range is the b segment, which due to its orientation, makes a negative contribution for the force-dependent barrier for forces lower than 10 pN and a positive contribution, for forces larger than 10 pN, which orient this segment along the direction of the applied force. These comparisons show that the apparent agreement between the GM or LM and experiment must be attributed to compensating errors in the distance and angular parameters used in the models.
Higher force regime
Although the experiments show that the rate of incorporation goes to zero at 25 pN, we extended the simulations to higher forces. A second domain of large force-induced change was found in the MD simulations extended for forces above the stalling force for polymerization. Results obtained for the polar angles, using conditions of Run I with 3-ns simulation intervals, are shown in Figs. 15 and 16 for forces up to 80 pN. In the closed complex, the segment b' undergoes a major excursion: ß' first being obtuse, reaching values around 160° by 30 pN, abruptly plunges by 80° to become acute by 40 pN. Moreover, in both the open and closed complex, at 40 pN the other segment also shifts significantly: and ' both become more acute, shifting closer to the direction of f by 10°.
View larger version (34K):
[in this window]
[in a new window]
|
FIGURE 15 Time series of the angles ß and ß' describing the orientation of the last ss DNA segment, b and b', respectively, during the pulling simulations at large forces (30–60 pN), for the open state (top panel) and closed state (bottom panel). For the open state the segment b has oriented toward the direction of the external force at smaller forces (between 8 and 10 pN; see Fig. 9). In the bottom panel of this figure, for the closed state, the evolution of the ß'-angle shows that b' experiences this "crossover" phenomenon at forces larger than 40 pN.
|
|
The striking polar angular shifts seen at 40 pN, especially for ', suggest a possible connection to exonucleolysis. In the pulling experiments on T7 DNAp, the polymerase activity was observed to stall at 34 ± 8 pN and a force-induced 100-fold increase in exonucleolysis activity began above 40 ± 3 pN (Wuite et al., 2000). Although the Taq polymerase has no detected exonuclease activity, this activity is well characterized for analogous family A polymerases. Exonucleolysis occurs at a different site that, in the Klenow fragment (which has the same overall structure as Klentaq1, 50% homology and 36% amino acid identity) is 30 Å away from the polymerization site (Beese et al., 1993). Thus, it is of interest to consider whether changes in the geometrical alignments at the polymerase site for Klentaq1 at high force may be akin to those that foster exonucleolysis. Improper orientation of the a' and b' segments in the closed conformation could trigger the 3'–5' proofreading activity of the enzyme. This is in accord with current views of the preservation of nucleotide insertion fidelity, which emphasize geometric selection by induced fit at the active site (Steitz, 1999; Kunkel and Bebenek, 2000). Particularly important for the chemical reaction of incorporating the incoming nucleotide is the orientation of the a' segment in the closed structure. A value of ' around 72° (distributions to the right in the lower panel of Fig. 16) corresponds to the correct positioning for the polymerization step, whereas values around 60° (distributions to the left in Fig. 16) are suboptimal and could trigger the exonucleolysis mode. Two-state behavior consistent with this bimodal hypothesis has been observed for T7 DNAp; near the stalling force, a competition between polymerization and exonucleolysis was observed (Fig. 5 in Wuite et al., 2000).
In summary, the MD simulations reveal how the two leading DNA segments (of which only one changes its geometry from an ss to a ds one) at the polymerase active site respond to competition between enzyme-DNA interactions and the external force. The segment lengths, corresponding to interbase distances between adjacent nucleotides, change by 7% between the open and closed conformations of the enzyme, but remain virtually unchanged by external forces up to 80 pN. The segment angular orientations with respect to the duplex axis direction (along which the external force is exerted) undergo large changes, in response to both the open to closed transition and to the external force. In contrast, the azimuthal orientations about the duplex axis change greatly in the open to closed transition, but are almost unaffected by the external force. As displayed in Fig. 17, the major effect of the applied force occurs for the outermost DNA segment. In the open complex, a force of 10 pN suffices to tilt b by 45° toward f, whereas in the closed complex, increasing force at first has no effect on b'. However, at 40 pN there is an abrupt change with the force reorienting b' by 80° toward f. Above 40 pN, in both the open and closed complexes, the inner segment a is also shifted further toward f by 10°, as seen in Fig. 16. These shifts, which change the geometry of the binding site, might inhibit base incorporation and foster exonuclease activity.
View larger version (22K):
[in this window]
[in a new window]
|
FIGURE 17 Representation of all polar angles for the open (in black and green) and the closed (in red and blue) states as a function of the force applied in the molecular dynamics simulations; error bars (standard deviations) are shown for several data points. Notice the abrupt decrease to acute values of ß at 10 pN and of ß' after 30 pN.
|
|
|
THE RESTRICTED-CONE LOCAL MODEL
|
|---|
Although in the LM all orientations of the DNA segments considered in the model are assumed to be accessible, the molecular dynamics results show that steric constraints exclude a significant amount of the space in both the open and the closed structures. To take account of this fact in a simple way, we introduce the restricted-cone local model. As in the LM, the RCLM focuses on the DNA vectors (a, b in the open and a' and b' in the closed states), which are treated as "dipoles" orienting in an external force field f and a heat bath of temperature T (as in conventional FJC models). The essential point of the RCLM, as implied by its name, is that, due to the presence of the protein, not all angular orientations 
(, ) are allowed for the DNA segments, as illustrated in Fig. 18. By writing down the partition function for the protein-restricted segments in the presence of the external force, we derive the corresponding enthalpy and free-energy barriers. The LM is recovered as a particular case of the RCLM by removing the angular restrictions.
View larger version (77K):
[in this window]
[in a new window]
|
FIGURE 18 Schematic representation of the RCLM components modeling the spatial restriction shown in Fig. 4. The shaded region depicts the volume excluded by the presence of the enzyme and that restricts the orientation available for the a and b vectors. Regions A and B contain the rest of the DNA degrees of freedom that the RCLM is neglecting: for example, region A contains the guanine base that, when the reaction proceeds toward the closed state, flips in to pair with the incoming cytosine nucleotide.
|
|
To introduce the RCLM, we assume that the potential energy of a segment d (d = a, b, a', b') in the presence of the protein
TOP
ABSTRACT
INTRODUCTION
) and Run II (X), 
) of DNA segments