Hamiltonian Replica Exchange Method Studies of a Leucine Zipper

Jun 19, 2009 - Address: Dept. of Chemistry, Michigan State University, E. Lansing, MI 48224-1322. Cite this:J. Phys. Chem. ... Robert I. Cukier. The J...
0 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2009, 113, 9595–9605

9595

Hamiltonian Replica Exchange Method Studies of a Leucine Zipper Dimer Li Su and Robert I. Cukier* Department of Chemistry and the QuantitatiVe Biology Modeling InitiatiVe Michigan State UniVersity, East Lansing, Michigan 48824 ReceiVed: January 12, 2009; ReVised Manuscript ReceiVed: April 8, 2009

A melting curve (stability versus temperature) of a leucine zipper dimer was obtained with the use of explicit solvent molecular dynamics (MD). Leucine zippers form stable dimers in solution at physiological temperatures by formation of coiled coil alpha helical structures, using a characteristic heptadic repeating unit. We simulated a 31-residue truncation of the 33 residue parallel two-stranded leucine zipper (GCN4-p1) of the yeast transcriptional activator GCN4 that has four such repeats. The dimer remains bound using conventional MD sampling at the lowest temperature used, while mainly breaking apart at the highest temperature used. However, to obtain a melting curve, the sampling efficiency must be improved, especially at the lower temperatures. Configurational sampling was enhanced with the introduction of a Hamiltonian replica exchange method (HREM), which will be referred to as HTREM, which scales the Hamiltonian in both potential and kinetic energies. The potential scaling is carried out only for the protein-protein and protein-solvent interactions and the kinetic scaling only for the protein degrees of freedom. By limiting the number of scaled degrees of freedom, a smaller number of systems can be used relative to temperature REM where all degrees of freedom are scaled. The HTREM does enhance the sampling especially at the lower temperatures and a melting curve is constructed. A mutant leucine zipper is also simulated where five glutamates are protonated and five lysines are deprotonated in order to investigate the role that electrostatic interactions and salt bridges play in dimer stability. The mutant is considerably less stable than the wild type based on the melting curves. The connection between dimer stability, monomer R helix unwinding, and salt bridge presence is investigated. Among the possible salt bridges of GCN4-p1, those found in experiments are also found in the simulation at the lowest simulation temperature, and the corresponding salt bridge fractions at higher temperature are much lower. A difference between the N terminal and C terminal halves of the monomers regarding their R helix stability is found, with the N terminal less stable than the C terminal parts, in accord with biophysical analyses of two monomeric 16-residue peptides derived from GCN4-p1. 1. Introduction Leucine zippers mediate the dimerization of proteins that regulate cellular activity, such as eukaryotic transcription factors.1,2 Generically, n repeating heptads of residues (abcdefg)n that dimerize (and oligomerize) are examples of coiled coils and form an important class of structural motifs.3-7 Crick8 proposed the “knobs-into-holes” mode of packing to form coiled coils. While packing is clearly essential to coiled coil formation, electrostatic effects also play an important role in their stability.9,10,4,11,12,7 The leucine zipper (GCN4-p1) comprising the dimerization domain of the DNA-binding protein GCN4 has attracted attention as an example of an alpha helical, parallel, double-stranded, coiled coil. The folding pathway and stability of GCN4-p1 and various mutants have been studied by methods such as calorimetry, circular dichroism, kinetics, hydrogenexchange, and NMR.13,11,12,14-17 Coiled coil formation introduces two intertwined concepts; that of the folding of each monomer and monomer-monomer association to form the dimer interface. Thus, coiled coil formation and stability has aspects of both protein folding and protein-protein interactions. From the perspective of molecular dynamics (MD), folding and protein interface formation are challenging to simulate due to barriers in the high-dimensional * To whom correspondence should be addressed. E-mail: cukier@ cem.msu.edu. Tel.: 517-355-9715 x263. Fax: 517-353-1793. Address: Dept. of Chemistry, Michigan State University, E. Lansing, MI 48224-1322.

