Vol. 53 No. 1/2006, 121–129 Regular paper

Molecular dynamics (MD) is, at present, a unique tool making it possible to study, at the atomic level, conformational transitions in peptides and proteins. Nevertheless, because MD calculations are always based on a more or less approximate physical model, using a set of approximate parameters, their reliability must be tested by comparison with experimental data. Unfortunately, it is very difficult to find a peptide system in which conformational transitions can be studied both experimentally and using MD simulations so that a direct comparison of the results obtained in both ways could be made. Such a system, containing a rigid alpha-helix nucleus stabilized by La(3+) coordination to a 12-residue sequence taken from an EF-hand protein has recently been used to determine experimentally the helix propagation parameters in very short polyalanine segments (Goch et al. (2003) Biochemistry 42: 6840-6847). The same parameters were calculated here for the same peptide system using the peptide growth simulation method with, alternatively, charmm 22 and cedar potential energy functions. The calculated free energies of the helix-coil transition are about two times too large for cedar and even three times too large for charmm 22, as compared with the experimental values. We suggest that these discrepancies have their origin in the incorrect representation of unfolded peptide backbone in solution by the molecular mechanics force fields.

Molecular dynamics (MD) methods are a valuable tool in studies of conformational transitions in peptides and proteins.Only global kinetic and/or thermodynamic parameters of these processes can be determined from experiments, whereas they are precisely described, at the atomic level, by MD simulations.Consequently, the MD methods can give, in principle, answers to all possible questions concerning such processes.Unfortunately, they have one fundamental weakness: energy and forces are calculated on the basis of a model (a so-called molecular mechanics, or MM, model), and while the form and parameters of the model are based on our best understanding of the underlying physical reality, it is always uncertain to what extent the results of the calculations accurately describe the physics of a particular system.To further confuse the issue, several implementations of the MM model have been developed, with equal claim to accuracy.
It is, therefore, essential to calibrate the MM model by establishing the reliability of MD results: only if the calculated global parameters are close enough to the experimental ones we can hope that the theoretical analysis is sound and that the detailed description given by simulations is close to reality.Many examples of calibration studies exist in the literature; however, the reliability of application of the MM model to helix-coil transitions of peptides remains inadequately documented.
An additional requirement of any test of reliability is that simulation results must be precise, in

