| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
* Beckman Institute, Department of Physics, and Howard Hughes Medical Institute, University of Illinois at Urbana-Champaign, Urbana, Illinois
Correspondence: Address reprint requests to K. Schulten, Tel.: 217-244-1604; E-mail: kschulte{at}ks.uiuc.edu.
 |
ABSTRACT |
|---|
|
INTRODUCTION |
|---|
–5). The helicases are involved in almost all aspects of DNA metabolism, including transcription, replication, and recombination. Defects in helicase functioning in humans can lead to genomic instability and predisposition to cancer (6). To achieve their functions, helicases move in a unidirectional manner along single-stranded DNA (ssDNA) (7); when a helicase continues its translocation along ssDNA encountering a junction formed by duplex DNA, the duplex DNA becomes unwound.
In their functional forms, helicases assemble as hexamers, tetramers, dimers, or monomers (1,3,4). PcrA helicase from Bacillus stearothermophilus (B. stearothermophilus) has been proposed to work as a monomer (8). Belonging to the superfamily 1 (SF1) helicases (5,9), monomeric PcrA (
80 kDa) is composed of four domains (1A, 2A, 1B, and 2B), resembling other SF1 helicases in their monomeric forms, e.g., Rep and UvrD, although Rep and UvrD are oligomers in their functional forms (10–12). The PcrA, Rep, and UvrD monomers exhibit 40% sequence identity, and all have been demonstrated in experiments as being capable of translocating progressively 3' to 5' on ssDNA (13–16).
PcrA helicase from B. stearothermophilus has been crystallized and resolved at high resolution (8), in both a substrate (with ATP bound) and a product (without ATP/ADP bound) state. The structures, as shown in Fig. 1, were crystallized in the presence of a DNA junction, i.e., duplex DNA flanked by a piece of 3' ssDNA, with the duplex bound to the side of domain 2B and an elongated piece of ssDNA crossing above the two RecA-like domains 2A and 1A. Domains 1A and 2A, conserved among a class of helicase-like proteins (2,9), are proposed to play the most essential role in the translocation. ATP binds into a cleft between domains 2A and 1A, the binding site being lined by amino acids that are highly conserved among SF1 helicases (5,9). Two of the motifs in the ATP binding site (Walker A and Walker B) are highly conserved among all ATPases. Indeed, superimposing the ATP binding pocket of PcrA helicase with that of F1-ATPase shows that the ATP binding sites have high structural identity, suggesting that a closely related ATP hydrolysis mechanism may be at work (17).
|
), Aksimentiev et al. (19), Ma et al. (20), and Wang and Oster (21) for molecular machines. PcrA helicase is selected here, since it is one of the smallest linear motors with full atomic scale structures available as well as experimental information on velocity and step size. The helicase system also involves protein-DNA interaction and recognition in a very confined space. Understanding the basic mechanism of PcrA helicase may facilitate understanding of more complex molecular motors.
Prior theoretical work investigated helicase function in a generic framework employing certain mathematical models that can describe helicase unwinding duplex DNA (22–24). This work proved very useful for this study, yet it postulated ad hoc the fundamental steps in helicase motor function, namely, different forward and backward translocation rates. The authors also addressed only generic helicases, not any particular one. In contrast, the present work seeks to establish the helicase motor mechanism from the structural and physical properties of a particular helicase, PcrA, the relevant physical properties being established through molecular dynamics (MD) simulations.
We based our study on the crystallographic structures reported in Velankar et al. (8) and mentioned above. We assume that the two structures, one with an ATP (analog) bound and one without ATP/ADP, present key states lying along the pathway of the PcrA translocation process. This is by no means certain since artifacts, in particular due to crystal packing, might shift the protein away from its mechanistically relevant conformations. The goal of our study was to identify through molecular dynamics simulations the mechanism of PcrA translocation. Two main translocation models had been discussed in the literature, an active rolling model, suggested actually for a dimeric helicase (25), and an inchworm model (26). As long as PcrA translocates as a monomer, the inchworm model is presently the only candidate. Our study assumes that PcrA works as a monomer and, hence, it focuses only on the inchworm model. This model was also proposed by the crystallographers who solved the structure of PcrA (8). According to the model, PcrA moves by alternating affinities between its translocation domains (2A and 1A) and ssDNA. So far, however, the inchworm model, while eminently insightful, was based on intuition, i.e., it is mainly qualitative and does not result from quantitative physical properties of PcrA derived from its crystallographic structures. Our study seeks to provide the missing physical basis for the inchworm model.
Experimental data showed that the translocation speed of PcrA along ssDNA is 50 nucleotides (nt) per second, presumably consuming one ATP for 1-nt distance (13). Based on this information, we propose a schematic model (see bottom left panels of Fig. 1, a and b), in which the substrate state exhibits lower energy barriers for domain movement of 2A along ssDNA (red curve) and higher ones for 1A (green curve), while the product state exhibits lower energy barriers for domain movement of 1A along ssDNA and higher ones for 2A; coupling these changes in mobilities to attraction (upon ATP binding) and repulsion (upon ADP+Pi dissociation) between the two domains leads to directed translocation. This model followed from the molecular dynamics (MD) simulations reported here, from which we derived the potential for individual domain (2A and 1A) motions along ssDNA as shown in Fig. 1. For this model we develop then a stochastic dynamics description of PcrA translocation. Our study suggests that the unidirectional translocation (3' to 5') of PcrA is a direct consequence of alternating high and low energy barriers experienced by 2A and 1A during each ATP hydrolysis cycle. Our work provides microscopic justification for the alternating affinities between domains and ssDNA proposed by the inchworm model (8).
In the following, we first introduce in Methods the stochastic model for the millisecond domain movements of PcrA. Then we describe a scheme to determine from MD simulations the potentials that govern the domain movements. The potentials obtained from MD simulations are provided in Results, followed by the demonstration that the potentials coupled to the ATP hydrolysis cycle lead to unidirectional movement. We finally identify the key amino-acid residues coordinating the unidirectional translocation in PcrA helicase, and provide further evidence for dynamic asymmetries of PcrA domains with bound ssDNA.
|
METHODS |
|---|
Stochastic dynamics description of PcrA translocation
Here we seek to describe the translocation of PcrA along ssDNA by means of stochastic dynamics theory. For this purpose, we will assume two limiting scenarios, expecting that the most realistic model falls between the two limits. We envision that the sliding of PcrA along ssDNA comes about through an inchworm motion involving separate, but coupled translocations of its 2A and 1A domains as suggested in Velankar et al. (8). Three factors govern the linked motion of domains 2A and 1A:
Below we introduce two scenarios that demonstrate how 1–3 can endow PcrA with unidirectional motion.
The motion of domains 2A and 1A translocating along ssDNA (coordinates x1 and x2, respectively) is a stochastic process, that, on the relevant timescale, can be described by means of a Langevin equation in the strong friction limit (27,28),
 | (1) |