potential energy surface and the large configuration space that must be traversed. Folding of the monomers and dimerization in GCN4-p1 and closely related structures have been addressed using lattice models,18-21 coarse grained models,22 implicit solvent,23,24 and explicit solvent25-27 simulations. One approach to obtaining information about the nature of coiled coils and the determinants of their stability is by the construction of melting curves. In principal, that can be accomplished by starting simulations from a dimer obtained from a crystal or NMR structure. Simulation of the melting of a dimer and comparison of melting curves for wild type and mutant forms should provide the desired insight. However, even this more limited task presents difficulties for simulation because the interactions responsible for the dimer stability, especially at lower temperatures, will most likely prevent dimer separation. In this work, we introduce an explicit solvent MD-based variant of the temperature replica exchange method (TREM)28-34 that is designed to accelerate sampling in the context of the determination of a melting curve. In TREM, multiple copies of systems are run independently by MD (or Monte Carlo) with different temperatures and, occasionally, two neighbors (different temperature systems) may be exchanged according to the Metropolis-Hastings algorithm.35 Sampling is improved for the lower temperature systems by higher temperature systems overcoming barriers in the potential energy surface.36 By construction, correct Boltzmann sampling is maintained at each temperature. TREM suffers from the deficiency that the number

10.1021/jp900309q CCC: $40.75  2009 American Chemical Society Published on Web 06/19/2009

9596

J. Phys. Chem. B, Vol. 113, No. 28, 2009

of system copies needed scales with the square root of the number of degrees of freedom of the system of interest.37 That is a serious problem when explicit solvent simulations are used, which should be important to the correct description of the monomer-monomer interface. To address this difficulty with TREM, Fukunishi et al.37 introduced Hamiltonian REM (HREM), in which the potential energy functions in the different systems differ only in a limited set of the total number of degrees of freedom required to characterize the system thereby, in principle, reducing the number of systems needed. We recently formulated HREM by modifying (weakening) the protein-protein and protein-solvent interactions while leaving the solvent-solvent interactions at their normal values and used it to explore the conformational space of the pentapeptide Met-enkephalin.38 We found that this HREM MD was capable of sampling both extended and closed conformers, while normal MD started from either extreme could not sample the other distinct conformer. Here, we reformulate the HREM by also scaling the kinetic energy of the protein (and designate the method as Hamiltonian Temperature REM (HTREM)) thereby introducing an effective temperature that permits construction of a melting curve. A virtue of the method is that the data from all the systems is of use for each effective temperature simulated. The HTREM is used to simulate a 31 residue (in each monomer) version of GCN4-p1 using explicit solvent MD. For purposes of comparison, we also create and simulate an in silico mutant that suppresses a great number of the electrostatic interactions that are in part responsible for the stability of a leucine zipper. The melting curves that are obtained do show that the wild type is considerably more stable than is the mutant. The obtained melting curve is defined in terms of monomermonomer contacts without regard for the alpha helical content of the monomers. Also explored is the fate of the alpha helices along the melting curve to see if there is a distinction between the processes of separation and helix destruction. There is evidence that the two monomeric 16 residue peptides obtained by splitting GCN4-p1 into its C terminal and N terminal parts have different properties.5,12 The C terminal part tends to form a stable helix while the N terminal part is only marginally stable. This observation suggests that the C-terminal part is a “trigger” sequence12 that drives coiled coil formation, and the conformations of these separate N and C terminus monomers has been simulated using MD.25 Thus, it is of interest to examine these aspects of helix unwinding and we find that unwinding is more likely for the N-terminal parts of the monomers for wild type while the mutant does not make this distinction between N and C-terminal behavior. The important role of salt bridges and electrostatic interactions in stabilizing the leucine zipper dimer is addressed by tracking their existence at the lowest and highest temperatures simulated. The existence of salt bridges for low temperature and their absence or great reduction at high temperature shows the important role they play in stabilizing the leucine zipper. 2. Methods 2.1. HTREM. To formulate the HREM with the inclusion of the kinetic energy the Hamiltonian is written as Hi(X,P) ) Ki(P) + Vi(X) where Ki(P) is the kinetic energy and Vi(X) is the potential function for the ith system with phase space coordinates X,P ≡ Γ. Between exchange attempts normal MD is run for each system i characterized by Hi. When system interchanges are to be attempted, detailed balance is enforced:

Su and Cukier