2006
M. Maciejczyk and others the sense that the effect of statistical noise has been eliminated (by sufficient averaging) from the model values.The results of experimental solution studies are averages over very large numbers of independent molecules, and hence have a negligible statistical error.Simulation studies are typically done with a single solute molecule; this gives a time-averaged behavior, which is then compared with the ensemble average obtained experimentally.The rate of convergence of a time average depends on the correlation over time of the averaged property: a property that fluctuates slowly requires a proportionally longer time for convergence to a given precision than one that fluctuates rapidly (with the RMSD approximately proportional to the square root of the simulation time).Thus, precision requires simulations extending over many times a property's correlation time (Straatsma et al., 1986).
Unfortunately, it is very difficult to find a peptide system in which conformational transitions can be studied both experimentally and using MD simulations, so that a direct comparison of the results obtained in both ways could be made.A major difficulty is the long relaxation time of the helix-coil transition, as is typical for any nucleation-growth process (Schellman, 1955).This time is much shorter for helix formation in more flexible peptides composed of β-amino acids, and this has permitted precise studies of such systems by van Gunsteren and coworkers (Daura et al., 1997;1998;1999).For peptides composed of α-amino acids, calculation results have frequently been compared with experimental data obtained for systems not identical with, but only more or less similar to, those studied by MD methods.The similarity is often so remote that the one conclusion that can be drawn from such comparisons is that the calculation results are not at odds with the experiment.
Because of its relevance for the protein folding process, the helix-coil transition has been studied very extensively.The experimental work concentrated mainly on determination of thermodynamic parameters of the transition in peptides of various lengths and sequences (Doig, 2002).Only peptides with a sufficiently large helix content (in practice >10%) and well soluble in water could be studied.Therefore, they must be long enough -built of at least thirteen, properly chosen amino-acid residues.Such long peptides present a major challenge for MD methods.
The fundamental problem met in MD simulations of helix-coil transitions is adequate probing of the peptide conformational space which is exceedingly large even for medium-sized peptides, built of about ten amino-acid residues.Moreover, the transition of each residue between the helical and coil states is a rare event because of a large en-ergy barrier separating the extended and the helical conformations corresponding, respectively, to the second and the third quadrant of the Ramachandran plot.
Therefore, MD methods have been used mainly to characterize the structural and dynamical properties of the peptide helical conformation and the mechanism of its unfolding (DiCapua et al., 1990;Daggett et al., 1991;Soman et al., 1991;Tirado-Rives & Jorgensen, 1991;Daggett & Levitt, 1992;Lee et al., 2000;Armen et al., 2003).The results of these studies were shown to be compatible with what was known from experimental investigations but quantitative comparisons with any experimental data were not possible.
The helix-coil transitions were simulated in very short peptides built of three (Tobias & Brooks, 1991) or five (Hummer et al., 2000;2001) amino-acid residues.Parallel experimental studies of such peptides are not feasible.Therefore, the calculated kinetic and/or thermodynamic parameters could be only quantitatively compared with information derived from experimental data obtained for much longer peptides.
A first attempt to estimate the free energy of the helix-coil transition in a peptide with a length comparable to that of the shortest peptide systems studied using experimental methods (thirteen amino-acid residues) was made by Daggett and Levitt (1992) from calculated waiting times for the h ® c and c ® h transitions.A similar approach was used recently by Armen et al. (2003).Because in both cases few such transitions were observed during the simulations, the free energy estimates can be treated, at best, as semi-quantitative.
Subsequently, Young and Brooks (1996) were able to calculate the free energy of helix propagation using the umbrella sampling method in a peptide containing as many as fifteen alanine residues.However, the authors met yet another problem: the calculation results strongly depended on the accepted, arbitrary definitions of the helix and coil states.Moreover, the calculated microscopic parameters were, as admitted by the authors, related to, but not identical with, the experimental parameters.
These two fundamental problems met in MD simulations of helix formation, i.e. adequate sampling of the peptide conformational space and the calculation of sufficiently precise thermodynamic parameters of reasonably long oligomers, were apparently solved in a satisfactory way by Wang and coworkers (1995) using the peptide growth simulation (PGS) method and, recently, by García and Sanbonmatsu (2002) using the replica exchange method (Sugita & Okamoto, 1999).Those two studies address the problem of convergence in a completely different manner: the PGS method does not depend on the (slow) equilibration between helix and coil states, while the replica exchange method exploits the much more rapid equilibration that occurs in the high-temperature replicas.Furthermore, the free energy of helix propagation in long polyalanine chains (up to fourteen (Wang et al., 1995) and twenty one (Garcia & Sanbonmatsu, 2002) Ala residues) calculated using these methods should, and was shown to be, close to the analogous parameter determined from experimental studies of polyalanine-based model peptide systems.Nevertheless, the compared theoretical and experimental systems, although close to each other, were still not identical.(These alanine oligomers are poorly soluble in water, and the experimental system is an extrapolation from experiments with molecules containing also amino-acid residues other than alanine.) A peptide system suitable for both experimental and MD studies of the helix-coil transition has been proposed by Siedlecka et al. (1999).The calcium-binding peptide Ac-Loop-NH 2 with the sequence AcDKDGDGYISAAENH 2 strongly coordinates lanthanum ions and as a result assumes a rigid conformation with the ϕ angles of Ala10, Ala11, and Glu12 and the ψ angles of Ser9, Ala10 and Ala11 confined to the third quadrant of the Ramachandran plot.The La 3+ -saturated Loop sequence provides, therefore, a helix nucleus for the polyalanine segments A n attached to its C-terminus in peptides PA n with the sequence AcDKDGDGYISAAE(A) n NH 2 .
The conformation of metal-free peptides Ac-Loop-NH 2 (alternatively denoted as PA 0 ) and PA n , where n = 1 ÷ 4, is close to random coil.Upon La 3+ coordination the A n part of the peptides (reduced to the NH 2 group in the case of PA 0 ) assumes, to a large extent (about 70%), the helical conformation.Recently, the free energies of these partial helixcoil transitions (∆G An ) have been determined from measurements of La 3+ -binding constants to PA n peptides, complemented by NMR measurements of proton exchange rates of their C-terminal NH 2 group (Goch et al., 2003).
In their metal-bound form, PA n peptides are excellent objects for MD studies.They are short, built of at most 16 residues, with twelve of them restricted to a very small fraction of conformational space.Moreover, their conformation is very compact, keeping at the minimum the number of solvent molecules that must be used in MD simulations.
Here we have calculated ∆G An values for a series of peptides PA n using the PGS method (see below) with, alternatively, charmm 22 and cedar potential energy functions (force fields) and compared them with exactly the same thermodynamic parameters determined experimentally.