Equation 1 holds in a particular state of PcrA, e.g., the s or p state. Here 
is the friction coefficient; 
is the fluctuating force, represented through so-called Gaussian white noise (29,30); W(x1, x2) is the potential governing the movement of domains 1A and 2A along ssDNA; xi, i = 1, 2, are defined through
 | (2) |
and 
are the coordinates of the centers of mass of domains 1A and 2A along ssDNA, 
p denotes the average in the p state, and l0 is defined as 1-nt distance (6.5 Å). We note that –x defines the forward, i.e., toward the DNA junction, direction. The Langevin equation, Eq. 1, ignores possible memory effects; at this exploratory stage of the investigation the neglect of memory effects seems to be justified. The approach, of course, would be wrong for extreme memory effects, e.g., if translocation in step 1, 3, ... follows a different mechanism from translocation in steps 2, 4, . The dynamics adopted ignores also hydrodynamic effects of the domain motions and translocation.
We define potentials, Ui
(xi), governing individual (i = 1,2 for 1A and 2A, respectively) domain motions in the s and the p states ( = s, p); the interaction potential between the two domains is written V'(x1, x2) (' = s, p). The values and ' are independent indices, i.e., below we will combine Ui and V' with different indices (states) , '. The value W(x1, x2) in Eq. 1, labeled by indices ', is then decomposed into three contributions:
 | (3) |
(xi) are derived from MD simulations (described below) and are shown in insets of Fig. 2.
|
'(x1, x2) = 0. Therefore, there are only two independent Langevin equations (Eq. 1) for the p and s states, respectively. Geometrical constrains (1) (see above) still apply, e.g., it holds that 0 < x1 – x2 
l0; the condition is satisfied through reflection boundaries. Starting from the p state, once 1A gets close enough to 2A, e.g., x1 – x2 < l0/3, PcrA can change rapidly to the s state (ATP bound); starting from the s state, once 2A moves far enough from 1A, e.g., x1 – x2 > 2l0/3, PcrA can change back rapidly to the p state (no ATP/ADP bound).
In the second scenario, referred to as the strong coupling case, the domains are pulled together when ATP binds and then are pushed apart when ADP and Pi unbind by means of the interaction potential V'(x1, x2). This potential is modeled through
 | (4) |