R(Γ, Γ′ f Γ′, Γ) Pi(Γ) Pj(Γ′) ) R(Γ′, Γ f Γ, Γ′) Pi(Γ′) Pj(Γ) (2.1) Here, R(Γ,Γ′fΓ′,Γ) is the acceptance probability (transition probability) that phase space coordinates Γ in the ith system and Γ′ in the jth system before exchange results in configuration Γ′ in the ith system and Γ in the jth system after exchange, and Pi(Γ) is the product of the Boltzmann configurational and Maxwellian distributions at temperature T ) 1/(kBβ) for the ith system. The kinetic energy is given by NP

K ) KP + KS )

pi2P

∑ 2mi

iP)1

NS

+ P

pi2S

∑ 2mi

iS)1

(2.2a) S

in terms of the momenta and masses for the NP protein and NS solvent atoms. The potential energy is decomposed as

Vλ(X) ) λVPP(xPP) + √λVPS(xPS) + VSS(xSS)

(2.2b) Thus, the phase space coordinates have been divided into those that are to be scaled, denoted with subscript “P” for protein, and those that are not scaled, denoted by subscript “S” for solvent. The protein-protein interactions are scaled with λ and the protein solvent with λ. This scheme associates scaled coordinates and associated momenta with the atoms, versus an approach that would consider the scaling of interactions. The atom-based scaling is required for use with Ewald-based methods38 and more natural when the kinetic energy is included in the exchange scheme. jλ ) It is convenient to define an effective temperature T j λ) that is based on the averaged kinetic energies 1/(kBβ

3 〈KP〉βλ ) NP / (kBβλ) ; 2

3 〈KS〉β ) NS / (kBβ) 2

(2.3a) according to

(

) (

)

NS NP 1 1 1 1 ) + ) [1 + (λ - 1)xS] ¯β NP + NS β NP + NS βλ βλ λ (2.3b) where xS ≡ NS/(NP + NS). Along with this kinetic definition, the configurational factor can be thought of as in terms of a j λV j λ ≡ βVλ. j λ, where β potential, V Consider that before exchange we have X,K:i and X′,K′:j and after exchange X′,K′:i and X,K:j. The Metropolis rule for exchange between two systems is

R(Γ, Γ′ f Γ′, Γ) ) min(1, exp(-∆(X, X' f X', X) ∆K(P, P' f P', P))) (2.4) For the potential part,

HREM Studies of a Leucine Zipper Dimer

J. Phys. Chem. B, Vol. 113, No. 28, 2009 9597

∆(X, X' f X', X) ) β[(Vi(X') + Vj(X)) (Vi(X) + Vj(X'))] ) β[(λi - λj)(VPP(x'PP) - VPP(xPP)) + (√λi - √λj)(VPS(x'PS) - VPS(xPS))] (2.5) using the scaling in eq 2.2b. For the kinetic energy, if, as is done for the TREM,39,31 the contribution to the exchange probability is to be set to zero, then, after exchange the rule

K':i f

β¯ λi β¯ λj

K':i;

K:j f

β¯ λj β¯ λi

(2.6)

K:j

can be used. Then,

(

∆K ) β¯ λi

β¯ λj β¯ λi

K' + β¯ λj

β¯ λi β¯ λj

)

K - (β¯ λiK + β¯ λjK') ) 0

(2.7)

The simplest way to, after exchange, reassign the momenta for the P and S atoms is to scale uniformly their respective momenta to obtain the desired temperature in eqs 2.3. This exchange scheme will guarantee that Boltzmann equilibrium will result in the extended ensemble of the product of all the systems’ ensembles for a sufficiently long trajectory. j λ) in eq 2.3, a property j λ ) 1/(kBβ The effective temperature T of the entire system of protein and solvent, is equal to the temperature T ) 1/(kBβ) for λ ) 1 and increasingly higher than this temperature as λ decreases. The effective temperature is always less than the temperature that would result if all degrees of freedom were scaled. Note that the TREM limit is regained by eliminating the solvent, xS ) 0, whereby all degrees of freedom would be scaled in eq 2.3b. For fixed values of λi and j λ /β j λ will increase as the fraction of protein (scaled) λj, the ratio β i j atoms increases, indicating that as the fraction of scaled atoms increases, the temperature separation between systems becomes larger and, other things being equal, will require more systems to produce acceptable overlap of neighboring distributions. To summarize, a HREM scheme can be formulated to satisfy the detailed balance condition in the total energy by using effective temperatures, as defined in eqs 2.3, when exchanges are attempted. Between attempts, MD is run for each system where the temperature is maintained at these effective temperatures. In this way, the scaled systems are simulated at higher temperatures than the normal (λ ) 1) system, gaining the advantage of overcoming barriers from the use of higher temperatures in addition to the weakening of the atom-atom interactions. To the extent that there is a substantial reduction in the number of degrees of freedom scaled, fewer systems will be required to provide adequate overlap between successive neighbor distributions. 2.2. Melting Curve. The HTREM method can only produce a “pseudo”-melting curve because different kinetic temperatures are used only for the protein degrees of freedom. Thus, a definition of a melting curve must account for this feature. A definition can be made as