MATEriALs AND METhoDs
PGs method.PGS (Wang et al., 1995) determine free energies of the transformation of shorter peptides into longer ones, in terms of the Thermodynamic Integration method (Allen & Tildesley, 1987).An atom (group of atoms) is slowly transformed into a different atom (group of atoms) and the free energy of this transformation is calculated according to the equation: where V is a potential energy function depending on the transformation parameter λ in the following way: V(λ = 0) is the potential energy function of the starting component, and V(λ=1) is the potential energy function of the final component.The brackets imply the averaging over a Boltzmann distribution.During the simulation the parameter λ is slowly changed either from 0 to 1 (growth process) or from 1 to 0 (annihilation process) and the integral described by Eqn. 1 is evaluated.For an infinitely long simulation process δG growth = -δG anihilation .In real calculations Eqn. 1 is approximated by: where W represents the work performed on the system in order to force its state to change.
The growth processes analyzed in this work are presented in Fig. 1.The molecule is transformed to one containing an additional Ala residue at the end of the chain, and this is done twice, once with the added residue in an extended conformation, and again with the residue in the α-helical conformation, the free energy difference being accumulated for both processes.The difference between the two free energies is then the free energy difference for the transition of the added residue from extended (i.e., coil) to helical conformation.One sees that no transition between helix and coil state is directly simulated; consequently, the relaxation time of the system is not the considerable relaxation time of the helix-coil equilibrium, but the (apparently much) shorter time associated with rearrangement of the solvent structure around the growing residue.
The starting, artificial peptide P* has the same sequence as peptide PA 0 but with the C-terminal CO-NH 2 group substituted by the neutral carbon atom.The conformation of the metal-binding loop coordinated to an Ln 3+ ion was taken from the PDB (code 1NKF).The lanthanum ion was modeled as a soft sphere with the Lennard-Jones parameters of the Ca 2+ ion and the charge set to +3 elementary charges. (2) (1) A more faithful representation of this system is not possible because the force field parameters for lanthanide ions are not known.For the cedar force field additional restraints on distances between the Ln 3+ ion and the groups coordinating it were applied, preventing their dissociation from the ion.
To keep the dihedral angles of the peptide backbone groups in the desired region of conformational space "flat-bottom" restraints (Hermans et al., 1992) were applied.Such restraints, when placed correctly, effectively divide the probed conformational space and allow one to simulate exclusively one particular conformation, but to not perturb it locally.
Positions of the alanine free energy barriers were located with the potential of mean force calculations of the dihedrals ϕ and ψ of the C-terminal alanine in peptide PA 2 (results not shown).Similar positions of energy barriers have been obtained from free dynamics simulations of a system modeling the alanine dipeptide (Smith, 1999;Hu et al., 2003).Each restraint consists of two ramp potentials placed on the free energy barriers on either side of the free energy minimum.These ramp consisted of a penalty energy changing linearly from 0 to 20 kcal/mol over a 10 o interval.For the helical and extended conformations the backbone dihedrals were restrained to the regions -180 o < φ < 0 o , -140 o < ψ < 20 o and 140 o < φ < 0 o , 20 o < ψ < -140 o of the Ramachandran plot, respectively.simulation procedure.All simulations were carried out with the program Sigma (Mann et al., 2002).Bond lengths were constrained to their equilibrium values with the SHAKE algorithm (Ryckaert et al., 1977).Electrostatic interactions were calculated with a multiple time step method (Tuckerman et al., 1992) with the dual cut-off set to 6 Å and 10 Å. Short range, medium range, and long range (via the particle-mesh Ewald summation method (Darden et al., 1993;Schlick et al., 1999)) interactions were calculated, respectively, at 2, 4, and 12 fs intervals.Simulations were carried out at constant pressure (1 atm) and temperature (300 K) (Berendsen et al., 1984).
Each peptide was placed in the center of a water-filled cubic box and periodic boundary conditions were applied.The box dimensions ranged from 38 up to 44 Å depending on peptide length and conformation (extended or helical) and always exceeded the peptide dimensions by at least 20 Å.
Each calculation of the growth free energy consisted of the following steps: 1) minimization of water energy with fixed solute coordinates; 2) 50 ps equilibration of water with fixed solute conformation; 3) minimization of water energy with fixed solute conformation; 4) minimization of solute energy with fixed solvent conformation; 5) energy minimization of the whole system; 6) 9.6 ps MD simulations at the temperature of 50 K with starting velocities assigned from Boltzmann distribution; 7) five 5.8 ps MD simulations at temperatures increasing from 100 to 300 K in 50 K steps; 8) two peptide growth/ annihilation cycles, each consisting of a) 29 ps equilibration at λ = 0, b) 576 ps growth simulation, c) 29 ps equilibration at λ = 1, d) 576 ps annihilation simulation.
Two force-fields were tested: 1) charmm 22 all-atom force field (MacKerell et al., 1998) distributed with the XPLOR program (Brunger, 1992) combined with the TIP3P (Jorgensen et al., 1983) water model and 2) cedar/gromos all-atom force field distributed in the Sigma package (Ferro et al., 1980;Hermans et al., 1984) combined with the SPC water model.Protein structure files were prepared with the XPLOR program (Brunger, 1992).
Error estimates.The calculated growth free energy, δG for each process is based on the estimates of energy changes, W calculated (with Eqn. 2) for at least two simulations in which a residue was inserted and two simulations in which the same residue was deleted.These values were used to give independent mean and mean-square deviation for the insertion and deletion processes.Because the simulation time is finite, the net energy change for a cyclic process of insertion followed by deletion is not zero, <W i > + <W d > > 0; however, the mean value, W i of the energy for insertion and of -W d (where W d is the mean energy for deletion) can be taken as a