' for ' = s, p is chosen as lp = l0 and ls < lp, e.g., ls = l0/3. The form of the potential in Eq. 4 is the simplest choice and has been adopted due to lack of more detailed information. To improve the description, one might sample intersubunit distances for both crystal structures, but such sampling is presently unfeasible given the expected microsecond-to-millisecond timescale of the domain motion.
As mentioned above, ' (in V') and (in Ui) are independent indices; hence combinations of them (in W') can represent four different states. In the equilibrium substrate (or product) state = ' = s (or p) holds, which will be labeled ss (or pp). For = s and ' = p, the system is in transition from the equilibrium substrate to the product state. We will refer to the intermediate state as sp; similarly we call the intermediate state from the product to the substrate state ps ( = p and ' = s). Correspondingly, there are four independent Langevin equations (Eq. 1) for these four states.
Transition from the pp state to the intermediate ps state is triggered by arrival of ATP and transition from the ss state to the intermediate sp state is triggered by the formation of ADP+Pi (through ATP hydrolysis). The two transitions happen at certain rates as specified below. For transitions from the intermediate ps or sp state to the equilibrium ss or pp state to happen, geometrical criteria were applied, such that ps transits to ss when the domain separation x1–x2 shrinks below ls while sp transits to pp as x1–x2 elongates beyond lp.
The Langevin equation (Eq. 1) can be solved numerically (28) when one assumes discrete time-steps 
t, adopting the scheme
 | (5) |
Here Z is a normal (Gaussian) random variable (29) (with mean 0 and variance 1); D is the diffusion coefficient, according to the fluctuation dissipation theorem (30) related to the friction coefficient through D = kBT/; and D is expected to assume a value 104 Å2/µs (corresponding to a typical diffusion coefficient of a 3 nm-radius protein in solution (27)). Coupling between equations is achieved through a random process which we describe now.
In case of weak coupling, the transitions between states s and p (right after the domain motion) are assumed to be fast compared to the domain motion described by Eq. 1. Once x1–x2 satisfies the criteria for state transitions mentioned above, transitions are induced through a Poisson process using rate constants specified below. In the strong coupling case, the transitions between states ss and sp or between states pp and ps are assumed to be slow compared to the stochastic motion in Eq. 1 and are also described through a Poisson process with rate constants specified below. Poisson processes are simulated by generating uniformly distributed random numbers Y (Y 
[0,1]) and adopting the transition in the case Y 
t (t << 1), where is the rate constant for the transition and t is the discrete time step in Eq. 5.
The stochastic method adopted is related to the so-called kinetic and dynamic Monte Carlo scheme widely adopted in physics and chemistry (31–34). The simple choice of reaction coordinate x (path along ssDNA) could be improved by determining a reaction coordinate accounting for the forward reptation of ssDNA in PcrA by means of a so-called reaction path method (for a review see (35)), e.g., the ones suggested in Elber (36) and Straub (37).
Nucleotide binding site energies and the site energy function
A key task in our study is the derivation of the potential Ui(xi) in Eq. 1, governing the movement of domains 1A and 2A. Since the translocation of helicase along ssDNA arises through ssDNA units (nucleotides) binding and unbinding sequentially to the domain surfaces, one strategy for determining the potential is to calculate the binding free energies Eb for individual nucleotide binding sites and then use these energies to obtain Ui(xi). The calculation of the binding free energy from MD simulations can be achieved by a method (38,39) which is based 1), on the linear response approximation for electrostatic forces; and 2), on linear scaling between solvation energies and average van der Waals (vdW) energies. The method evaluates the absolute free energy solely from the difference between the average interaction energy between nucleotide and protein+water+ion (bound state) and the average interaction energy between nucleotide and water+ion (free state). The value Eb, for a given binding site located at xi, is given as a weighted sum of the electrostatic and the vdW interaction energy,
 | (6) |
and ß discussed further below. Here denotes the difference between average energies (stemming from simulations) of the bound and the free states, i.e.,
 | (7) |
In the present case, we will need only relative energies, not absolute ones, e.g., Eb(xi) – Eb(xj). Since, in the free state (nucleotides with water + ion, but without protein), the nucleotide-solvent interaction energies Eele/vdW(free, xi) should not vary along ssDNA (poly-thymine), i.e., be independent of xi, it is not necessary to actually calculate Eele/vdW(free, xi) and we set these energies to zero. We introduce a further approximation in evaluating Eele/vdW(bound, xi). For this purpose we split the energy terms as follows:
 | (8) |