〈I〉(βλ) )

∫ I(xPP) exp(-βλ[VPP(xPP) + wλ(xPP)]) (2.8a)

Figure 1. Itineration of the low T (black) and high T (red) number of bonds broken (defined by the a′-a and d′-d distances) for part of the HTREM trajectory. The data is typical and only this limited time range is shown for clarity. The more-or-less random walk among the number of bonds broken for both temperatures shows the repeated coming together and separation of the dimer that is a signature of proper sampling of the configuration space.

where I(xPP) is an indicator function for the dimer state that is to be specified and

exp(-βλwλ(xPP)) ) 〈exp(-β√λV(xPS))〉β,S

(2.8b)

introduces a λ-dependent potential of mean force wλ(xPP). The notation in eq 2.8b denotes an average over the solvent (S) degrees of freedom at temperature 1/(kβ) for a fixed protein (P) configuration and generalizes the definition of a potential of mean force to account for the system temperatures, 1/(kBβλ). The definition of the indicator function will be adapted to the current problem based on the use of nine intermonomer distances defined in section 3. We will refer to these distances as “bonds” and track various combinations of broken bonds. The definition of “broken” for a particular bond is made by examining its distance histogram over the trajectory at the lowest temperature used and picking a cutoff distance from the histogram. Several possible choices for an indicator function are (1) the number of broken bonds, irrespective of which particular bonds are broken, averaged over the trajectory at each temperature, (2) the numbers (from none to all) of broken bonds averaged over the trajectory at each temperature, or (3) the same as definition 1 but only using some of the interior bonds, because the ends of the dimer tend to fray. For definition 2, defining a separated dimer by adding the average number of broken bonds over, for example, 0+1+2+3 broken bonds is a less restrictive definition of separation into monomers than demanding that more than three bonds be broken for dimer separation. In essence, once fewer than the total of nine bonds have broken, the dimer is considered to have separated. Definition 3 will tend to, especially for the lower temperatures where few bonds are broken, imply that the dimer is more stable than when definition 1 is used. As relative measures for a given melting curve and for comparison between the wild and mutant, there is no change in conclusions between using definitions 1 and 3. Thus, we mainly use the most straightforward definition 1, but also use 2. 2.3. Molecular Dynamics Simulations. The CUKMODY protein molecular dynamics code, with the GROMOS9640 force field, was modified to incorporate the HTREM, based on our previous replica exchange codes.41,38 The systems are run

9598

J. Phys. Chem. B, Vol. 113, No. 28, 2009

Figure 2. WT melting curves based on the native fraction defined as the average number of broken bonds and 0+1+2+3 defined as averages for the possibility of having no, one, two, or three broken bonds for any snapshot. The bonds refer to the nine possible a′-a and d′-d distances between monomer CA atoms. The 0+1+2+3 allows for the intrinsic fraying of the N terminus of GCN4-p1. The three traces for each dimer stability definition are for the third, fourth, and fifth ns of the HTREM simulation and show that the results have converged.