Helix propagation in polyalanine
best estimate of the growth free energy.In addition, an error estimate follows by taking the mean of the mean-square deviations of the insertion and deletion energies (Hermans et al., 1992;Hu et al., 2002).

rEsULTs simulation results
The calculated growth free energies δG i,e , cumulative free energies ΔG i,e , and free energies of single residue transitions from the extended to the helical conformation ΔΔG i,e = ΔG i,e -ΔG i-1,e defined in Fig. 2 are listed in Table 1.Cumulative free energies ΔG i,e are sums of growth free energies δG i,e along the appropriate path (see Fig. 2).

Free energies of helix-coil transitions
The helix-coil transitions observed in experimental studies do not correspond precisely with the transitions between the helical and the extended conformations simulated by the PGS method.Although the coil state is dominated by the extended conformation it includes also some residues with the ϕ and/or ψ angles characteristic for the helical conformation (Wang et al., 1995), if these residues are not a part of a helix nucleus and the amide protons of the peptide groups following their C α carbons do not form the helical (i → i -4 or i → i -3) hydrogen bonds.These states are excluded from the extended conformation by the restraints imposed during the simulations.
The non-helical, coil state of PA 0 peptide (Pc) contains only one, extended conformation (Pe) with the ψ angle of Glu12 kept outside the helical region, because otherwise the terminal NH 2 group would be involved in a helical hydrogen bond.Consequently, the coil state is identical with the extended conformation simulated by the PGS procedure.This is not the case for longer peptides.E.g., the coil state of PA 2 peptide (Pccc), consists not only of the extended conformation Peee but also of three additional conformations: Pehe, Peeh, and Pehh.Consequently, the free energy of this state ΔG 0,c < ΔG 0,e .Similarly, ΔG 1,c < ΔG 1,e because not only the Phee conformation but also the Pheh one is populated in the Phcc state.
Corrections needed to transform ΔG i,e parameters into ΔG i,c ones are small (see Table 1).They can be calculated if populations of all conformations including "h" alanines in the respective coil states are known.We estimated them from the data obtained by Hu et al. (2003) from MD simulations of the alanine dipeptide.For both the charmm 22 and cedar force fields they calculated the equilibrium constants K between the alanine residue conformations with helical and non-helical ϕ and ψ dihedrals: K = 1 and 0.29 for charmm 22 and cedar, respectively.
Free energy levels corresponding to additional h conformations are shifted down by the value equal to B ln l k T K below the extended conformation.l is the number of 'e' alanines replaced by 'h' ones and k B is Boltzmann's constant.From statistical thermodynamics (Hill, 1986) it then follows that Free energies of single residue helix-coil transitions ΔΔG i,c = ΔG i,c-ΔG i-1,c , analogous to the ΔΔG i,e parameters, were calculated from the ΔG i,e values using Eqn. 4 and are listed in Table 1.Extended, helical, and partly helical conformations of C-terminal segments of PA n peptides are described by chains of e and/or h subscripts.E.g., the symbol P hhee denotes peptide PA 3 with its loop conformation restricted by Ln 3+ ion binding, Ala13 and Ala14 in the helical conformation, and Ala15 and the terminal NH 2 group in the extended conformation.Growth simulations are indicated by thick arrows and growth free energies corresponding to them are denoted by δG i,e where n is the number of alanine residues created by the PGS procedure, i denotes the total number of peptide C-terminal groups, including NH 2 , in the helical conformation, and e indicates that the remaining groups are in the extended conformation, even if their number is zero, as in the case of δG n+1,e symbols.
Similarly, the cumulative free energy of a series of growth steps transforming P* into peptide PA n with i terminal groups in the helical conformation is denoted by ΔG i,e symbol.Thin vertical arrows indicate the conformational transitions of a single alanine residue or NH 2 group from the extended to the helical conformation.The free energies of these transitions ΔΔG i,e = ΔG i,e -ΔG i,e .Since the free energies of the coil conformation of A n segments of apo-PA n peptides equal ΔG 0,c the free energy changes ΔG An related with the conformational transitions induced in them by Ln 3+ coordination are described by the equation: The ΔG An values calculated from Eqn. 6 and the ΔΔG i,c parameters listed in Table 1 are presented in Table 2.