For nucleotides buried inside PcrA, we assume that the contributions 
are independent of xi; therefore, we also set these energies to zero and count only nucleotide-protein interactions in evaluating Eb(xi). On the side of the protein all amino acids were added in calculating the relevant energies.
The weight coefficient in Eq. 6 was shown to be 0.5 (38,39) while ß could adopt values from 0.15 to 1.0, depending on the hydrophobicity of the binding sites (the more hydrophobic a site, the larger ß) (40). Since the (ssDNA) nucleotide binding sites are buried inside the protein, i.e., are less exposed to water, we adopted ß = 1.0 (slight variation of ß does not affect our major conclusions). The energies Eele and EvdW were calculated (with a cutoff distance of 60 Å) from MD trajectories, sampling every 100 ps and averaging over 2 ns, beginning after the first 1 ns of equilibration.
The weighted sums of Eele and EvdW for individual ssDNA nucleotide binding sites, calculated according to Eq. 6, are presented in Fig. 2 by discrete site energies Eb(xi), where xi, i = 15/14 (in s/p), ..., 21, are positions of corresponding nucleotides (numbered in the same way as in Fig. 1). To relate the energies Eb(xi) to a process in which ssDNA nucleotides move collectively, in a continuous way, across the surface of domain 1A or 2A, one needs a continuous site energy function Eb(x) connecting the discrete Eb(xi) values, x being the path length along an imaginary line going along the backbone of the ssDNA bound by PcrA. Since the Eb(xi) should represent energy minima or at least marginally stable points for the ssDNA nucleotides, the functional values of a continuous energy function Eb(x) between points xi represent the roughness of the energy surface, the roughness determining the overall speed of ssDNA motion and, hence, of PcrA. The roughness effect can be adjusted effectively at a later stage of the calculation (through a parameter 
, see Eq. 13) and, hence, we assume Eb(x) to be represented by a smooth curve interpolating the Eb(xi) values. Accordingly, Eb(x) was constructed using a third-order polynomial interpolation; the construction scheme is detailed in Appendix A. The resulting curve is shown in Fig. 2. As one can see, the construction assigns to the boundary sites (exposed to solvent) x15/x14 (in s/p) and x21 equal energy values.
The interpolation scheme constitutes a significant assumption in our description. One may suggest to replace Eb(x) by a potential of mean force derived through umbrella sampling (41) or steered MD (SMD) (42) linked to use of the Jarzynski identity (43,44). However, such an approach is unfeasible because of the complex degrees of freedom orthogonal to x; these degrees of freedom might participate also very specifically in the translocation, e.g., through base flipping—requiring then a new reaction coordinate and in any case, would be slow to relax.
Potentials governing individual domain motions along ssDNA
With the knowledge of Eb(x) one can estimate the potentials Ui(xi) governing individual motions of individual domains. We assume that PcrA translocates during one ATP hydrolysis cycle one nucleotide (nt) as suggested by experiment (13). Then the total energy of PcrA moving along ssDNA (to which Ui is but one contribution) should have period-one (1-nt distance); we will impose, however, the more stringent condition that both U1 and U2 have period-one.
We want to estimate now Ui from the binding energies Eb(x) considering first the s state. When domain 2A moves forward (to the left in Fig. 1 a), nucleotides 19 and 20 (and those beyond 20) remain in the same binding sites; however, nucleotide 15 will move toward 16, 16 toward 17, then 17 toward a new site, A (which would be occupied by 17 in the p state; see Fig. 1 b), and 18 will move toward another new site, B (which would be occupied by 18 in the p state; see Fig. 1 b), anchored into the pocket formed by side chains Tyr-257 and Phe-64. The potential energy governing this motion of 2A, U2s, is equal to the sum of site energy differences connected with moving nucleotides 15, 16, 17 and 18 (and those beyond 15, which along with 15 are assumed to be of equal energy) backward, as indicated by colored arrows in Fig. 2. This energy is
 | (9) |
In our description we assume that U2s(x) adopts a symmetrical form around x = 1/2, the position of an energy barrier separating the states before (x = 0) and after (x = 1) an ATP hydrolysis half-cycle. This assumption also applies to other Ui(x) derived below.
Now we can determine the boundary energy at site x15. The emergence of the two new sites at x = 1 (as the system moves to the p state), A and B, arises effectively from the disappearance of the middle site at x18 and the boundary site at x15; therefore, it holds that 
, with 
and 
being the (unknown) binding energies at sites A and B, respectively. Since A and B are sites close to x17 and x18, 
might be approximated by Eb(x17) + Eb(x18); the periodicity condition U2s(0) = U2s(1) = 0 thus yields Eb(x15) = Eb(x17).
One can apply a similar reasoning to the backward (to the right in Fig. 1 a) motion of domain 1A and derive Eb(x21) = Eb(x17). In this case, nucleotides 16 and 17 (and those to the left beyond 16) will occupy their original binding sites while nucleotide 21 moves toward 20, 20 toward 19, 19 toward the new site B, and 18 toward the new site A. Thus, one can conclude
 | (10) |
Similarly, we consider the energies for the p state, i.e., Uip. When domain 2A moves backward (to the right in Fig. 1 b), nucleotides 19, 20 ... remain at the same binding sites while nucleotide 15 moves toward 14, 16 toward 15, 17 toward 16, and 18 toward the new site C (which was occupied by 18 in the s state, see Fig. 1 a). Following the reasoning for the s state above, the energy U2p(x) is
 | (11) |
The periodicity condition 
stipulates 
. With 
approximated by Eb(x17), this gives Eb(x14) = Eb(x18). For a similar reason this condition applies also to the case that domain 1A moves forward with its energetics described by U1p. In that case, nucleotides 15, 16 occupy their original binding sites while nucleotide 17 moves toward the new site C, 18 toward 19, 19 toward 20, and 20 toward 21. The energy is
 | (12) |