independently on different nodes of a Linux cluster computer and, when exchanges are attempted, information is passed using the message passing interface technique implemented as MPICH. SHAKE42 is used to constrain bond distances enabling a 2 fs time step and temperature is separately controlled with a Berendsen thermostat43 with relaxation time of 0.2 ps for the solvent and 0.002 ps for the solute degrees of freedom. Tight temperature control for the solute was used to be able to define the different solute temperatures for this small number of solute degrees of freedom. For the evaluation of the electrostatic and the attractive parts of the Lennard-Jones energies and forces, the PME method44 was applied with a direct-space cutoff of 8.98 Å, an Ewald coefficient of 0.45, and a 60 × 60 × 60 reciprocal space grid. The electrostatic and Lennard-Jones interactions are uniformly scaled; therefore, as λ decreases, softer Lennard-Jones and reduced electrostatic interactions are obtained. Eight systems are simulated with λ values: 1, 0.944444, 0.888889, 0.833333, 0.777778, 0.722222, 0.666667, and 0.6111110. The lowest temperature used (λ ) 1) corresponds to T ) 275 K for all the degrees of freedom, protein and solvent, and highest (λ ) 0.6111110) corresponds to T ) 450 K (for the protein degrees of freedom). The simulations were carried out in a box with 59.1851 Å sides, having 6912 SPC waters initially to produce the normal water density, and 535 were removed to accommodate the leucine zipper dimer. The starting GCN4-p1 leucine zipper configuration was obtained from its X-ray structure45 (PDB accession code 2ZTA). The residue charge states were set conventionally for pH 7 with the glutamic acid ionized (-1) and lysine and arginine protonated (+1) and the histidines (delta protonated) neutral. Two Cl anions were added to neutralize the system. The in silico mutant protonates five glutamates and deprotonates the five lysines to not change the charge, and thus the mutant can be simulated without changing the number of counterions. Before the HTREM, conventional MD was run on the wild type and mutant at 275 K for 20 ps with one body restraints and then for another 20 ps without one body restraints to relax the structures. All eight systems were then run at their given temperatures, spanning 275-450 K, for 1 ns without exchanges. Subsequently, HTREM was started from the final configurations of the eight systems, with exchanges attempted

Su and Cukier every 40 MD steps. The average exchange acceptance percentage was about 15% for all the exchanges with the mutant showing a fall off to about 10% for the highest temperature system exchanges. This fall off is most likely due to the great destruction of dimer and monomer structure found for the mutant at high temperature, as discussed in section 3. The first 2 ns of HTREM simulation were considered as equilibration time and subsequent data were analyzed in 1 ns intervals to make sure that the considered properties were stable. Protein coordinates were written out every 1 ps. The details are presented in section 3. It is important to verify that during the trajectory the dimer is repeatedly and, on average, uniformly separating and reforming in order to simulate successfully the melting curve. Much like in the construction of a potential of mean force, the values of the reaction coordinate under consideration must be repeatedly sampled to collect sufficient statistics for reliable predictions. For the dimer separation measure, we will use the nine a′-a and d′-d intermonomer distances. Figure 1 displays the itineration of the low T (black) and high T (red) number of a′-a and d′-d distances that are larger than their cutoffs (see section 3 for details) for part of the HTREM trajectory. The data is typical of the entire trajectory; only this limited time range is shown for clarity. The more-or-less random walk for both temperatures shows the repeated coming together and separation of the dimer, which is a signature of proper sampling of the configuration space. In contrast, for normal MD especially at lower temperatures the bond breaking itineration is poor. Thus, obtaining a reliable melting curve for normal MD would require much longer simulation times. 3. Results and Discussion 3.0. Melting Curves and Dimer Separation. Simulations are carried out on the wild type (WT) dimer with the monomer consisting of the first 31 residues of the 33 residue crystal structure, whose last two residues were not resolved and considered to be in a random-coil configuration.45 The residue order is R[M(a)KQL(d)EDK][V(a)EEL(d)LSK][N(a)YHL(d)ENE][V(a)ARL(d)KKL]V(a)G where the four heptads (abcdefg) are bracketed and the a and d positions are indicated. The typical heptad repeat tends to use aliphatic hydrophobic residues at the a and d sites.7 The next closestin-space residues are the g (of one monomer) and e (of the other monomer’s next heptad repeat), and they often correlate with pairing specificity.4 These typically are charged residues, potentially forming salt bridges. For purposes of comparison, we also created and simulated an in silico mutant, referred to henceforth as MUT, that suppresses a great number of the electrostatic interactions by protonating five of the six glutamate residues (6, 7, 10, 11, and 22) and deprotonating the five lysine residues. An advantage of this mutation is that the total charge has not changed; thus, the same number of counterions could be used to make the system overall neutral as for the WT, as must be the case for MD with an Ewald-method-based evaluation of the electrostatics.38 Also, because the van der Waals interactions are not modified by the protonation/deprotonation, the packing abilities of the residues are the same as those of the WT. Furthermore, this MUT also permitted startup of the simulation from the same heavy atom configuration as for the WT. The nine (five a′-a and four d′-d) intermonomer (R carbon-R carbon) distances are used as a measure of stability. The result of a 1 ns T ) 275 (the lowest temperature used) simulation was used to set the cutoff value for each of the nine

