J. Mol. Bid. (1992) 227, 1244-1252
Stability of a Model B-Sheet in Water Douglas
J. Tobias?,
Scott F. Sneddonl
and Charles L. Brooks III
Department of Chemistry Carnegie Mellon University Pittsburgh, PA 15213, U.S.A. (,Received 18 February
1992; accepted 29 June 1992)
We have used molecular dynamics simulations to determine the stability in water of a model p-sheet formed by two alanine dipeptide molecules with two intermolecular hydrogen bonds in the closely spaced antiparallel arrangement. In this paper we describe our computations of the binding free energy of the model sheet and a portion of the free energy surface as a function of a reaction co-ordinate for sheet formation, We used the free energy surface to identify stable conformations along the reaction co-ordinate. To determine whether or not the model sheet with two hydrogen bonds is more stable than a single amide hydrogen bond in water, we compared the results of the present calculations to results from our earlier study of linear hydrogen bond formation between two formamide molecules (the formamide “dimer”). The free energy surfaces for the sheet and formamide dimer each have two minima corresponding to locally stable hydrogen-bonded and solvent,-separated configurations. The binding free energies of the model sheet and the formamide dimer are -55 and --034 kcal/mol, respectively. Thus, the model sheet with two hydrogen bonds is quite stable while the simple amide hydrogen bond is only marginally stable. To understand the relative stabilities of the model sheet and formamide dimer in terms of solute-solute and solute-water interactions, we decomposed the free energy differences between hydrogenbonded and solvent-separated conformations into energetic and entropic contributions. The changes in the peptide-peptide energy and the entropy are roughly twice as large for the sheet as they are for the formamide dimer. The magnitude of the peptide-water energy difference for the sheet is less than twice (by about 3.5 kcal/mol) that for the formamide dimer, and this accounts for the stability of the sheet. The presence of the side-chains and/or blocking groups apparently prevents the amide groups in the sheet from being solvated as favorably in the separated arrangement as in the formamide dimer, where the amide groups are completely exposed to the solvent.
Keywords: P-sheets; conformation
stability; thermodynamic
1. Introduction The a-helices and P-sheets are the two types of prevalent, repetitive secondary structures in folded proteins (Richardson, 1981). In both a-helices and b-sheets, all of the backbone NH and CO groups are involved in hydrogen bonds. However, in b-sheets the hydrogen bonds are between one strand and another, in contrast to the situation in a-helices, where the hydrogen bonds are between groups in the same strand. The strands that comprise the t Present address: Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104, U.S.A. i Present address: Central Research, Pfizer, Inc., Groton, CT 06379, U.S.A. 002%2836/92/201244-09
$08.00/0
hydrogen bonding; simulations
secondary
structure;
sheet are extended, with 4, $ values in the /? region (the upper left quadrant) of a Ramachandran-type plot, near - 120”,140” (Richardson $ Richardson, 1989). There are two possible arrangements for the strands in p-sheets, antiparallel and parallel (Richardson & Richardson, 1989). fn the antiparallel arrangement, the strands FLUI in opposite directions, with the interstrand hydrogen bonds perpendicular to the strands; alternating between a closely spaced pair (Fig. l(a)) and a widely spaced pair (Fig. l(b)). In the parallel arrangement, the strands run in the same direction, and the hydrogen bonds are evenly spaced but alternately slant forward and backward with respect to the strand direction (Fig. l(c)). 1244
0 1992 Acadr,mic I’ress Limited
Stability
1245
of a Model P-Sheet
(4
a” 7 \,+fh;& i
m % HR
y
P
/N~~C~#---.
I o
I 0
Figure
7 Abbreviations used: BPTI, bovine pancreatic trypsin inhibitor; DMISO, dimethyl sulfoxide; MD, molecular dynamics; TPT, thermodynamic perturbation theory.
i-I
p&!-Q e
A
1. The 3 possible b-sheet hydrogen bonding arrangements: (a) closely antiparallel; (c) parallel. Arrows indicate the strand directions (N to C).
In addition to being prevalent secondary structures in folded proteins, B-sheets are also found in folding intermediates. For example, the entire P-sheet in ribonuclease A is present in a folding intermediate, and it appears to be formed rapidly and co-operatively at an early stage in folding (Udgaonkar & Baldwin, 1988). Another example is the main j-sheet of bovine pancreatic trypsin inhibitor (BPTI?), which appears to be intact in a substantial population ( > 20%) of early folding intermediates containing a single (30-51) disulfide bond (States et al., 1987). A hairpin loop of the antiparallel b-sheet also forms in a synthetic analog of the 30-51 intermediate from the BPTI folding pathway (Oas & Kim, 1988). Lewis et al. (1971) and Wright et al. (1988) have proposed that certain types of secondary structures, namely reverse turns and helical structures, and hydrophobic clusters are potential folding initiation structures. These proposals are supported by the observations that such structures are detected in aqueous solutions of small peptides derived from native sequences (Wright et al., 1988). Thus far, there is little evidence to back a similar proposal for b-sheets. One problem is that it is difficult to study small peptide models of P-sheets because intermolecular aggregation often occurs (Oas & Kim, 1988). Isolated P-sheet formation has been observed in model polymer systems, for example polylysine and polytyrosine (Hartman et al., 1974; Auer & Patton, 1976). However, it is not clear that these polymer systems are good models for sheet formation in proteins, in part because the timescale for sheet formation in the model systems is many orders of magnitude longer than the timescales for the complete folding reactions of sheet-containing proteins (Oas & Kim, 1988). To date, there are no known small peptide sequences, taken from /l-sheet
\
/C’N’KYN’t H
spaced antiparallel;
0
(b) widely spaced
regions in native proteins, which form stable sheets when isolated from the remainder of the protein, and hence, sheet-stabilizing tertiary interactions. In the only known case where a piece of a sheet forms in a protein fragment, the PclPp model for the 30-61 intermediate of BPTI, the existence of the sheet seems to rely on the proximity of the C-terminal helical region, brought forth by the 30-51 disulfide bond (Oas & Kim, 1988). If small peptide sequences can be found for which sections of isolated P-sheets are stable in solution, then, using the line of reasoning leading to the proposal that reverse turns, helical structures and hydrophobic clusters are likely folding initiation structures, sheets are also possible folding initiation structures. There is some experimental evidence that suggests sheets might be stable in various solvents: alanine dipeptide molecules (the alanine dipeptide is the blocked alanine residue, Ac-AlaNHMe; AC is the N-terminal blocking group: GOGH,, and NHMe is the C-terminal blocking group, NHCH,) appear to associate in antiparallel b-sheet arrangements in the non-polar solvent carbon tetrachloride (Cung et al., 1976, 1977) and in dilute dimethyl sulfoxide (DMSO)/chloroform (Asakura et al., 1979; Asakura, 1981). The observation that stable sheet structures can form in a good hydrogen bond accepting solvent like DMSO implies to us that they might also form in aqueous solution. The results of the present study support this hypothesis. To gain a quantitative understanding of the contribution of P-sheets to the stability of folded proteins and folding intermediates, and to investigate the possible role of small sections of P-sheets as folding initiation structures, we have begun to use molecular dynamics (MD) simulations to determine the stabilities of P-sheets in water. In this preliminary investigation, we consider a simple model p-sheet formed by two alanine dipeptide molecules with two intermolecular hydrogen bonds in the closely spaced antiparallel arrangement (Fig. 2). We arbitrarily made this choice from the three possible P-sheet arrangements (Fig. 1).
1246
D. J. Tobias
et al
hydrogen bonds is more stable than a, single amide hydrogen bond in water, we compared the results of the present calculations to results from our earlier study of linear hydrogen bond formation between two formamide molecules (the formamidle dimer). Furthermore, we decomposed the free energy differences between hydrogen-bonded and solventseparated conformations into energetic and entropic contributions to understand the stabilities of the model sheet and formamide dimes in terms of solute-solute and solute-water interactions.
2. Materials Figure 2. Reaction co-ordinate for formation model P-sheet, T(H,-0,) = r(O,-H,) E r.
of the
Following our previous work on secondary structural stability (Tobias & Brooks, 1991; Tobias et al., 1991a), we have chosen to focus initially on what we call the “intrinsic stability” of the model b-sheet by studying peptides composed solely of alanine residues. By intrinsic stability we mean the stability imparted by the intrinsic factors, those which involve only the backbone, and hence, are not strongly sequence-dependent, e.g. backbone sterics, backbone hydrogen bonds, backbone-solvent interactions and backbone configurational entropy (Tobias et al., 1991a). With regard to these intrinsic factors, alanine is representative of all of the common amino acids containing /? carbon atoms. For the purpose of thermodynamic calculations, we define a reaction co-ordinate for the formation of the model b-sheets as the Hz-O5 and 0,-H, distances; these two distances are constrained to be equal and fixed during our simulations (Fig. 2). This reaction co-ordinate is convenient for thermodynamic calculations, but it may not be the actual path for the intermolecular association process in proteins. In studies of intermolecular association, the reaction co-ordinate, chosen as the distance between atoms in different molecules, is unbounded since it ranges from zero to infinity. In this case we obviously cannot use simulations to determine the free energy surface over the entire range of the reaction co-ordinate. Instead, we determine the surface over a small, int,eresting range, compute the binding free energy (the free energy difference between the associated and infinitely separated species), and use the binding free energy to place the surface relative to the free energy of the unassociated species. In this paper we describe our computations of a portion of the free energy surface as a function of the reaction co-ordinate for formation of the model b-sheet. In addition to calcula.ting the surface, we calculated the binding free energy of the sheet, and we used the free energy surface to identify stable conformations along the reaction co-ordinate. To determine whether or not the model sheet with two
(a) Calculation
and Methods
of the free
energy
profile
To calculate the free energy profiles for P-sheet formation between two alanine dipeptide molecules, we used MD simulations with holonomic distance constraints and thermodynamic perturbation theory (TE’T: Tobias & Brooks, 1987, 1988). In this approach. we perform a simulation with the distances, r(H,-OS) = r(O,-H,) = r (Fig. 2), constrained at ri, and compute the Helmholtz free energy difference between systems where r = ri and r = ri + dri using the following formula: AA = A(r,+d%,)-A(rj) = -p-’ In(exp(-fi[ll(r,fdr,)--
U(ri)]))r,.
(1)
In eqn (l), /I = l/&T, where k, is the Boltzmann constant, T is the temperature, U(r) is the potential energy of a. system where the values of the distances equal r, and the (_ .),i denotes a statistical mechanical average over the canonical ensemble (constant number of atoms, volume and temperature) where T = yi. In the present work, we used eqn (1) to compute 4 free energy differences from a single simulation with dri = 10.125 8. _+025 A (1 a = 0.1 nm). Hence, in each simulation we obtained a piece of a free energy surface spanning a 0.5 A range of r values. We used 7 simulations with ri = l&5, 2.35, ., 485 A to map out the free energy surface over the range I.6 A < T i 51 A. Since at r z 4.5 a there is enough space for a water molecule between the peptides, this range of the reaction co-ordinate is sufficiently large to enable us to study the equilibrium between hydrogen-bonded and solvent-separated species in detail.
(b) Calculation
of the binding
jree energy
We defined the binding free energy of the model P-sheet as the Helmholtz free energy change. AAbind, resulting from the following process: 2[Ac-Ala-NHMe](aq)
-+ [AC-Ala-NHMe],(aq),
where [Ac-Ala-NHMe], denotes the peptides in a hydrogen-bonded, antiparallel b-sheet arrangement, and 2[AcAla-NHMe] denotes the peptides at infinite separation. Following Jorgensen et al. (1988), we computed dAbi, as the difference between AA( 1) and AA(Z): the free energy changes accompanying the annihilation processes (1) and (2), respectively: (1) Ac-Ala-NHMe(aq) --t O(aq); (2) (Ac-Ala-NHMe],(aq) --t Ac-Ala-NHMe(aq). The notation “O(aq)” means that the peptide Imolecule is removed entirely from solution in process (1). To obtain AA(l) and AA(2), we defined the following hybrid potential function (Fleischman & Brooks, 1987): U(A) = (l--1)U,+AU,+
CJ,,,,
(2)
Stability
qf a Model
where Us and U, are the potential energies of the reactant (e.g. [Ac-Ala-NHMe], in process (2)) and product (e.g. Ac-Ala-NHMe in process (2)), respectively, U,,, is the potential energy of the environment, the part of the system whose identity does not change during the process (e.g. the solvent), and 1 is a coupling parameter (1 = 0 and ;I = 1 correspond to the reactant and product systems, respectively). Then we used TPT to compute incremental Helmholtz free energy differences, AA,, for small changes, dl, in the coupling parameter: AA = A(L,+dL)-A(&) = -p-l ln(exp{
-j[U(l,+dl)-
U(L)]})n,.
(3)
In eqn (3), the (. .)ii denotes an average over the canonical ensemble, where 1 = Ri. Finally, we summed the AA, from a series of simulations with the 3Li and dl chosen to span the full range (0 to 1) of I to obtain the free energy difference for the transformation of the reactant to the product. We used dl = kO.1 to compute two A.Aj from a single simulation, and 10 simulations with Izi = 605, 915, 0.95 for each of the free energy differences, AA (1) and ii(2).
(c) Thermodynamic
decompositions free energies
of relative
In order to understand the relative stabilities of the hydrogen-bonded model P-sheet and formamide dimer (e.g. compared to their respective solvent-separated species) in solution, in terms of microscopic interactions, we have decomposed the free energy differences between hydrogen-bonded and solvent-separated species into energetic and entropic contributions arising from differences in peptide-peptide and peptide-water interactions. The Helmholtz free energy difference between 2 conformations (e.g. 2 points on the reaction co-ordinate) may be written (Tobias et al., 1991a): AA = (AU,,>+(AU,,>-TAS,,..,
(4)
where (AU”,) and (AU,,,) are differences in the average peptide-peptide and peptide-water interaction energies, respectively, and AX,, uv is the sum of the differences in the configurational and peptide-water entropies. For the potential energy function we use, the peptide-peptide interaction energy differences contain contributions from bond, angle, torsion and non-bonded (van der Waals and electrostatic) energies, while the peptide-water differences contain only non-bonded contributions. It is straightforward to compute the interaction energies but difficult to calculate the entropy directly from the simulation data. Since AS,, uv is the only missing term in eqn (4), we obtain it indirectly as the difference between AA and the average interaction energies. (d) Backbone
constraints
In order to keep the peptides in extended tions during all of the reported simulations, following harmonic constraint potential:
8 conformawe used the
,!i-Sheet
1247
(e) Details
of the simulations
To prepare a system for a solution simulation, we placed either a single peptide molecule or a pair of and peptides in the center of a box of water molecules, removed the water molecules that overlapped with the peptide atoms. The box dimensions were chosen to yield approximate agreement with the experimentally observed room temperature water density (1.0 g/‘cm3) after the solute volume was subtracted from the box volume. For the thermodynamic surface calculations. the systems contained 222 to 230 water molecules in rectangular boxes with dimensions 31.1032 A x 155515 A x 155516 A. For the binding free energy calculations, the systems contained either 196 ([AC-Ala-NHMeJz + Ac-Ala-NHMe) or 205 (Ac-Ala-NHMe + 0) water molecules in cubic boxes with edge length 18856 A. We also carried out 2 additional simulations with r fixed at 1.97 A and 3.9’7 A in the rectangular boxes containing 224 to 230 water molecules for the purpose of calculating average interaction energies (see below). All of the simulations were performed with periodic boundary conditions. We used the CHARMM (Brooks et al., 1983) PARAM parameters for the peptides and the Y-site TIPSP model of Jorgensen et al. (1983) for water. The non-bonded energies and forces were smoothly truncated at either 775 A (for the binding free energy calculations) or 975 A (for the free energy surface and interaction energy calculations), using a van der Waals’ switching function and an electrostatic shifting function (Brooks et al., 1983), based on atomic centers, and according to the minimum image convention (Allen & Tildesley, 1989). The Verlet (1967) algorithm was used to integrate Newton’s equations of motion with a time-step of 1.5 fs. Each of the simulations consisted of either 6000 steps (for the binding free energy calculations) or 5000 steps (for the free energy surface and interaction energy calculations) of equilibration and 10,000 steps of data collection. The non-bonded interactions were processed using a list-based algorithm (Verlet, 1967), and the lists were updated every 10 steps. The velocities were peroditally reassigned from a Maxwell-Boltzmann distribution to maintain temperatures of approximately 300 K. The SHAKE constraint algorithm (Ryckaert et al., 1977) was used to keep the water molecules rigid, and to maintain rigid N-H bonds in the peptide molecules, the internal coordinate constraint algorithm (Tobias & Brooks, 1988) was used to keep r(O,-H,) and r(H,-0,) fixed, and the constraint potential U* (eqn (5)) was used to keep the peptides extended. All of the remaining degrees of freedom were allowed to fluctuate. The potential energy of each system was stored every time-step for the calculation of the free energy differences, and the co-ordinat’es were saved every 10 time-steps for the calculation of other average structural and energetic properties which are presented and discussed below. The simulations were carried out using the CHARMM program (Brooks et ak., 1983) on the Gray YMP computer at the Pittsburgh Supercomputing Center. The details of the formamide dimer simulations are given by Sneddon et al. (1989).
3. Results and Discussion where the @i and $i (i = 2,5) are the backbone dihedral angles, Cpaand Go are the reference values at the minimum in U*, and K is the force constant. We used 40 = - 139”, 9a = 140”, and K = 50 kcal/mol rad’ (1 cal = 4184 J) in all of the simulations. We comment below on the effects of this constraint potential on properties computed from our simulations.
(a) E8ect.s
of the backbone
constraints
Before we present our results for the thermodynamics of P-sheet formation, we briefly examine the effects of the harmonic backbone dihedral angle constraints on the perturbation free energy differences. The presence of constraints can alter the
1248
D. J. Tobias et al.
Determination
Table 1 of the binding free energy for the model antiparallel
P-sheet in w&r --
PIWCIZB
AA
(I) AC-Ala-NHMe(aq) + O(aq) (2) [Ac-Ala-NHMe],(aq) + AC-Ala-NHMe(aq) Binding: Z[Ac-Ala-NHMe](aq) + [Ac-Bla-NHMe],(aq) Binding: 2[formamide](aq) --f [formamide],
--
58.64&@14 6412&0.16t -5.48&02Ej -@34&@28§
Free energy differences are given in kcal/mol. The uncertainties were computed as error propagated standard deviations of 100 averages of 100 datasets. t Peptides were constrained so that T = 1.97 A (the minimum corresponding to the hydrogen-bonded B-sheet on the free energy surface, see below) during the calculation of AA(2). :: AAbind = AA(l)-AA(2). 5 Sneddon et al. (1989).
probability distribution of an observable quantity from what it would be in an unconstrained system. If the distribution is altered significantly, then properties computed from statistical mechanical averages of the quantity, such as free energy differences, should be corrected to account for the constraints. The correct statistical mechanical average of the quantity 0 is (Valleau & Torrie, 1977): (0)
= 10 exP(BU*)>* (exp(BU*)>*
’
(6)
where U* is the constraint potential and the notation (. .)* denotes an average over the constrained ensemble. In our case, the quantity of interest is 0 = exp( -a[ U(r, +dr,) - U(ri)]} (see eqn (1)) a,nd U* is the harmonic potential constraining the backbone dihedral angles (see eqn (5)). We can see in equation (6) that if 0 and exp(fiU*) are statistically independent, then the average of the product in the numerator factors into the product of the averages, and therefore (0) =
Stability of a Model B-Sheet in Water Douglas
J. Tobias?,
Scott F. Sneddonl
and Charles L. Brooks III
Department of Chemistry Carnegie Mellon University Pittsburgh, PA 15213, U.S.A. (,Received 18 February
1992; accepted 29 June 1992)
We have used molecular dynamics simulations to determine the stability in water of a model p-sheet formed by two alanine dipeptide molecules with two intermolecular hydrogen bonds in the closely spaced antiparallel arrangement. In this paper we describe our computations of the binding free energy of the model sheet and a portion of the free energy surface as a function of a reaction co-ordinate for sheet formation, We used the free energy surface to identify stable conformations along the reaction co-ordinate. To determine whether or not the model sheet with two hydrogen bonds is more stable than a single amide hydrogen bond in water, we compared the results of the present calculations to results from our earlier study of linear hydrogen bond formation between two formamide molecules (the formamide “dimer”). The free energy surfaces for the sheet and formamide dimer each have two minima corresponding to locally stable hydrogen-bonded and solvent,-separated configurations. The binding free energies of the model sheet and the formamide dimer are -55 and --034 kcal/mol, respectively. Thus, the model sheet with two hydrogen bonds is quite stable while the simple amide hydrogen bond is only marginally stable. To understand the relative stabilities of the model sheet and formamide dimer in terms of solute-solute and solute-water interactions, we decomposed the free energy differences between hydrogenbonded and solvent-separated conformations into energetic and entropic contributions. The changes in the peptide-peptide energy and the entropy are roughly twice as large for the sheet as they are for the formamide dimer. The magnitude of the peptide-water energy difference for the sheet is less than twice (by about 3.5 kcal/mol) that for the formamide dimer, and this accounts for the stability of the sheet. The presence of the side-chains and/or blocking groups apparently prevents the amide groups in the sheet from being solvated as favorably in the separated arrangement as in the formamide dimer, where the amide groups are completely exposed to the solvent.
Keywords: P-sheets; conformation
stability; thermodynamic
1. Introduction The a-helices and P-sheets are the two types of prevalent, repetitive secondary structures in folded proteins (Richardson, 1981). In both a-helices and b-sheets, all of the backbone NH and CO groups are involved in hydrogen bonds. However, in b-sheets the hydrogen bonds are between one strand and another, in contrast to the situation in a-helices, where the hydrogen bonds are between groups in the same strand. The strands that comprise the t Present address: Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104, U.S.A. i Present address: Central Research, Pfizer, Inc., Groton, CT 06379, U.S.A. 002%2836/92/201244-09
$08.00/0
hydrogen bonding; simulations
secondary
structure;
sheet are extended, with 4, $ values in the /? region (the upper left quadrant) of a Ramachandran-type plot, near - 120”,140” (Richardson $ Richardson, 1989). There are two possible arrangements for the strands in p-sheets, antiparallel and parallel (Richardson & Richardson, 1989). fn the antiparallel arrangement, the strands FLUI in opposite directions, with the interstrand hydrogen bonds perpendicular to the strands; alternating between a closely spaced pair (Fig. l(a)) and a widely spaced pair (Fig. l(b)). In the parallel arrangement, the strands run in the same direction, and the hydrogen bonds are evenly spaced but alternately slant forward and backward with respect to the strand direction (Fig. l(c)). 1244
0 1992 Acadr,mic I’ress Limited
Stability
1245
of a Model P-Sheet
(4
a” 7 \,+fh;& i
m % HR
y
P
/N~~C~#---.
I o
I 0
Figure
7 Abbreviations used: BPTI, bovine pancreatic trypsin inhibitor; DMISO, dimethyl sulfoxide; MD, molecular dynamics; TPT, thermodynamic perturbation theory.
i-I
p&!-Q e
A
1. The 3 possible b-sheet hydrogen bonding arrangements: (a) closely antiparallel; (c) parallel. Arrows indicate the strand directions (N to C).
In addition to being prevalent secondary structures in folded proteins, B-sheets are also found in folding intermediates. For example, the entire P-sheet in ribonuclease A is present in a folding intermediate, and it appears to be formed rapidly and co-operatively at an early stage in folding (Udgaonkar & Baldwin, 1988). Another example is the main j-sheet of bovine pancreatic trypsin inhibitor (BPTI?), which appears to be intact in a substantial population ( > 20%) of early folding intermediates containing a single (30-51) disulfide bond (States et al., 1987). A hairpin loop of the antiparallel b-sheet also forms in a synthetic analog of the 30-51 intermediate from the BPTI folding pathway (Oas & Kim, 1988). Lewis et al. (1971) and Wright et al. (1988) have proposed that certain types of secondary structures, namely reverse turns and helical structures, and hydrophobic clusters are potential folding initiation structures. These proposals are supported by the observations that such structures are detected in aqueous solutions of small peptides derived from native sequences (Wright et al., 1988). Thus far, there is little evidence to back a similar proposal for b-sheets. One problem is that it is difficult to study small peptide models of P-sheets because intermolecular aggregation often occurs (Oas & Kim, 1988). Isolated P-sheet formation has been observed in model polymer systems, for example polylysine and polytyrosine (Hartman et al., 1974; Auer & Patton, 1976). However, it is not clear that these polymer systems are good models for sheet formation in proteins, in part because the timescale for sheet formation in the model systems is many orders of magnitude longer than the timescales for the complete folding reactions of sheet-containing proteins (Oas & Kim, 1988). To date, there are no known small peptide sequences, taken from /l-sheet
\
/C’N’KYN’t H
spaced antiparallel;
0
(b) widely spaced
regions in native proteins, which form stable sheets when isolated from the remainder of the protein, and hence, sheet-stabilizing tertiary interactions. In the only known case where a piece of a sheet forms in a protein fragment, the PclPp model for the 30-61 intermediate of BPTI, the existence of the sheet seems to rely on the proximity of the C-terminal helical region, brought forth by the 30-51 disulfide bond (Oas & Kim, 1988). If small peptide sequences can be found for which sections of isolated P-sheets are stable in solution, then, using the line of reasoning leading to the proposal that reverse turns, helical structures and hydrophobic clusters are likely folding initiation structures, sheets are also possible folding initiation structures. There is some experimental evidence that suggests sheets might be stable in various solvents: alanine dipeptide molecules (the alanine dipeptide is the blocked alanine residue, Ac-AlaNHMe; AC is the N-terminal blocking group: GOGH,, and NHMe is the C-terminal blocking group, NHCH,) appear to associate in antiparallel b-sheet arrangements in the non-polar solvent carbon tetrachloride (Cung et al., 1976, 1977) and in dilute dimethyl sulfoxide (DMSO)/chloroform (Asakura et al., 1979; Asakura, 1981). The observation that stable sheet structures can form in a good hydrogen bond accepting solvent like DMSO implies to us that they might also form in aqueous solution. The results of the present study support this hypothesis. To gain a quantitative understanding of the contribution of P-sheets to the stability of folded proteins and folding intermediates, and to investigate the possible role of small sections of P-sheets as folding initiation structures, we have begun to use molecular dynamics (MD) simulations to determine the stabilities of P-sheets in water. In this preliminary investigation, we consider a simple model p-sheet formed by two alanine dipeptide molecules with two intermolecular hydrogen bonds in the closely spaced antiparallel arrangement (Fig. 2). We arbitrarily made this choice from the three possible P-sheet arrangements (Fig. 1).
1246
D. J. Tobias
et al
hydrogen bonds is more stable than a, single amide hydrogen bond in water, we compared the results of the present calculations to results from our earlier study of linear hydrogen bond formation between two formamide molecules (the formamidle dimer). Furthermore, we decomposed the free energy differences between hydrogen-bonded and solventseparated conformations into energetic and entropic contributions to understand the stabilities of the model sheet and formamide dimes in terms of solute-solute and solute-water interactions.
2. Materials Figure 2. Reaction co-ordinate for formation model P-sheet, T(H,-0,) = r(O,-H,) E r.
of the
Following our previous work on secondary structural stability (Tobias & Brooks, 1991; Tobias et al., 1991a), we have chosen to focus initially on what we call the “intrinsic stability” of the model b-sheet by studying peptides composed solely of alanine residues. By intrinsic stability we mean the stability imparted by the intrinsic factors, those which involve only the backbone, and hence, are not strongly sequence-dependent, e.g. backbone sterics, backbone hydrogen bonds, backbone-solvent interactions and backbone configurational entropy (Tobias et al., 1991a). With regard to these intrinsic factors, alanine is representative of all of the common amino acids containing /? carbon atoms. For the purpose of thermodynamic calculations, we define a reaction co-ordinate for the formation of the model b-sheets as the Hz-O5 and 0,-H, distances; these two distances are constrained to be equal and fixed during our simulations (Fig. 2). This reaction co-ordinate is convenient for thermodynamic calculations, but it may not be the actual path for the intermolecular association process in proteins. In studies of intermolecular association, the reaction co-ordinate, chosen as the distance between atoms in different molecules, is unbounded since it ranges from zero to infinity. In this case we obviously cannot use simulations to determine the free energy surface over the entire range of the reaction co-ordinate. Instead, we determine the surface over a small, int,eresting range, compute the binding free energy (the free energy difference between the associated and infinitely separated species), and use the binding free energy to place the surface relative to the free energy of the unassociated species. In this paper we describe our computations of a portion of the free energy surface as a function of the reaction co-ordinate for formation of the model b-sheet. In addition to calcula.ting the surface, we calculated the binding free energy of the sheet, and we used the free energy surface to identify stable conformations along the reaction co-ordinate. To determine whether or not the model sheet with two
(a) Calculation
and Methods
of the free
energy
profile
To calculate the free energy profiles for P-sheet formation between two alanine dipeptide molecules, we used MD simulations with holonomic distance constraints and thermodynamic perturbation theory (TE’T: Tobias & Brooks, 1987, 1988). In this approach. we perform a simulation with the distances, r(H,-OS) = r(O,-H,) = r (Fig. 2), constrained at ri, and compute the Helmholtz free energy difference between systems where r = ri and r = ri + dri using the following formula: AA = A(r,+d%,)-A(rj) = -p-’ In(exp(-fi[ll(r,fdr,)--
U(ri)]))r,.
(1)
In eqn (l), /I = l/&T, where k, is the Boltzmann constant, T is the temperature, U(r) is the potential energy of a. system where the values of the distances equal r, and the (_ .),i denotes a statistical mechanical average over the canonical ensemble (constant number of atoms, volume and temperature) where T = yi. In the present work, we used eqn (1) to compute 4 free energy differences from a single simulation with dri = 10.125 8. _+025 A (1 a = 0.1 nm). Hence, in each simulation we obtained a piece of a free energy surface spanning a 0.5 A range of r values. We used 7 simulations with ri = l&5, 2.35, ., 485 A to map out the free energy surface over the range I.6 A < T i 51 A. Since at r z 4.5 a there is enough space for a water molecule between the peptides, this range of the reaction co-ordinate is sufficiently large to enable us to study the equilibrium between hydrogen-bonded and solvent-separated species in detail.
(b) Calculation
of the binding
jree energy
We defined the binding free energy of the model P-sheet as the Helmholtz free energy change. AAbind, resulting from the following process: 2[Ac-Ala-NHMe](aq)
-+ [AC-Ala-NHMe],(aq),
where [Ac-Ala-NHMe], denotes the peptides in a hydrogen-bonded, antiparallel b-sheet arrangement, and 2[AcAla-NHMe] denotes the peptides at infinite separation. Following Jorgensen et al. (1988), we computed dAbi, as the difference between AA( 1) and AA(Z): the free energy changes accompanying the annihilation processes (1) and (2), respectively: (1) Ac-Ala-NHMe(aq) --t O(aq); (2) (Ac-Ala-NHMe],(aq) --t Ac-Ala-NHMe(aq). The notation “O(aq)” means that the peptide Imolecule is removed entirely from solution in process (1). To obtain AA(l) and AA(2), we defined the following hybrid potential function (Fleischman & Brooks, 1987): U(A) = (l--1)U,+AU,+
CJ,,,,
(2)
Stability
qf a Model
where Us and U, are the potential energies of the reactant (e.g. [Ac-Ala-NHMe], in process (2)) and product (e.g. Ac-Ala-NHMe in process (2)), respectively, U,,, is the potential energy of the environment, the part of the system whose identity does not change during the process (e.g. the solvent), and 1 is a coupling parameter (1 = 0 and ;I = 1 correspond to the reactant and product systems, respectively). Then we used TPT to compute incremental Helmholtz free energy differences, AA,, for small changes, dl, in the coupling parameter: AA = A(L,+dL)-A(&) = -p-l ln(exp{
-j[U(l,+dl)-
U(L)]})n,.
(3)
In eqn (3), the (. .)ii denotes an average over the canonical ensemble, where 1 = Ri. Finally, we summed the AA, from a series of simulations with the 3Li and dl chosen to span the full range (0 to 1) of I to obtain the free energy difference for the transformation of the reactant to the product. We used dl = kO.1 to compute two A.Aj from a single simulation, and 10 simulations with Izi = 605, 915, 0.95 for each of the free energy differences, AA (1) and ii(2).
(c) Thermodynamic
decompositions free energies
of relative
In order to understand the relative stabilities of the hydrogen-bonded model P-sheet and formamide dimer (e.g. compared to their respective solvent-separated species) in solution, in terms of microscopic interactions, we have decomposed the free energy differences between hydrogen-bonded and solvent-separated species into energetic and entropic contributions arising from differences in peptide-peptide and peptide-water interactions. The Helmholtz free energy difference between 2 conformations (e.g. 2 points on the reaction co-ordinate) may be written (Tobias et al., 1991a): AA = (AU,,>+(AU,,>-TAS,,..,
(4)
where (AU”,) and (AU,,,) are differences in the average peptide-peptide and peptide-water interaction energies, respectively, and AX,, uv is the sum of the differences in the configurational and peptide-water entropies. For the potential energy function we use, the peptide-peptide interaction energy differences contain contributions from bond, angle, torsion and non-bonded (van der Waals and electrostatic) energies, while the peptide-water differences contain only non-bonded contributions. It is straightforward to compute the interaction energies but difficult to calculate the entropy directly from the simulation data. Since AS,, uv is the only missing term in eqn (4), we obtain it indirectly as the difference between AA and the average interaction energies. (d) Backbone
constraints
In order to keep the peptides in extended tions during all of the reported simulations, following harmonic constraint potential:
8 conformawe used the
,!i-Sheet
1247
(e) Details
of the simulations
To prepare a system for a solution simulation, we placed either a single peptide molecule or a pair of and peptides in the center of a box of water molecules, removed the water molecules that overlapped with the peptide atoms. The box dimensions were chosen to yield approximate agreement with the experimentally observed room temperature water density (1.0 g/‘cm3) after the solute volume was subtracted from the box volume. For the thermodynamic surface calculations. the systems contained 222 to 230 water molecules in rectangular boxes with dimensions 31.1032 A x 155515 A x 155516 A. For the binding free energy calculations, the systems contained either 196 ([AC-Ala-NHMeJz + Ac-Ala-NHMe) or 205 (Ac-Ala-NHMe + 0) water molecules in cubic boxes with edge length 18856 A. We also carried out 2 additional simulations with r fixed at 1.97 A and 3.9’7 A in the rectangular boxes containing 224 to 230 water molecules for the purpose of calculating average interaction energies (see below). All of the simulations were performed with periodic boundary conditions. We used the CHARMM (Brooks et al., 1983) PARAM parameters for the peptides and the Y-site TIPSP model of Jorgensen et al. (1983) for water. The non-bonded energies and forces were smoothly truncated at either 775 A (for the binding free energy calculations) or 975 A (for the free energy surface and interaction energy calculations), using a van der Waals’ switching function and an electrostatic shifting function (Brooks et al., 1983), based on atomic centers, and according to the minimum image convention (Allen & Tildesley, 1989). The Verlet (1967) algorithm was used to integrate Newton’s equations of motion with a time-step of 1.5 fs. Each of the simulations consisted of either 6000 steps (for the binding free energy calculations) or 5000 steps (for the free energy surface and interaction energy calculations) of equilibration and 10,000 steps of data collection. The non-bonded interactions were processed using a list-based algorithm (Verlet, 1967), and the lists were updated every 10 steps. The velocities were peroditally reassigned from a Maxwell-Boltzmann distribution to maintain temperatures of approximately 300 K. The SHAKE constraint algorithm (Ryckaert et al., 1977) was used to keep the water molecules rigid, and to maintain rigid N-H bonds in the peptide molecules, the internal coordinate constraint algorithm (Tobias & Brooks, 1988) was used to keep r(O,-H,) and r(H,-0,) fixed, and the constraint potential U* (eqn (5)) was used to keep the peptides extended. All of the remaining degrees of freedom were allowed to fluctuate. The potential energy of each system was stored every time-step for the calculation of the free energy differences, and the co-ordinat’es were saved every 10 time-steps for the calculation of other average structural and energetic properties which are presented and discussed below. The simulations were carried out using the CHARMM program (Brooks et ak., 1983) on the Gray YMP computer at the Pittsburgh Supercomputing Center. The details of the formamide dimer simulations are given by Sneddon et al. (1989).
3. Results and Discussion where the @i and $i (i = 2,5) are the backbone dihedral angles, Cpaand Go are the reference values at the minimum in U*, and K is the force constant. We used 40 = - 139”, 9a = 140”, and K = 50 kcal/mol rad’ (1 cal = 4184 J) in all of the simulations. We comment below on the effects of this constraint potential on properties computed from our simulations.
(a) E8ect.s
of the backbone
constraints
Before we present our results for the thermodynamics of P-sheet formation, we briefly examine the effects of the harmonic backbone dihedral angle constraints on the perturbation free energy differences. The presence of constraints can alter the
1248
D. J. Tobias et al.
Determination
Table 1 of the binding free energy for the model antiparallel
P-sheet in w&r --
PIWCIZB
AA
(I) AC-Ala-NHMe(aq) + O(aq) (2) [Ac-Ala-NHMe],(aq) + AC-Ala-NHMe(aq) Binding: Z[Ac-Ala-NHMe](aq) + [Ac-Bla-NHMe],(aq) Binding: 2[formamide](aq) --f [formamide],
--
58.64&@14 6412&0.16t -5.48&02Ej -@34&@28§
Free energy differences are given in kcal/mol. The uncertainties were computed as error propagated standard deviations of 100 averages of 100 datasets. t Peptides were constrained so that T = 1.97 A (the minimum corresponding to the hydrogen-bonded B-sheet on the free energy surface, see below) during the calculation of AA(2). :: AAbind = AA(l)-AA(2). 5 Sneddon et al. (1989).
probability distribution of an observable quantity from what it would be in an unconstrained system. If the distribution is altered significantly, then properties computed from statistical mechanical averages of the quantity, such as free energy differences, should be corrected to account for the constraints. The correct statistical mechanical average of the quantity 0 is (Valleau & Torrie, 1977): (0)
= 10 exP(BU*)>* (exp(BU*)>*
’
(6)
where U* is the constraint potential and the notation (. .)* denotes an average over the constrained ensemble. In our case, the quantity of interest is 0 = exp( -a[ U(r, +dr,) - U(ri)]} (see eqn (1)) a,nd U* is the harmonic potential constraining the backbone dihedral angles (see eqn (5)). We can see in equation (6) that if 0 and exp(fiU*) are statistically independent, then the average of the product in the numerator factors into the product of the averages, and therefore (0) =