helix populations in PA n peptides
The population of helical segments built of j residues is given by the equation: Table 1.Results of PGS obtained with charmm 22 and cedar force fields.δG i,e , ΔG i,e , and ΔΔG i,e and parameters are defined in Fig. 2. ΔΔG i,e parameters have been transformed into the free energies of the respective helix-coil transitions (ΔΔG i,c ) as described in the text.The conformation of C-terminal segments of apo-PA n peptides is close to random coil (Goch et al., 2003).Upon La 3+ binding a large fraction of their residues assume the helical conformation.The helix builds up from its nucleus formed in the Loop sequence by coordination of the metal ion.The conformational state of A n segments can be described as a dynamical equilibrium between the coil conformation and an assembly of helical segments containing from 1 up to n + 1 backbone groups (n alanine residues and the C-terminal NH 2 group).The free energy of such a multiple-equilibrium system ΔG met (see Fig. 3) is given by the equation (Poland & Scheraga, 1970):   (2003).(b) Obtained from ΔΔG tn values determined from La 3+ -binding constants to peptides PA n given in the last column of Table 1 in Goch et al. (2003).ΔG An = ΔG A0 + ΔΔG tn .The error limits have been estimated from the errors of binding constant determinations listed in the first column of Table 1 in Goch et al. (2003).Using Eqn. 6 the molar fractions of full helices including the C-terminal NH 2 groups (x nNh = p j+1 ) and the average populations of amino-acid residues in the helical conformation were calculated and compared with the experimental values in Table 2.