HREM Studies of a Leucine Zipper Dimer

J. Phys. Chem. B, Vol. 113, No. 28, 2009 9599

Figure 3. Histograms for the WT of the number of bonds broken (the a′-a and d′-d distances) for the lowest (275 K) and highest (450 K) temperature systems over 1 ns of trajectory data.

Figure 4. The MUT melting curve based on the native fraction defined as the average number of broken bonds (excluding none broken). The bonds refer to the nine possible a′-a and d′-d distances between monomers CA atoms. The two curves correspond to the third and fourth ns of HTREM simulation.

distances, which we will refer to as bonds. The nine crystal structure/MD cutoff distances (in Å) are 6.31/8.0, 6.34/6.5, 5.47/ 6.0, 6.11/6.5, 5.60/6.0, 6.20/6.5, 5.87/6.0, 6.43 /6.5, 6.19/7.0. The cutoffs were picked from the histograms, some quite Gaussian in shape and others with a large-distance tail, of these distances. The melting curves will shift somewhat based on the cutoff definitions, but some experimentation shows that they give a reliable qualitative indication (basically just vertical shifting) especially for comparing WT and MUT curves. The WT melting curve is displayed in Figure 2. The curves labeled native uses as a definition of dimer separation the average number of broken bonds (distances larger than the above listed cutoffs) at each temperature. The curves labeled as 0+1+2+3 are averages for the possibility of having no, one, two, or three broken bonds for any snapshot. There are three curves to each group corresponding to the third, fourth, and fifth ns of HTREM, after 2 ns of HTREM (preceded by 1 ns of no HTREM). The plots indicate that with a minor deviation at the lowest temperature the simulation is well converged. Examination of histograms of the distribution of none-broken to nine broken bonds at the various temperatures also shows that the results are quite reproducible over the three 1 ns intervals. The native data curves imply that over the 3 ns of data the average number of bonds broken at the lowest temperature is 1.8 and at the highest temperature 4.6. Thus, even at the highest temperature used the dimer is not completely separated, while at

the lowest T there is some helix separation by this definition. Histograms of the number of bonds broken for the lowest and highest temperatures are presented in Figure 3. The helicity of GCN4-p1 is approximately 90% at 273 K13 indicating that even at this temperature there is some fraying of the dimer that is attributable (see section 3.1) to the N terminal end. As noted in section 2.3, the only “true” temperature in the HTREM method is the lowest one used, here T ) 275, because it corresponds to the fully coupled (λ ) 1) system. The lack of complete stability at this temperature is consistent with the experimental result. The MUT melting curve is displayed in Figure 4. The curves shown are for the third and fourth ns after 1 ns of no exchange and 2 ns of HTREM simulation. The data are extremely stable; the second HTREM ns (data not shown) is already quite close to the results from the third and fourth ns HTREM simulations. The fraction bound of the MUT is reduced relative to the WT, the reduction increasing with increasing temperature. The average of bonds broken at the lowest temperature is 3.4 (1.8 for the WT) and at the highest temperature 7.8 (4.6 for the WT), indicating the increased separation of the MUT dimer relative to the WT dimer. Histograms of the bonds broken for the lowest and highest temperatures are displayed in Figure 5. At the highest T, the largest proportion of the number of bonds broken is nine (all) while there are no zero or one bonds broken contributions. As reflected in the melting curves, the MUT stability is reduced at all temperatures relative to the WT. The high T data for the mutant shows that it is essentially separated. These results indicate the importance of electrostatic interactions in dimer stability. As noted above, the mutation, by protonating five glutamates and deprotonating five lysines, does not change the van der Waals volume of the dimer but produces a very large electrostatic perturbation. 3.1. Monomer Unwinding. The nature of the monomer conformations and fluctuations along the dimer separation path is also of interest. Figure 6 displays the rmsf (root-mean-square fluctuations) of the WT phi and psi angles for res 3-29 (of one monomer) and 34-60 (of the other monomer) for the lowest temperature of the normal MD data. The excluded phi/psi angles are residues at the monomer ends that do not participate in persistent alpha helices. This WT lowest temperature data provides a basis for the analysis since it is essentially a stable structure from the perspective of a dimer. The fluctuations in the phi and psi angles are quite limited. The corresponding phi and psi trajectory data for selected residues are shown in Figure 7. The standard R helix values for phi angles of -57 deg and