Combining the definitions of Ui(x) above, i.e., Eqs. 9–12, it can be recognized that each barrier Ai, defined as Ui(x) at x = 1/2, can be written
 | (13) |
i is a combination of site energy terms, stated explicitly in Eq. 14. The site energies terms are determined from MD simulations as explained above. We add here a term 4 that controls the speed of PcrA translocation, i.e., 1 nt per 20 ms, consistent with observation (13). The relationship of , which is chosen identical for both domains (i = 1, 2) and for both states ( = s, p), to the energy surface roughness (energy undulations) of Eb(x), is stated in Appendix A.
We emphasize that the construction and later use of the potentials Ui(x) hinges on particular motions of ssDNA along the domains as specified above. If the actual translocation would involve different behavior of ssDNA, the description adopted here would be called into question. The results of our study should be interpreted as proving the feasibility of the inchworm model once it is assumed.
Molecular dynamics simulations
Starting from the crystal structures of the PcrA helicase in the s (PDB code 3PJR) and p (PDB code 2PJR) state (8), we added missing residues in the protein as well as elongated both duplex DNA and ssDNA (poly-T). The ssDNA was elongated by first opening, i.e., swiveling by 30°, the 2B domain (45), adding the corresponding piece of ssDNA (five nucleotides) from the structurally homologous Rep complex (PDB code: 1UAA) (45) to the end of the original ssDNA, then swiveling the 2B domain back to its initial position by means of SMD (42), and letting the elongated ssDNA segment relax to fit to the PcrA domains.
The structures of the PcrA-DNA complex were solvated in a box of explicit water. Sodium, magnesium, and chloride ions were added to neutralize the negative charges of the PcrA-DNA complex and to control the ionic strength (0.1 M). The program Delphi (46) was used to locate positions of minimal electrostatic energies around the complex, where the water molecules were then replaced by ions. The whole simulated system of protein, DNA, water, and ions contained 110,000 atoms.
All simulations used the program NAMD2 (47), the CHARMM27 force field (48), with an integration time step of 1 fs, and periodic boundary conditions. VdW energies were calculated using a smooth (10–12 Å) cutoff. The particle-mesh Ewald method (49) was employed for full electrostatics, with the density of grid points at least 1/Å in all cases. The particle-mesh Ewald electrostatic forces were computed every four time steps. The simulations were performed in the NpT ensemble, using the No
e-Hoover Langevin piston method (50,51) for pressure control (1 atm), with an oscillation period of 100 fs and a damping time of 50 fs; Langevin forces (52) were applied to all heavy atoms for temperature control (310 K) with a coupling coefficient of 5/ps.
To verify the relationship between barrier heights of individual domain motions, i.e., A2s < A1s and A2p > A1p, which were derived from equilibrium MD simulations, we carried out an independent set of SMD simulations pulling ssDNA one half-step (3 Å, corresponding to x = 1/2 in Eqs. 10–13) forward and backward along the ssDNA binding interfaces of 2A/1A in both the s and p states. For this purpose, we attached 10 harmonic springs (force constants of 2 kcal/mol Å2) to 10 phosphorous atoms of ssDNA, and pulled the ends of the springs with a constant velocity of 3 Å/ns. The SMD methodology is explained and reviewed in Isralewitz et al. (42), Izrailev et al. (53), and Sotomayor et al. (54). Each pulling simulation was repeated four times and measurements of SMD forces were averaged over the four trajectories. Sample movies of the SMD simulations are provided in Supplementary Material.
To estimate the force constant k arising in the interaction potential, Eq. 4, we monitored in our MD simulations the fluctuations of the distance L between the center of mass of two groups of atoms that represent the edges of domains 2A and 1A. We chose for these groups a helix segment (residues 369–374) from domain 2A and a helix segment (residues 76–81) from domain 1A. The MD simulations revealed a standard deviation L 1Å in both s and p states. According to the Brownian oscillator relationship k = kBT/L2 we estimate k 1 kBT /Å2. The simulation data are provided in Supplementary Material.
|
RESULTS AND DISCUSSION |
|---|
3 ns. The root mean-square deviation curves for backbone atoms of protein and DNA during each equilibration are shown in Fig. 3, demonstrating that the system was more or less stabilized after 1 ns.
|
By monitoring an angle formed by two protein helices, residues 360–374 from 2A and residues 66–81 from 1A, during equilibration, it was found that the ATP binding cleft in-between domain 2A and 1A is 10° wider in the p state than in the s state. The opening and closing of the ATP binding cleft in different states (binding of ATP favors a narrower separation between domains 2A and 1A while unbinding of ADP and Pi favors a wider separation between the domains) turn out to be the most important structural feature directly supporting the PcrA translocation, as discussed below.
Potentials Ui from MD simulation
Utilizing the equilibrium MD trajectories, we obtained the binding energies Eb(xi) of individual ssDNA nucleotides. The results, based on Eq. 6, are shown in Fig. 2 for the s and the p state, specifying both the electrostatic and the vdW contributions. It can be recognized that electrostatic energies, only due to the protein, differ significantly for individual nucleotides between s and p states, whereas the vdW energies do not. The standard deviation of the average value of the electrostatic (vdW) energy for each nucleotide is 2–4 kcal/mol (0.5 kcal/mol). For each state, the continuous site energy function Eb(x) (also shown in Fig. 2) connecting the discrete binding free energy values was constructed.
Eb(x) permits one to derive the potential Ui(xi) governing the motion of individual domains along ssDNA (see Eqs. 9–12). The potential Ui(x) is shown in the insets of Fig. 2. One can recognize that in the s state the barrier height (A2s) in the potential U2s for 2A motion is lower than the barrier (A1s) in U1s; the opposite is true in the p state. The barrier heights Ai can be expressed explicitly in terms of the site energy differences of nucleotides
 | (14) |