DisCUssioN
The importance of a "direct" comparison of MD results with experimental data was stressed in the introduction.Let us specify what it actually means.Raw data of any experimental or theoretical study must be analyzed and interpreted to give physically meaningful results.Such a procedure is based on accepted theoretical models and on more or less solid assumptions.Although indispensable, such an analysis in some way always reduces and distorts the information contained in the data.Therefore, the data obtained using different methods should be compared at the lowest possible level of interpretation.
The crucial merit of our peptide system is that after the very first steps of data analysis the experimental and theoretical results can be compared.Only one assumption, i.e., that the interaction of the An and Loop segments with each other in metalfree PA n peptides is negligible, was necessary to determine the ΔG exp values from the experimental data (Goch et al., 2003) without using any explicit helix-coil transition theory.On the other hand, the model's ΔG An parameters (see Table 1) have been calculated using Eqn.6, based only on fundamental An An principles of statistical thermodynamics, from δG i,e the values obtained from PGS, with only a slight correction being needed to transform these into the δG i,c parameters.
The statistical errors of δG i,c (see Table 1) lead to errors in the determined ΔG An values that can be evaluated at about 0.3 kcal/mol (PA 0 ) to 1 kcal/mol (PA 4 ), which is of the same order of magnitude as the experimental values.Clearly, the peptide systems studied by us, although very small, need much longer simulation times, which, at present, are unacceptable because of the limits on available computing power.However, the systematic errors due to the use of imperfect force field parameters are much larger: the ΔG An values calculated using cedar and charmm 22 exceed the experimental ones by as much as 3 kcal/mol.
The free energy is, of course, the fundamental parameter describing the helix-coil transition.Nevertheless, the precision of its determination does not correlate linearly with the precision of the prediction of the helix population which is often used, in practice, to estimate the value of theoretical calculations.Whereas the helix populations (x nNh and x nh ) determined using the cedar force field are very close to the experimental data, much larger deviations are observed when charmm 22 is used (see Table 2).
We suggest that the serious discrepancies between the results of MD simulations and experimental data have their origin in the incorrect representation of unfolded peptide backbone in solution (Hu et al., 2003) by the MM force fields.Hu and coworkers presented results of MD simulations of Ace-Ala-Nme and Ace-Gly-Nme "dipeptides" in solution, using several force fields including charmm22 and cedar, and also with a mixed QM/MM force field, which they performed with the SCC-DFTB method (Elstner et al., 1998).One result was the great variation of the conformational distributions obtained with the dif- ΔG i,c parameters were calculated from the corresponding ΔG i,e values (see Fig. 2 and Table 1) using Eqn. 3. ΔG met describe free energies of multiple equilibrium systems composed of n + 2 conformational states.

Figure 1 .
Figure 1.illustration of PGs.a) In P* → PA 0 transformation the starting, artificial peptide is terminated with a neutral carbon which is slowly replaced by CONH 2 group; b) The C-terminal amide NH 2 group of peptide PA 0 is replaced by terminally amidated alanine in the PA 0 → PA 1 transformation.The segment built of amino acids 1-11 is denoted as aa 1-11 .The transformed atoms are shown in bold characters.Flat-bottom restraints placed on the indicated ϕ and ψ dihedrals keep the growing groups in the desired (helical or extended) conformation.
Figure 2. schematic presentation of PGs performed in this work.Extended, helical, and partly helical conformations of C-terminal segments of PA n peptides are described by chains of e and/or h subscripts.E.g., the symbol P hhee denotes peptide PA 3 with its loop conformation restricted by Ln 3+ ion binding, Ala13 and Ala14 in the helical conformation, and Ala15 and the terminal NH 2 group in the extended conformation.Growth simulations are indicated by thick arrows and growth free energies corresponding of simulated and experimental results: free energies of the helix-coil transitions induced in C-terminal segments of peptides PA n by La 3+ binding calculated using charmm 22 and cedar force fields (ΔG An ) and determined from experiments (ΔG An ), full-helix populations in An segments of La 3+ -saturated peptides PA n calculated (x nNh ) and determined experimentally (x nNh ), fractional population of helical hydrogen bonds calculated from the MD simulations (x nh ) and the experimental data (x nh

Figure 3 .
Figure 3. schematic presentation of free energy changes (ΔG An ) related to conformational transitions induced in the C-terminal segments of peptides PA n by La 3+ ion coordination.ΔG i,c parameters were calculated from the corresponding ΔG i,e values (see Fig.2and Table1) using Eqn. 3. ΔG met describe free energies of multiple equilibrium systems composed of n + 2 conformational states.

n Free energies of conformational transitions induced in C-terminal segments of PA n peptides by La 3+ binding
).
(a)Determined from NMR measurements of the exchange rates of C-terminal NH 2 protons in peptides PA n .See Table2inGoch et al.