9600

J. Phys. Chem. B, Vol. 113, No. 28, 2009

Su and Cukier

Figure 5. Histograms for the MUT of the number of bonds broken (the a′-a and d′-d distances) for the lowest (275 K) and highest (450 K) temperature systems over 1 ns of trajectory data.

Figure 6. WT rmsf of the phi and psi angles for res 3-29 (of one monomer) and 34-60 (of the other monomer) for the lowest temperature of the normal MD data.

for psi angles of -47 deg1 are drawn on the figure. The alpha helical structure of the monomers is well maintained. The phi and psi results for the HTREM low T simulation are similar to the displayed data. The HTREM-obtained rmsf for the WT at high temperature is plotted in Figure 8. Clearly, the phi and psi fluctuations are considerably larger than for the low temperature, and some are very large. The data in Figure 8 are truncated for clarity. Examining the trajectory data displayed in Figure 9 shows that some of the psi angles are going out of their standard (-180,+180) range indicating that there are “flips” of residues. Evaluating and plotting dihedral angles in this manner is useful to be able to distinguish between fluctuations that destroy alpha helices and those that preserve them. Note that transitions between different nominal dihedral angle values occur quickly, indicating that they are activated processes. The barriers to transition are dominated by the main chain “1-4” interactions, that is, the repulsive and electrostatic interactions centered on CA atoms i and i + 4. Thus, there is a distinction between unwinding and large fluctuations of R-helical structures. The average and standard deviations for the five psi and phi dihedrals are collected in Table 1A for the low T normal MD (Figures 6 and 7) and in Table 1B for the high T HTREM (Figures 8 and 9) data. They confirm that for low T the alpha helices are well maintained, while for high T some are well maintained with just larger phi and psi angle fluctuations and others are going around. The rmsf plots for the mutant are displayed in Figure 10 for the lowest and highest temperature data. The data are again

truncated to clarify the fluctuation scale for most of the residues. The scale of fluctuations for the MUT at low T is similar to that of the WT low T displayed in Figure 6. While the melting curve shows that the mutant is considerably less stable in terms of dimer stability, the monomers are roughly as stable, indicating that the mutant separates but does not unwind at low T. The high T data shows that there are considerable phi and psi fluctuations including considerable going around for the mutant, though not more so, qualitatively, than for the WT. The stability of the C terminus 16-residue monomer and the marginal stability of the 16 residue N-terminus monomer obtained from GCN41-p1 led to the hypothesis that the C-terminus acts as a trigger sequence to nucleate the formation of a coiled coil.5,12 While our simulations are not directed to the kinetics of monomer folding and their correlation with dimer formation, we can investigate if there is a distinction between the C and N terminal halves from the trajectory data. For low temperature, the dimer fraction is large and, as noted above, the R-helical structure is well maintained. In contrast, Figure 11 shows 10 snapshots from the WT high T HTREM data, displayed in ribbon view. The bottom left dimer is drawn in black for the C and gray for the N terminal halves and the C and N termini indicated on the dimers. The snapshots show that the N terminal half is distinctly less stable than the C terminal half. As discussed above, even at the highest T simulated, the wild type GCN4-p1 breaks, on average, 4.6 bonds; thus, it does not completely separate. The snapshots in Figure 11 show that dimer separation may be correlated with unwinding of the N-terminus. In Figure 12, we display a more detailed view of a particular snapshot that shows the N terminal half unwinding, while the C terminal half is maintained both from the perspective of the dimer not separating and the monomers maintaining their R-helical structure. The N terminal half could display residual structure and to investigate the issue the data was analyzed with the clustering algorithm, g_cluster.46 For each input structure, the algorithm counts the number of neighbors within an rmsd cutoff, takes the structure with the largest number of neighbors, along with all its neighbors, as a cluster, and eliminates all the structures within the cluster from the pool of structures. These steps are repeated for the remaining structures in the pool until no structures are left. The 1000 WT high T HTREM snapshots were first superposed on the C terminal stable halves to provide a common framework. Clustering these 1000 snapshots with this method showed that for a rmsd cutoff of 3 Å there is 1 cluster, with a 2 Å cutoff 2 clusters, and with a 1 Å cutoff >