is a tunable parameter (see Appendix A). Although the barrier heights Ai individually depend on , the difference between barrier heights of competitive (1A vs. 2A) domain motions in each state is independent of , i.e., the value A1s – A2s 9 kcal/mol in the s state and A2p – A1p 12 kcal/mol in the p state are independent of . The difference between barrier heights actually governs the timescale separation between the competing domain motions. From Eq. 14, one can recognize that A2 – A1, after canceling identical terms in A2 and A1, is determined by energy imbalances (combined energy differences between neighboring sites) among some specific nucleotide binding sites, which are indicated in Fig. 2.
This difference between barriers A2 and A1 for domain motions is a key result of our study and was verified by us. In our analysis we accounted only for protein-DNA interactions without explicit solvent contributions. To partially verify the result we carried out an independent set of SMD simulations (sample movies of the SMD simulations are provided in Supplementary Material), pulling ssDNA one half-step forward and backward, along the protein-ssDNA interface across the helicase domains. The forces needed to pull the relevant nucleotides were monitored and are reproduced in Fig. 4. The results show that in the s state, the average force needed to move nucleotides 15–18 rightward (corresponding to 2A motion toward the left) is smaller than that needed to move nucleotides 18–21 leftward (corresponding to 1A motion toward the right); vice versa, in the p state, the average force needed to move nucleotides 15–18 leftward (corresponding to 2A motion toward the right) is larger than the force to move nucleotides 17–20 rightward (corresponding to 1A motion toward the left). This is consistent with the Ui(x) potentials shown in Fig. 2. We note that the SMD simulations took into account interactions of ssDNA with both protein and solvent (water/ions) as well as the short time part of the entropic effect.
|
20 ms (13). Our MD simulations cover, however, only nanoseconds. As explained above, the timescale gap can be overcome through a stochastic, i.e., Langevin, dynamics description that builds on the potentials Ui (see Eq. 3 and 9–12) estimated from the MD simulations. As explained in Methods, we simulated the millisecond domain movements of PcrA through Eq. 1 in corresponding states employing the potential function defined in Eq. 3. In case of the weak coupling scenario we choose for V' a vanishing potential; in case of the strong coupling scenario, we describe V' through Eq. 4. Our description includes random transitions between states. The results of this description are provided in Fig. 5.
|
|
70 Å, reflecting 11 translocation steps with variable durations (due to the stochastic nature) in between. One can recognize that the domains move in a nearly synchronous fashion along the –x-direction (see movie provided in Supplementary Material). Underlying this motion is the scenario shown in the right panel of Fig. 5 a. The system traverses sequentially configurations 1' 
2' 3' 4' 1'' as identified in the figure. The configurations are also identified through state labels p and s, introduced earlier. In configuration 1', PcrA is in the p state, domains 2A and 1A are separated by 1-nt distance and move according to potentials U2p and U1p, respectively; 1A experiences a low barrier and can move readily, while 2A experiences a high barrier and is essentially stuck. When 1A has moved forward (to the left) close enough to domain 2A (e.g., by l0/3, see Methods) it reaches configuration 2'. In this configuration the system has a high probability to transit to the s state (in reality this corresponds to the approach of domains 2A and 1A, leading to an optimal binding geometry for ATP), reaching configuration 3'. In the s state, however, the potentials U2s and U1s differ qualitatively, in that now 2A is easy to move and 1A becomes stuck. When 2A moves forward (to the left) far enough from domain 1A (e.g., by 2l0/3), PcrA reaches configuration 4'. In this configuration the system has a high probability of transiting to the p state (this corresponds to ATP hydrolysis/dissociation of ADP and Pi), reaching the configuration 1''. Since PcrA has now moved 1-nt distance to the left compared to the original configuration 1', we denote the new configuration by 1''.
In the weak coupling scenario, the barrier-crossing domain movements (governed solely by potential Ui) are assumed rate-limiting. The barrier heights in configurations 1' and 2' assume the values A2p = 24.7 kcal/mol, A1p = 12.8 kcal/mol and in configurations 3' and 4' assume the values A2s = 10.7 kcal/mol, A1s = 19.4 kcal/mol (see Eq. 14). The values were chosen through their -dependence in Eq. 14 ( = –0.62 kcal/mol) such that the simulated speed of translocation agrees with the observed speed (13).
Strong coupling scenario
The random motion of PcrA along ssDNA in the strong coupling scenario is presented in Fig. 5 b. The left panel of Fig. 5 b presents the corresponding stochastic motions of domains 2A and 1A. The motions are qualitatively similar to those in the weak coupling case. Over the period of 10 hydrolysis cycles the positions of domains 2A and 1A change by 65 Å, reflecting 10 translocation steps. One can recognize again that the domains move in a nearly synchronous fashion toward the –x-direction (see movie provided in Supplementary Material). The scenario underlying this motion is shown in the right panel of Fig. 5 b. The system traverses again sequentially configurations 1' 2' 3' 4' 1'' as identified in the figure. The configurations are also characterized through state labels pp, ps, ss, and sp introduced earlier. In configuration 1', PcrA is in the pp state, domains 2A and 1A are separated by 1-nt distance (lp in Vp, see Eq. 4), and are both essentially "stuck" in potential Uip exhibiting high barriers. Arrival of ATP for binding leads to configuration 2' and the intermediate ps state. In this state, domains 2A and 1A experience the interaction potential Vs (corresponding to a harmonic spring, with ls lp) while individual potentials Uip maintain the pp state form. As a result, 1A experiences a lower barrier (combination of U1p and Vs) and can move readily; 2A also experiences a reduced barrier, but a higher one than 1A and still remains stuck. As Vs quickly draws the domains 2A and 1A close to each other (through the motion of 1A, but not 2A), PcrA relaxes to the 3' configuration that corresponds to the ss state. In this configuration, both domains are stuck by Uis exhibiting high barriers. In the 3' configuration, ATP hydrolysis can take place. Once this happens and the hydrolysis products ADP and Pi start dissociating, configuration 4', corresponding to the sp state, is reached. In configuration 4' the interaction potential changes from Vs back to Vp while individual potentials Uis retain the prior form. As a result, domain 2A experiences a lower barrier (combination of U2s and Vp) and can move readily; 1A also experiences a reduced barrier, but a higher one than experienced by 2A, and remains stuck. The relaxation of PcrA under Vp moves 2A to the left and brings the system back to configuration 1'. However, since PcrA has now moved 1-nt distance to the left compared to the original configuration 1', we denote the new configuration by 1''.
In the strong coupling scenario, the barrier-crossing domain movements (under potential Ui and V') are not rate-limiting. In the present case, the transitions 1' 2' and 3' 4' corresponding to the arrival of ATP and occurrence of ATP hydrolysis are assumed to be much slower than the domain motions and, therefore, are rate-limiting for the overall translocation process. This provides us with more leeway in the choice of the barriers Ai. In fact, -values in the range between –1.2 and 0.4 kcal/mol give reasonable Ai values such that translocation occurs with a speed consistent with observation (13). We used = 0 kcal/mol, which gives barrier heights A2s = 13.2 kcal/mol, A1s = 21.9 kcal/mol, A2p = 27.2 kcal/mol, and A1p = 15.3 kcal/mol (see Eq. 14).
The intrinsic mechanism of unidirectional translocation
Fig. 5 demonstrates that both suggested PcrA translocation mechanisms, involving weak coupling and strong coupling, lead to unidirectional motion along ssDNA. This behavior is described through the solution of coupled Langevin equations. Here we seek an alternative description in terms of rate equations that provide a more systematic explanation for the unidirectional PcrA translocation as well as a convenient mathematical formulation to represent the associated translocation velocity. Our approach is closely related to the generic description developed in Betterton and Jülicher (23). We account for the motion of PcrA in terms of the average position x of its two domains 1A and 2A, defining x = (x1 + x2)/2. The description needs to also attribute to the moving PcrA the various states and configurations underlying the translocation process. These states and configurations are defined for both scenarios in Fig. 5 (right panels). The linking between x and states/configurations is presented in the schematic diagram in Fig. 6, again separately for the two scenarios. We begin with Fig. 6 a, which shows the weak coupling scenario (see Fig. 5 a). The figure places the states/configurations along the position axis. The states correspond to discrete values of position x, namely 
. For the sake of better presentation, the states/configurations are arranged in two tiers, I and II. The correspondence with the configurations 1', 2', 3', 4', 1'' in Fig. 5 a is indicated in Fig. 6 a: 1' corresponds to point (j, I), 2' to point (
, I), 3' to point (, II), 4' to point (j – 1, II), and 1'' corresponds to point (j – 1, I). The configurations capture the cyclic translocation dynamics of PcrA and repeat themselves along the position axis. For example, configuration 1'', located at (j – 1, I), is shifted by one unit along the position axis from configuration 1'. The translocating PcrA cycles through the states/configurations by undergoing transitions between them. The transitions are indicated in Fig. 5 a and include, for example, 
, 
, 
, and 
. The transitions will be described mathematically by linear rate equations and associated rate constants. The rate constants r1, r2 (s1, s2) correspond to the forward (backward) motion of domains 1A and 2A, respectively (we note that the forward direction in Fig. 6 is to the left); transitions starting from tier I and II to the other are described by rate constants I, II; reverse transitions are denoted by primes, e.g., r'1 and 'II. In the weak coupling scenario the transitions denoted by I, II, i.e., transitions that do not translocate PcrA, are fast relative to the transitions that do translocate PcrA; this is also indicated in Fig. 6 a. The latter transitions are governed by the potentials Ui. The corresponding energy barriers are drawn in Fig. 6 a to indicate which transitions are feasible (only those connected with low barriers).
As stated in Table 1, the average time for the transition 1' 2', which corresponds to domain 1A moving forward and is governed by potential U1p, i.e., 1/r1, is estimated in Appendix B to be 19 ms; the average time for the transition 1' 2, which corresponds to domain 2A moving backward, i.e., 1/s2, is estimated likewise to be 106 s. We note that 1A, in principle, can move both forward and backward facing the same low barrier height; however, due to the geometric constraint (x1 – x2 l0, see Methods), 2A and 1A cannot be separated too far and since the initial separation between 2A and 1A in configuration 1' is large (x1 – x2 l0), only the forward motion of 1A is allowed. Values for transition rates (times), 1/r2 and 1/s1, are also stated in Table 1. One can recognize that 1/r1 + 1/r2 assume a value of 20 ms as we fit it to the observed speed (13). The values r'i and s'i (i = 1, 2) are not listed and they assume the same values (in the weak coupling scenario) as ri and si, respectively. The transition rates between tiers I and II are unknown in this case; however, since the transitions denoted by I and II are relatively fast, i.e., 1/i << 20 ms (i = I, II), and should be much faster than the reverse transitions (primed), we concluded i >> 'i and used 1/i = 100 ns and 'i = 0, accordingly.
|
, I), and 4' to point (, II). However, in this case the domains experience an interaction described by V'(x1, x2), which leads to an important variation. After transition 1' 2', corresponding to the arrival of ATP, domains 2A and 1A are still separated (by lp = l0) while the potential Vs(x1, x2) with a short equilibrium length (ls = l0/3) is switched on, placing the system in an energized state. This is depicted in Fig. 6 b by placing configuration 2' upwards on the energy profile by E k(lp – ls)2/2. Likewise, configuration 4', corresponding to PcrA with ATP just hydrolyzed, is depicted in an upward position on the energy axis again by E, where domains 2A and 1A are close while the potential Vp(x1, x2) with a large equilibrium length (lp = l0) is switched on. Therefore, the transitions 2' 3' and 4' 1'', which correspond, respectively, to domain 1A and 2A moving forward, are energetically favorable and happen fast.
The associated rate constants (times) for the transitions involving domain movements in Fig. 6 b are also provided in Table 1, with the values estimated in Appendix B. The primed rate constants, not shown in the table, are estimated as in Betterton and Jülicher (23): r'i = ri exp(–E/kBT) and s'i = si exp(–E/kBT), thus r'i and s'i (i = 1, 2) assume much smaller values than ri and si. For transitions between tiers I and II (not involving domain movements), 'I = I is used assuming an isoenergetic transition (55). We use 'II = I for simplicity; 'II and II are related by II = 'II exp(µ – 2E/kBT) due to the free energy changes, where µ is the chemical energy generated in one ATP hydrolysis cycle (with µ = 20 kBT used). Since now the transitions between tiers I and II are rate-limiting, 1/I + 1/II has to assume a value of 20 ms according to the observed translocation speed (13).
The transitions shown in Fig. 6 capture the coupling of ATP binding and hydrolysis/dissociation to the translocation of PcrA. The transitions can be cast into a rate equation that describes the probabilities of finding PcrA in the configurations 
shown in Fig. 6 . For example, in the strong coupling scenario, the rate equation of the associated probabilities P1', P2', etc., reads (
),
 | (15) |
with µ, = 1', 2', etc., are probability fluxes defined through
 | (16) |
= – Jµ. Here 
are rate coefficients that can be identified from Fig. 6. For example, in the weak coupling scenario 
and 
(= 0).
The rate equation can be cast into the form
 | (17) |
TOP
ABSTRACT
INTRODUCTION
