J. Phys. Chem. B 2001, 105, 12601-12608
12601
Brownian Dynamics Simulating the Ionic-Strength Dependence of the Nonspecific Association of 434 Cro Repressor Binding B-DNA Fan Yang,† Igor V. Ouporov,† Coretta Fernandes,§ Dagmara Motriuk,‡ and Kathryn A. Thomasson*,† Department of Chemistry, UniVersity of North Dakota, Grand Forks, ND 58202-9024, Lansing Community College, P.O. Box 40010, Lansing, Michigan 48901-7210, and DiVision of Physical Science, Casper College, Casper, Wyoming 82601 ReceiVed: June 5, 2001; In Final Form: August 29, 2001
The spatial and temporal interplay of the association, dissociation, and translocation of proteins along DNA is of central importance to genome regulation. Brownian dynamics (BD) is employed to study the energetics of nonspecific binding of 434 Cro repressor protein (Cro) to model B-DNA as a function of ionic strength. BD simulates the diffusional dynamics as Cro encounters the DNA surface including the steric effects of encounter between the irregular surfaces of the protein and DNA and the electrostatic effects based on a finite difference numerical solution of the Poisson-Boltzmann (PB) equation. The free energy of a CroDNA encounter is determined by computing the potential of mean force versus the radial distance from the protein center to the DNA helix axis. The PB equation is solved by two approximations, the linearized form (LPB) and the full PB form (FPB); periodic boundary conditions are implemented for both solutions. The effect of the solution of the PB equation on the predicted free energy curve shows that both methods give qualitatively similar results but statistically the best results are achieved using the full PB. Both methods show that the depth of the free energy curve dramatically decreases as ionic strength increases from 0.01 to 0.50 M. For example, the full PB gives average depths of the free energy curves at -7 ( 2 and -0.51 ( 0.09 kcal/mol for ionic strengths of 0.01 and 0.50 M, respectively. The lifetimes of nonspecifically docked states, estimated from the free energy BD profile, depend on the ionic strength and are approximately 1000 and 0.03 µs for ionic strength of 0.01 and 0.50 M, respectively.
Introduction Genome regulation is a critical step in reproduction for all organisms. Protein-DNA interactions play a key role in such processes. Many proteins have evolved to activate or deactivate genes for transcription. The spatial and temporal interplay of the association, dissociation, and translocation of proteins along DNA is of central importance to regulation of genome expression1. One goal is the elucidation of the control of gene expression at the molecular level. A piece of this puzzle involves how the protein binds DNA nonspecifically. Two of the common models proposed for how nonspecific binding contributes to the location of specific sites on the DNA are the sliding model and the intersegment transfer model. The sliding model is a one-dimensional diffusion mechanism in which nonspecifically bound proteins play a major role in facilitating guidance to specific sites on the DNA2 through a reduction-in-dimensionality effect3. The enzymatic functions of nucleic acid polymerases and helicases4 are facilitated by their ability to rapidly slide along DNA. The study of Escherichia coli RNA polymerase on DNA provided direct evidence for sliding as a mechanism for relocation of the enzyme on DNA.5 Photolyase uses sliding to scan DNA for damage sites.6 Other examples in which the sliding mechanism is observed in binding * To whom correspondence should be addressed. E-mail: kthomasson@ chem.und.edu. Fax: (701)777-2331. † University of North Dakota. § Lansing Community College. ‡ Casper College.
DNA include the FLP recombinase of Saccharomyces cereVisiae,7 BAL 31 nuclease from Alteromonas espejiana binding to single-stranded DNA,8 the restriction enzymes endonuclease EcoRI9, EcoR124II,10 and EcoRV,11 and the bacterial restriction enzyme BamHI.12 The intersegment transfer model is the one in which the protein binds nonspecific DNA and searches for a specific site by binding a second strand of DNA before disassociating from the first.13 Analysis of the kinetics of the interaction of lac repressor with lac operator found that sliding was not involved but the increased association rate to the operator was due to intersegment transfer involving lac-operator-like sequences.14 Furthermore, evidence supports that steroid receptors also prefer the intersegment transfer model because their association with DNA is most favorable under high concentrations of DNA.13 Further study of energetic and dynamic aspects of nonspecific binding is necessary, including the determination of the interplay between the dynamics of radial and lateral diffusion (“sliding”) and the mean lifetime and sliding lengths of individual encounters. Such information is essential to understanding the competition between direct diffusion, intersegment transfer, and sliding mechanisms in guidance to the specific site and how this competition is modulated by salt concentration. Brownian dynamics (BD), coupled to rigorous numerical methods for computing electrostatics, provides an excellent opportunity to determine a quantitative understanding of not only the influences of the long-range electrostatic forces between a spatially complicated array of charges on the protein and DNA
10.1021/jp012122z CCC: $20.00 © 2001 American Chemical Society Published on Web 11/22/2001
12602 J. Phys. Chem. B, Vol. 105, No. 50, 2001 but also the orientational aspects required for protein association with DNA. BD can characterize the energetics of the nonspecific interaction of proteins with the electrostatic force field of B-DNA and the way in which the interaction is affected by changing ionic strength. BD is a powerful tool that generates a Boltzmann population of the protein conformational states in proximity of DNA. This in turn provides a potential of mean force for the interaction between the protein and DNA along a radial coordinate defining the nonspecific docking of the two molecules. From the potential of mean force, the lifetimes of nonspecific complexes and the possible sliding lengths of a protein on DNA during a single encounter are estimated as a function of ionic strength. Specifically, BD simulations examine the energetics of nonspecific binding of 434 Cro repressor protein (Cro) to model B-DNA. The nonspecific Cro-DNA interaction has been shown to be salt-sensitive, suggesting that the interaction is predominantly electrostatic.15 Complementary electrostatic interactions between the negatively charged B-DNA and these proteins have been shown already in theoretical studies to substantially contribute to the formation of specific and nonspecific proteinDNA complexes.16-18 BD reproduces the ionic strength dependence of the interactions of Cro with DNA over a range of ionic strengths (0.01-0.50 M). Methods Atomic Models. Cro 434 Models: Crystal Model (X-ray). Atomic scale modeling of the Cro homodimer began with the 2.5 Å resolution crystal structure of Mondragon and Harrison16 found in the RCSB Protein Data Bank19 (accession code 3CRO); this structure is the phage 434 Cro dimer and is complexed with a 20 base-pair piece of DNA. The DNA was removed, and the protein part was saved; henceforth, this model will be referred to as the X-ray model. The crystal structure provided coordinates for the first 66 residues of the monomer designated L and the first 63 residues of the monomer designated R. The missing C-terminal residues were unresolved because of disorder.20 Two models that included these residues were built as follows to determine whether these missing C-terminal residues were critical for the simulations of nonspecific interactions. Cro 434 Models: First Homology Model (Cro1). The 2.35 Å resolution crystal structure of the uncomplexd 434 Cro monomer21 was retrieved from the RCSB Protein Data Bank19 (accession code 2CRO); in this structure, the first 65 residues were ordered enough to provide coordinates so that only five of the C-terminal residues were missing. A structural alignment with the L monomer of 3CRO using Homology in InsightII (MSI, San Diego, CA) showed no major structural differences. A related crystal structure, the Cro insertion mutant (accession code 1ORC), was found in the RCSB Protein Data Bank;19 this protein was engineered to have the same function and a similar helix-turn-helix motif as Cro. The ORC structure by Albrigth, Mosing, and Matthews22 was a monomer and had a resolution of 1.54 Å. Through the use of Homology in InsightII, the sequence of a monomer of 434 Cro was aligned with the sequence of ORC (Figure 1). Then the coordinates for the missing 434 Cro residues were borrowed from the overlapping residues of ORC; the resulting monomer structure was subjected to 1000 steps of steepest descents energy minimization with the DISCOVER consistent valence force field (MSI, San Diego, CA). To create the dimer, the minimized monomer was superimposed onto the first 66 residues of the L monomer of 3CRO and another copy of the minimized monomer was superimposed onto the first 63 residues of the R monomer of
Yang et al.
Figure 1. Sequence alignment of Cro with ORCsthe basis for model Cro1. Each box includes the sequence that Homology identified to be similar. The residues in boldface were used to assigned coordinates to the missing 2CRO residues; the missing 2CRO residues are listed in small letters. All other coordinates assigned to model Cro1 came directly from 2CRO, the monomer X-ray structure. The EnDashhyphens refer to gaps in the sequence alignment identified by Homology.
3CRO. This placed the two monomers in the same orientation as the original 3CRO structure. Finally, the new dimer model was subjected to 1000 steps of steepest descents minimization. Henceforth, this resulting model is referred to as Cro1. Cro 434 Models: Second Homology Model (Cro2). Because residues of ORC that were used to build the model C-terminus for Cro1 did not exactly match the sequence of 434 Cro, a second model was built using a different repressor protein that still showed a helix-turn-helix motif and provided residues that overlapped with the C-terminal residues of 434 Cro. The theoretical model for the lactose repressor (accession code 1LTP) built by Nichols et al.23 was retrieved from the RCSB Protein Data Bank19 and aligned with the monomer sequences of 2CRO and 1ORC (Figure 2). As before, the coordinates for the first 64 residues were obtained from the 2CRO structure. Coordinates for residue Lys-64 were obtained from 1ORC and the coordinates for rest of the C-terminus were obtained from 1LTP. This second monomer model was minimized and transformed into a dimer in the same way as Cro1. This second model is referred to as Cro2. DNA Model. The model of nonoperator B-DNA based on the repeating duplex unit d(CGCGAATTCGCG) was the same as that described previously.24 The model had exactly nine DNA turns. The B-DNA helix was centered on a cubic grid placing the helix axis along the Z axis. Molecular Charges. Charge assignments for the DNA and all three Cro models were computed by the Tanford-Kirkwood method with the surface accessibility modification25 at pH 7, at a temperature of 298.15 K, and at various ionic strengths: 0.01, 0.05, 0.10, 0.20, 0.30, 0.40, and 0.50 M. Partial charges were assigned to best represent the details to be used for determining the electrostatic field. The resulting charges under the different conditions are found in Table 1. Formal charges were also assigned to the Cro models for use in the BD simulations; these formal charges were found on the surface of the protein and fully exposed to water. Electrostatic Field. Once charges were assigned, the electrostatic field about the DNA or the Cro was found by solving the Poisson-Boltzmann equation as implemented in the program MacroDox.26 For DNA, both the full Poisson-Boltzmann (FPB) and the linearized Poisson-Boltzmann (LPB) were solved with periodic boundary conditions; for details, see Thomasson et al., 1997.24 For the Cro models, only the linearized PoissonBoltzmann was solved to aid in visualizing the effect of adding the C-terminus. During BD simulations, the formal charges on each of the Cro models were used as test charges imbedded in
Brownian Dynamics
J. Phys. Chem. B, Vol. 105, No. 50, 2001 12603 coordinate was monitored throughout each trajectory of 107 steps. For each Cro model at each ionic strength, a series of 10 single-trajectory simulations were performed and averaged. The potential of mean force (PMF) was determined from the distribution F(Rc) of the residence times in radial concentric bands of 1 Å thickness by the statistical mechanical formula
PMF ) -kBT ln(F(Rc))
(1)
The free energy, A, was calculated from the PMF by
A ) -kBT ln(F(Rc)) + C ) PMF + C
Figure 2. Sequence alignment of Cro with LTP and ORCsthe basis for model Cro2. Each box includes the sequence that Homology identified to be similar. The residues in boldface were used to assigned coordinates to the missing 2CRO residues; the missing 2CRO residues are listed in small letters. All other coordinates assigned to the Cro2 model came from 2CRO, the monomer X-ray structure. The hyphens refer to the gaps in the sequence alignment identified by Homology. Note that the LTP sequence shown starts with residue Gly-170 because the preceding residues did not align with any of the other sequences. The ‚‚‚ represents the LTP segment from Pro-188 to Met-232 that did not align with the other sequences.
TABLE 1: Total Charges of Models at pH 7 and Temperature 298.15 K with Varying Ionic Strengtha 0.01
0.05
0.1
0.2
B-DNA -177.8 -179.9 -180.0 -180.0 X-ray +12.6 +12.8 +12.9 +13.0 Cro1 +16.2 +16.6 +16.8 +16.8 Cro2 +16.2 +16.6 +16.8 +16.8
where the constant C was chosen to define where the free energy A returns to zero, i.e., beyond Rc > 53 Å. This value was chosen because it was clear that the DNA electrostatic forces had dissipated by this point and all simulations regardless of the method used to solve the electrostatic potential or the ionic strength converged by this point. The actual value of C was determined by averaging the potential of mean force from Rc between 58 and 65 Å to minimize the noise in the data. By this distance, not only had the potential flattened to a clear zero point but also the noise in the simulation was minimal. The noise only became an issue in the lower ionic strength simulations because the protein spent more time around the DNA than it spent in the region beyond which the electrostatic forces were zero. Association Times. The radial potential of the mean force profile from the inner reflecting wall to 53 Å was used to estimate the mean time that Cro associates nonspecifically with DNA before escape. The mean first passage time, τ, for a particle to escape from a well was determined by Szabo, Schulten, and Schulten.27
τ)
∫aR dRc[D(Rc)Peq(Rc)]-1[∫RR
cm
dyPeq(y)]2
(3)
where
Peq(r) ) {
ionic strength (M) molecule
(2)
0.3
0.4
0.5
-180.0 +13.0 +16.8 +16.8
-180.0 +13.0 +16.9 +16.9
-180.0 +13.1 +16.9 +16.9
a Charges are given in electron charge units. X-ray refers to the 3CRO crystal structure. Cro1 and Cro2 are the two homology built models.
the high-dielectric medium of water representing the Cro surface charges, which are fully exposed to a water environment. Brownian Dynamics (BD). Spatially constrained BD simulations were performed using the same method as described previously.24 The direct force between DNA and Cro and the torque operating on Cro were determined at 25 ps time steps by immersing the test charges of Cro in the high-dielectric solvent region surrounding the low-dielectric DNA cavity. This approximation prevented the need to iterate the PoissonBoltzmann electrostatic potential at every mutual orientation of two low-dielectric cavities; this was an indispensable approximation for the consideration of the 107 steps in the simulation. This procedure allowed for the rotation and translation of Cro in the electrostatic field of the DNA. Each simulation began with Cro at random orientations with a central Cro atom, (γC Leu-45) 70 Å away from the DNA helix axis. The docking coordinate, Rc, was defined as the distance from the central Cro atom to the DNA helix axis. This docking
∫aR dyyd-1 e-βf(y)}-1rd-1 e-βf(r)
(4)
Peq(Rc) is the equilibrium distribution in the radial coordinate Rc, d is the dimensionality of the diffusion (d ) 2 in this case), f(Rc) is the potential of mean force calculated in the BD simulation, D(Rc) is the effective diffusion constant along the radial coordinate (taken to be spatially uniform), and β ) (kBT)-1. The escape distance, a, was set to the distance at which the potential energy returns to zero (53 Å in this case). The position of the inner reflecting wall, R, was determined by the simulation from the distance of closest approach (between 17 and 20 Å depending on the ionic strength with 18 Å being the most common value). The position of the deepest part of the well, Rcm, was also determined by the simulation (23-28 Å depending on ionic strength and the protein model). Sliding Lengths. On the basis of the estimate of the association time, τ, and the assumption that the protein can freely diffuse laterally (“slide”) along the DNA during an encounter, the distance of lateral sliding was approximated from the diffusion half-width, x2Dτ, assuming Cro follows Stokes law. The number of base pairs, N, covered in an encounter was then estimated by
N)
half-width 3.4
where 3.4 Å is the repeat distance per base pair of DNA.
(5)
12604 J. Phys. Chem. B, Vol. 105, No. 50, 2001
Figure 3. Atom by atom comparison of Cro1 to Cro2. The atom shift is the rms difference between identical atoms on the two different models of Cro. Each atom is identified by its number in the file of the model.
Figure 4. Atom by atom comparison of the X-ray structure to Cro1. The setup is the same as Figure 3 but comparing the X-ray and Cro1 models.
Yang et al.
Figure 5. Electrostatic potential around the Cro1 model. Red represents an isopotential surface where a +1e charge possesses the electrostatic energy equivalent to -1 kcal/mol. The black lines are a heavy atom trace of the protein. The blue surfaces are for the electrostatic energy of +1 kcal/mol. Note that there are no holes in the electrostatic surface because the C-termini are included.
Figure 6. Electrostatic potential around the X-ray model. Colors and lines have the same meaning as in Figure 5 but for the X-ray model. Note the wide hole in the potential where the C-terminus is missing.
Visualization of Dynamics. The dynamics of the Cro searching the DNA surface was followed by the newly developed BD helical algorithm.28 When using this algorithm, instead of choosing atom C as the monitor atom against the DNA helix axis, the center of mass of Cro was used. Also, when Cro reached the top or bottom of the cylinder, Cro was transformed by a translation distance of 29.7 Å and a rotation angle of -48.6° to place it in a similar configuration with respect to the DNA. As done previously, Cro started at a random position 70 Å away from the DNA helix axis. Single trajectories with a 10 ps time step were tracked at an ionic strength of 0.05 M for 106 steps. Results Structural Comparison of Cro Models. The two homologybuilt models, Cro1 and Cro2, were quite similar to each other; typically, the rms deviations between the atoms were less than 0.25 Å (Figure 3). There were regions that were exceptions where the rms deviations between the atoms were between 1 and 3 Å. The first region of major difference was the C-terminus (Asp-62 to Ala-71) where the structural coordinates were borrowed from different proteins. The second region of major difference included residues 43-47 that were located where the monomers interfaced. Comparing the X-ray structure to the two Cro models showed larger deviations between the atoms (between 1 and 20 Å) with the largest deviation occurring at the C-terminus where the X-ray structure terminated (Figure 4). In observing the electrostatic potential around the different Cro models, the two homology models appeared very similar as a surface that is predominately positive with small negative patches (Figures 5). The electrostatic potential of the X-ray
Figure 7. Radial free energy profile for Cro interacting with DNA predicted by the FPB. The profiles are for the Cro2 model. The full Poisson-Boltzmann was used to solve the electrostatic potential around the DNA. The box contains the region where the constant C was determined.
structure, however, shows clear holes where the missing C-termini would have been located (Figure 6). From a visual standpoint, the electrostatic potential of the X-ray structure is substantially different from the homology-built models because of the missing charges at the C-terminus. The missing charges also clearly reduced the total charge of the X-ray model compared to the two homology-built models, but the total charges of the two homology-built models were identical at all ionic strengths (Table 1). Free Energy Profiles. The free energy curves have the same shapes for both the LPB and FPB methods and for the different Cro models; all gave qualitatively similar results (Figures 7 and 8). These profiles show a deep well in the vicinity of 23-28 Å for all ionic strengths. The physiological ionic strengths (0.050.2 M) show a shallower well around 35-41 Å; this shallower
Brownian Dynamics
J. Phys. Chem. B, Vol. 105, No. 50, 2001 12605 TABLE 3: Free Energy Well Depths and Mean Lifetimes for Model Cro1-DNA Encountersa free ionic energy strength (kcal/mol) (M) 0.01 0.05 0.1 0.2 0.3
Figure 8. Radial free energy profile for Cro interacting with DNA predicted by the LPB. Parameters are the same as for Figure 7 but applying the linearized Poisson-Boltzmann to solve the electrostatic potential around DNA.
TABLE 2: Free Energy Well Depths and Mean Lifetimes for X-ray Cro-DNA Encountersa free ionic energy strength (kcal/mol) (M)
encounter time (µs)
no. of base pairs
Rcmin electrostatic (Å) method
0.4 0.5
-8(2) -7.6(5) -7(2) -7.1(4) -7(2) -4.5(3) -5(1) -2.6(1) -2.8(3) -1.51(4) -1.8(1) -1.00(6) -1.2(1) -0.62(5)
encounter time (µs)
no. of base pairs
2500(630) 2200(550) 1140(75) 1500(100) 1900(540) 1900(540) 240(13) 670(40) 240(70) 670(190) 2.2(1) 64(4) 4(1) 82(16) 0.25(1) 22(1) 0.33(4) 25(3) 0.089(2) 13.1(3) 0.12(1) 15(1) 0.059(4) 11(1) 0.072(6) 12(1) 0.046(4) 9(1)
Rcmin electrostatic (Å) method 24 25 25 25 24 24 25 25 25 26 26 26 26 26
LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB
a The value of the free energy is given for each method of solving the electrostatic potential with periodic boundary conditions (full Poisson-Boltzmann, FPB; linearized Poisson-Boltzmann, LPB). Rcmin is the inner reflecting wall that is the closest that the Cro central atom γC Leu-45 comes to the DNA helix axis. The error in the free energy represents the standard deviation of the 10 different simulations and the error in the constant C. The numbers in parentheses represent this error in the last significant digit. The error is extrapolated to the derived quantities using relative errors.
0.01
-8(2) -7.6(4)
2500(630) 750(40)
2200(550) 1200(60)
23 23
LPB FPB
0.05
-7(2) -6.1(3)
660(190) 27(1)
1100(320) 230(11)
23 23
LPB FPB
0.1
-7(2) -4.2(1)
24(7) 1.52(3)
210(60) 54(1)
23 23
LPB FPB
TABLE 4: Free Energy Well Depths and Mean Lifetimes for Model Cro2-DNA Encountersa
0.2
-4(1) -2.0(1)
1.6(4) 0.17(1)
60(15) 18(1)
24 24
LPB FPB
-2.4(2) -1.15(2)
0.25(2) 0.080(1)
22(2) 12.4(2)
24 24
LPB FPB
free ionic energy strength (kcal/mol) (M)
-0.9(1) -0.36(2)
0.07(1) 0.044(2)
11(1) 9(1)
25 25
LPB FPB
0.3 0.5
a The value of the free energy is given for each method of solving the electrostatic potential with periodic boundary conditions (full Poisson-Boltzmann, FPB; linearized Poisson-Boltzmann, LPB). Rcmin is the inner reflecting wall that is the closest that the Cro central atom γ C Leu-45 comes to the DNA helix axis. The error in the free energy represents the standard deviation of the 10 different simulations and the error in the constant C. The numbers in parentheses represent this error in the last significant digit. The error is extrapolated to the derived quantities using relative errors.
well represents the contribution of an end-on collision of Cro with the DNA helix.24 The major differences between the LPB and FPB solutions occurred in the value of free energy well depth (Tables 2-4) and scatter in the zero region (Figures 7 and 8). For all three models of Cro and for both the LPB and FPB methods, there is a marked ionic strength dependence. The minimum free energy dramatically decreases as ionic strength increases from 0.01 to 0.50 M (Tables 2-4) so that the well becomes shallower as ionic strength increases. The LPB method always produces deeper wells than the FPB method (Tables 2-4), but statistically the best results are achieved using the FPB, especially at lower ionic strength (Figure 7). Sliding times are affected by changing the ionic strength as well; they decrease as ionic strength increases (Tables 2-4). On the basis of the sliding times, Cro associates with B-DNA for longer periods at the lower ionic strength. Visualization of the Dynamics. Figure 9 demonstrates that the Cro throughout 106 steps does manage to circle around the DNA during an entire trajectory. From the Cro starting position, the protein clearly finds its way to the DNA quickly and spends most of its time near the DNA (Figures 10 and 11). Furthermore, once the Cro reaches the DNA, it steps predominately in one
0.01 0.05 0.1 0.2 0.3 0.4 0.5
-8(2) -7(2) -7(2) -6(2) -6(2) -3.6(4) -4.1(4) -2.2(1) -2.4(1) -1.3(1) -1.5(1) -0.78(8) -1.01(7) -0.54(9)
encounter time (µs) 2400(460) 400(110) 700(200) 65(15) 100(30) 0.80(9) 1.6(2) 0.165(8) 0.216(13) 0.069(6) 0.085(6) 0.041(4) 0.052(4) 0.031(6)
no. of Rcmin electrostatic base pairs (Å) method 2100 800 1200 400 400 40 60 18 20 12 13 9 10 8
25 25 25 25 24 25 25 25 25 26 26 28 26 28
LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB LPB FPB
a The value of the free energy is given for each method of solving the electrostatic potential with periodic boundary conditions (full Poisson-Boltzmann, FPB; linearized Poisson-Boltzmann, LPB). Rcmin is the inner reflecting wall that is the closest that the Cro central atom γC Leu-45 comes to the DNA helix axis. The error in the free energy represents the standard deviation of the 10 different simulations and the error in the constant C. The numbers in parentheses represent this error in the last significant digit. The error is extrapolated to the derived quantities using relative errors.
direction along the DNA (Figures 10 and 11) providing very strong evidence for a sliding mechanism and supporting the sliding assumption used to calculate how many base pairs are covered during an encounter. On the basis of the two trajectories observed, the overall direction that Cro traveled depended on its initial position. In Figure 10, the Cro clearly prefers to move toward the left fairly quickly. In Figure 11, the Cro prefers to move predominately in the opposite direction. Discussion BD simulations characterize the energetics of the nonspecific interaction of 434 Cro repressor with B-DNA by using BD as
12606 J. Phys. Chem. B, Vol. 105, No. 50, 2001
Figure 9. Track of a single trajectory of Cro about DNA. The view is down the helix axis. Each + is the center of mass of Cro1 in a single (10 ps) step around the DNA during a single trajectory of 106 steps.
a tool to generate a Boltzmann ensemble of the protein conformational states in the proximity of DNA.24 The free energy of interaction provides a crude estimate of the lifetimes of nonspecific encounters and possible sliding lengths of the Cro on DNA in a single encounter. The three different Cro models explore the influence of minor changes in protein structure on the interaction. The primary results are the ionic strength dependence of the free energy profile of Cro nonspecifically docking with B-DNA along the radial coordinate for various solutions of the PB equation. The free energy profiles are Boltzmann-averaged over a statistically representative set of mutual orientations of the protein and DNA at a range of distances so that the profiles are truly the free energy of docking resulting from nonspecific interaction.24 This provides novel quantitative information about the effect that the energy well experiences during nonspecific binding between Cro and DNA at a range of separations.24 The choice of method for the electrostatic potential had a large effect on the predicted free energy well depth. LPB simulations produce deeper free energy wells than FPB simulations. It is well-known, however, that the LPB overestimates the electrostatic potential for a molecule with a high charge density such as DNA.29 As a result, the estimate of the sliding
Yang et al. length is quite sensitive to the choice of method used to generate the electrostatic potential around DNA. The LPB also produces more noise in the radial free energy curves making the zero determination less precise. This effect becomes worse as ionic strength decreases causing the error in a simulation to increase as ionic strength decreases. The reason for this phenomenon is that when there exists a strong interaction between Cro and DNA, the Cro spends most of its time in the inner region near DNA and very little time in the outer region where zero is defined. Thus, there are fewer steps to provide better quality statistics in the outer region. One way to combat this problem would be to increase the number of steps taken during a single trajectory or to perform a larger number of these long trajectories. Even when doing this, however, the number of steps in the outer region is still small at the lowest ionic strength (0.01 M) resulting in a less precise evaluation of zero. To make a significant improvement in the statistics in the zero region will become computationally very expensive. Considering that the physiological ionic strength is 0.1-0.2 M, where the Cro is more stable,30 this limit is not a serious problem for this study. The disadvantage to increasing the number of steps or the number of trajectories is the dramatic increase in CPU time for the simulation. Regardless, the free energy well depth remains fixed after 107 steps (i.e., the well depth remains essentially the same even when 108 steps are taken, but 108 steps increase the CPU time by a factor of 8). An experimental study of λ Cro repressor, which is very similar to 434 Cro repressor, measured the free energy of nonspecific interaction of -10 kcal/mol at 0.1 M ionic strength.15 Our simulation using the FPB electrostatics predicts -4.1 ( 0.6 kcal/mol. The difference agrees with experiment to the same power of 10, but an exact match is not possible without taking into account effects such as desolvation and the release of counterions neither of which were included in the simulation. Experimental results suggest that nonspecific CroDNA interactions are entropically driven because of the release of counterions and the Coulombic electrostatic forces play a major role in the association.15 The BD simulations do a good job with the electrostatic interactions but assume that the binding is entropically inhibited as two independent rotors dock to form a complex. For these simulations, the counterions were only considered as a mean field treatment, which ignores the ionion correlations. To obtain better entropic results, it will be
Figure 10. Track of a single trajectory that moves predominantly to the left along the DNA. The first 555 steps of a single trajectory beginning with the Cro1 center of mass at the position x ) -69.4, y ) 9.0, and z ) 3.1. Each + is the center of mass of Cro1 in a single (10 ps) step along the DNA. A few representative steps are labeled by their number in the trajectory. Step number 555 was the first encounter with the end of the cylinder about the DNA model so that subsequent steps, which are not shown, were transformed back into the cylinder in a similar orientation on the DNA.
Brownian Dynamics
J. Phys. Chem. B, Vol. 105, No. 50, 2001 12607
Figure 11. Track of a single trajectory that moves predominantly to the right along the DNA. The first 1713 steps of a single trajectory beginning with the Cro1 center of mass at the position x ) -34.9, y ) 60.7, and z ) -0.1. Note that the first three steps are not visible in this orientation, which is identical to the orientation in Figure 10. Each + is the center of mass of Cro1 in a single (10 ps) step along the DNA. A few representative steps are labeled by their number in the trajectory. Step number 1713 was the first encounter with the end of the cylinder about the DNA model so that subsequent steps, which are not shown, were transformed back into the cylinder in a similar orientation on the DNA.
necessary to incorporate the counterions explicitly, which will require significant changes and testing of the existing BD code. Exploring the ionic strength influence on the interaction, all three Cro models agree that the strongest interaction between Cro and DNA occurs at 0.01 M. The interaction decreases as ionic strength the increases. For nonspecific protein-DNA binding, the interaction is salt-sensitive. This conclusion agrees with the result of Takeda et al. that demonstrates that λ Cro repressor binding of DNA nonspecifically drastically decreases as ionic strength increases.31 By the time an ionic strength of 0.07 M is reached, λ Cro is 50% dissociated, and λ Cro stops binding upon reaching an ionic strength of 0.3 M.31 The shallow wells predicted for the ionic strengths of 0.3 M and above indicate that the 434 Cro is not essentially binding appreciably by this point. Furthermore, the number of base pairs that could be covered in a single encounter is fewer than the operator sequences for 434 Cro, which are 14 base pairs each.32 When considering what individual trajectories actually look like, the assumption of sliding can be supported. Following the trajectory displayed in Figure 10, the Cro first encounters the DNA in an end-on fashion (see Figure 12 for an example), which corresponds to the outer well observed in the free energy curves (Figures 7 and 8). As the trajectory in Figure 10 continues, the center of mass of Cro comes closer to the DNA (see Figure 12 for an example) and it continues to move to the left along the DNA; these more intimate interactions correspond to the deeper well observed in the free energy curves. A different initial position was chosen for the trajectory displayed in Figure 11; just as in Figure 10, the Cro finds its way to the DNA quickly, but instead of moving left, it moves across to the other side of the DNA and then proceeds to the right. Both Figures 10 and 11 show that the Cro continues to choose predominately one direction for its search process supporting the assumption that sliding is occurring. They also demonstrate that the Cro quickly finds the DNA and spends most of its time within close proximity to the DNA. The lifetimes calculated from free energy well depths are very sensitive to the ionic strength. When the encounter time is longer, there is a greater possibility for the protein to slide farther along the DNA before detaching. At the lowest ionic strength (0.01 M), the protein has the potential to cover an average of
Figure 12. Two encounter complexes of Cro with nonoperator B-DNA at different orientations. All molecules are represented as trace ribbons. The complex on the left is representative of complexes with Rc ) 2328 Å, which corresponds to the deepest part of the free energy wells portrayed in Figures 7 and 8. The complex on the right is representative of complexes with Rc ) 35-41 Å, which corresponds to the shallow well in the free energy profiles (Figures 7 and 8). Although the end-on kinds of complexes are not the most stable or intimate possible, they contribute to the region of stability of nonspecific complexes and potentially to the sliding mechanism.
1200 base pairs, which is over 80 operators in length. At the physiological ionic strength of 0.1 M, it has the potential to cover an average of 54 base pairs, which is nearly the length of four operators. This means that Cro is associating long enough to actually search a piece of DNA with a possibility of encountering an operator before dissociating. This is not true for the ionic strengths of 0.3 M or higher. Concerning the protein models and how they individually affect the predicted free energy and sliding times, within the error of the simulations, the three models predict approximately the same results even though the X-ray structure is missing a few charged residues. Although the X-ray structure has a significantly lower charge and a hole in its electrostatic field, the BD simulation does not seem very sensitive to these differences indicating that the missing residues are not critical for nonspecific binding. Finally, BD has successfully reproduced the ionic strength influence of nonspecific Cro-DNA binding. The average free energy predicted using the FPB decreases from -7 ( 2 kcal/ mol at an ionic strength of 0.01 M to -0.51 ( 0.09 kcal/mol at
12608 J. Phys. Chem. B, Vol. 105, No. 50, 2001 an ionic strength of 0.5 M. The free energy dependence on ionic strength results in a strong dependence in the average lifetimes of nonspecifically docked states; when they are averaged over all three Cro models, the lifetimes are 1000 and 0.03 µs for ionic strengths of 0.01 and 0.50 M, respectively. Future studies include further examination of the dynamics of single trajectories, incorporation of counterions to correctly predict the entropic effects, and examination of the temperature, pH, and sequence dependence. Acknowledgment. This research was funded by Grant CHE9322231 from the National Science Foundation and grants from North Dakota EPSCoR. We would like to thank Scott H. Northrup for the MacroDox package and his advice in using it. We also thank Martin Gruebele for his suggestions and advice. References and Notes (1) Lohman, T. M. CRC ReV. Biochem. 1986, 19, 191. (2) Mazur, S. J.; Record, M. T., Jr. Biopolymers 1989, 28, 929. (3) Berg. O. G.; von Hippel, P. H. Annu. ReV. Biophys. Biophys. Chem. 1985, 14, 131. (4) Geider, K.; Hoffmann-Berling, H. Annu. ReV. Biochem. 1981, 50, 233. (5) Kabata, H.; Kurosawa, O.; Arai, I.; Washizu, M.; Margarson, S.; Glass, R.; Shimamoto, N. Science 1993, 262, 1561. (6) van Noort, S. J. T.; van der Werf, K. O.; Eker, A. P. M.; Wyman, C.; de Grooth, B. G.; van Hulst, N. F.; Greve, J. Biophys. J. 1998, 74, 2840. (7) Pan, G.; Sadowski, P. D. J. Biol. Chem. 1993, 268, 22546. (8) Lu, T.; Gray, H. B., Jr. Biochim. Biophys. Acta 1995, 1251, 125. (9) Jeltsch, A.; Alves, J.; Wolfes, H.; Maass, G.; Pingoud, A. Biochemistry 1994, 33, 10215. (10) Ramsden, J. J.; Dreier, J. Biochemistry 1996, 35, 3746. (11) Jeltsch, A.; Pingoud, A. Biochemistry 1998, 37, 2160. (12) Viadiu, H.; Aggarwal, A. K. Mol. Cell 2000, 5, 889. (13) Lieberman, B. A.; Nordeen, S. K. J. Biol. Chem. 1997, 272, 1061. (14) Fickert, R. Muller-Hill, B. J. Mol. Biol. 1992, 226, 59.
Yang et al. (15) Takeda, Y.; Ross, P. D.; Mudd, C. P. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 8180. (16) Warwicker, J.; Engleman, B. P.; Steitz, T. A. Proteins 1987 2, 283. (17) (a) Matthew, J. B.; Ohlendorf, D. H. J. Biol. Chem. 1985, 260, 5860. (b) Zacharias, M.; Luty, B. A.; Davis, M. E.; McCammon, J. A. Biophys. J. 1992, 63, 1280. (c) Jayaram, B.; diCapua, F. M.; Beveridge, D. L. J. Am. Chem. Soc. 1991, 113, 5211. (18) Salemme, F. Proc. Natl. Acad. Sci. U.S.A. 1982, 79, 5263-5267. (19) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Sindyalov, I. N.; Bourne, P. E. Nucleic Acids Res. 2000, 28, 235. RCSB PDB http://www.rcsb.org/pdb/ (accessed 6/25/97). (20) Mondragon, A.; Harrison, S. C. J. Mol. Biol. 1991, 219, 321. (21) Mondragon, A.; Wolberger, C.; Harrison, S. C. J. Mol. Biol. 1989, 205, 179. (22) Albright, R. A.; Mossing, M. C.; Matthews, B. W. Biochemistry 1996, 35, 735. (23) Nichols, J. C.; Vyas, N. K.; Quiocho, F. A.; Matthews, K. S. J. Biol. Chem. 1993, 268, 17602. (24) Thomasson, K. A.; Baumgartner, T.; Czlapinski, J.; Kaldor, T.; Northrup, S. H. J. Phys. Chem. B 1997, 101, 9127. (25) (a) Tanford, C.; Kirkwood, J. G. J. Am. Chem. Soc. 1957, 79, 5333. (b) Tanford, C. Roxby, R. Biochemistry 1972, 11, 2192. (c) Shire, S. J.; Hanania, G. I. H.; Gurd, F. R. N. Biochemistry 1974, 13, 2967. (d) Matthew, J. B. Annu. ReV. Biophys. Biophys. Chem. 1985, 14, 387. (26) (a) Northrup, S. H.; Thomasson, K. A.; Miller, C. M.; Barker, P. D.; Eltis, L. D.; Guillemette, J. G.; Inglis, S. C.; Mauk, A. G. Biochemistry 1993, 32, 6613. (b) Andrew, S. M.; Thomasson, K. A.; Northrup, S. H. J. Am. Chem. Soc. 1993, 115, 5516. (c) Northrup, S. H.; Laughner, T.; Stevenson, G. MacroDox macromolecular simulation program; Tennessee Technological University, Department of Chemistry: Cookeville, TN, 1997. (27) Szabo, A.; Schulten, K.; Schulten, A. J. Chem. Phys. 1980, 72, 4350. (28) Ouporov, I. V.; Knull, H. R.; Thomasson, K. A. Biophys. J. 1999, 76, 17. (29) Gilson, M. K.; Sharp, K. A.; Honig, B. H. J. Comput. Chem. 1988, 3, 327. (30) Brennan, R. G.; Roderick, S. L.; Takeda, Y.; Matthews, B. W. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 8165. (31) Takeda, Y.; Kim, J. G.; Caday, C. G.; Steers, E. J., Jr.; Ohlendorf, D. H.; Anderson, W. F.; Matthews, B. W. J. Biol. Chem. 1986, 251, 8608. (32) Donner, A. L.; Paa, K.; Koudelka, G. B. J. Mol. Biol. 1998, 283, 931.