HREM Studies of a Leucine Zipper Dimer

J. Phys. Chem. B, Vol. 113, No. 28, 2009 9601

Figure 7. WT phi/psi angles versus time of the no-HTREM lowest temperature data for selected residues: 3, 8, 13, 18, 23, and 28. The drawn horizontal lines correspond to the R helix average values of phi (-57 deg) and psi (-47 deg).

Figure 8. WT rmsf of the phi and psi angles for res 3-29 (of one monomer) and 34-60 (of the other monomer) for the highest temperature of the HTREM data. The data is truncated to clarify the fluctuation scale.

500 clusters. This indicates that the N terminal half of reach monomer is sampling a more-or-less random set of conformations, at least on the 1 Å scale of resolution. The mutant is less stable at the lowest temperature (about 60% native) than is the wild type and, as shown in Figure 13, snapshots of the mutant configuration show that the R-helical structure is maintained. Thus, even though there is more dimer separation of the mutant, at this temperature the monomers are essentially unperturbed. For the mutant high-temperature data, the dimer separation is essentially complete and, as is clear from Figure 14, there is a variety of monomer configurations, from maintained alpha helical structure to the alpha helices completely unwound. From the configurations shown in Figure 14 and an examination of the complete data set there does seem to be more stability for this mutant in the N terminus half than in the C terminus half, in contrast with the wild type. 3.2. Electrostatic/Salt-Bridge Interactions. The role of electrostatic interactions, in particular salt-bridge formation, as contributors to the stability of coiled coils has not been well established.11,12,47 The compensation between salt bridge formation and solvation of the participating charged side chains is known to be quite subtle, and the relative free energy magnitude and even the sign of it are under debate. There are also differences in conclusions that are due to differing experimental

measurements, that is, pKa measurements versus double mutant cycle thermodynamic analysis.10 The geometry of coiled coils suggests that salt-bridges can form between the side-chains of residues i and i + 3, i + 4 on the same monomer and between residues i′(g) and i + 5(e) (for parallel coiled coils) on different monomers. In the crystal structure of GCN4-p1, there are intramonomer, Glu22-Arg25; Glu11-Lys15; Lys8-Glu11, and intermonomer, Glu22′-Lys27; Lys15′-Glu20, salt bridges. However, not all these salt-bridges are present in both monomers, and some distances are quite long. NMR NOE data on GCN4-p111 did not find evidence for the Glu11-Lys15 and Lys15-Glu20 salt bridges that are only present in one monomer of the crystal structure. The role of the Glu22-Arg25 intramonomer salt bridges in stabilizing the dimer has been stressed.9,48 There are a number of ways of defining a salt-bridge from simulation data. For example, one could measure the distance between a glutamate CG (the carboxylate group’s carbon) and a Lysine NZ (the amino group’s nitrogen) and consider ∼4.5 Å as a cutoff for a salt bridge. The problem with this type of measure is that if the salt-bridge distance is larger than the cutoff then one does not know if there is some shorter, perhaps relevant distance, between these two side-chains. Thus, we track throughout a trajectory the shortest distance between an atom in the side chain of glutamate/glutamic acid (when protonated as in the mutant) and any of the three (two in the deprotonated mutant) hydrogen atoms in the side chain of lysine or any of the arginine side-chain hydrogen atoms. From the record of atom types and these (shortest) distances, the presence of salt-bridges as well as other side-chain electrostatic interactions can be monitored. Table 2 summarizes the results of this analysis for the three intramonomer and two intermonomer salt bridges for the lowand high-temperature wild type and low-temperature mutant trajectory data. The wild type low-temperature HTREM data is in good agreement with the conclusions based on the crystal structure and NMR NOE data.11 Salt bridges for the Glu11Lys15 intramonomer and Lys15′-Glu20 intermonomer interactions are not present in significant amounts (