Temperature and Pressure Dependence of Alanine Dipeptide Studied

12038

J. Phys. Chem. B 2008, 112, 12038–12049

Temperature and Pressure Dependence of Alanine Dipeptide Studied by Multibaric-Multithermal Molecular Dynamics Simulations Hisashi Okumura* and Yuko Okamoto† Department of Physics, School of Science, Nagoya UniVersity, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan ReceiVed: December 27, 2007; ReVised Manuscript ReceiVed: May 6, 2008

We applied the multibaric-multithermal (MUBATH) molecular dynamics (MD) algorithm to an alanine dipeptide in explicit water. The MUBATH MD simulation covered a wide range of conformational space and sampled the states of PII, C5, RR, RP, RL, and Cax 7 . On the other hand, the conventional isobaric-isothermal simulation was trapped in local-minimum free-energy states and sampled only a few of them. We calculated the partial molar enthalpy difference ∆H and partial molar volume difference ∆V among these states by the MUBATH simulation using the AMBER parm99 and AMBER parm96 force fields and two sets of initial conditions. We compared these results with those from Raman spectroscopy experiments. The Raman spectroscopy data of ∆H for the C5 state against the PII state agreed with both MUBATH data with the AMBER parm96 and parm99 force fields. The partial molar enthalpy difference ∆H for the RR state and the partial molar volume difference ∆V for the C5 state by the Raman spectroscopy agreed with those for the AMBER parm96 force field. On the other hand, ∆V for the RR state by the Raman spectroscopy was consistent with our AMBER-parm99 force-field result. All the experimental results fall between those of simulations using AMBER parm96 and parm99 force fields, suggesting that the ideal force-field parameters lie between those of AMBER parm96 and parm99. 1. Introduction There are difficulties in discussing theoretically the temperature and/or pressure dependence of biomolecules. The first problem is that they have complicated free energy surfaces with many local minima. Conventional molecular dynamics (MD) and Monte Carlo (MC) simulations in physical ensembles, such as the canonical1–4 and isobaric-isothermal5,6 ensemble, tend to get trapped in these local-minimum states. It causes the simulations to sample only a limited range of conformational space and give inaccurate results. One of the powerful techniques to avoid this difficulty is generalized-ensemble algorithms.7–9 One of the most well-known generalized-ensemble algorithms is the multicanonical algorithm.10–13 In the multicanonical ensemble, a free one-dimensional random walk is realized in the potential-energy space. Thus, a simulation with this algorithm does not get trapped in free-energy-minimum states and is able to sample a wide range of the conformational space. It is also possible to obtain various canonical-ensemble averages in a wide range of temperature from a single simulation run by the reweighting techniques.14 However, because the multicanonical simulation is performed in a fixed volume, the pressure cannot be specified as in experimental environment. To investigate the temperature and/ or pressure dependence, not only temperature but also the pressure must be specified. We have recently proposed multibaric-multithermal (MUBATH) MC15–17 and MD18–20 algorithms to overcome this difficulty. This is an extension of the multicanonical simulation, which can perform twodimensional random walks not only in the potential-energy space but also in the volume space. The MUBATH algorithm has the following advantages: (1) It allows the simulation to escape from * Corresponding author. E-mail: [email protected]. † E-mail: [email protected].

local-minimum-energy states and to sample the conformational space more widely than the conventional isobaric-isothermal method. (2) One can obtain various isobaric-isothermal ensembles not only at any temperature as in the multicanonical algorithm but also at any pressure from only one simulation run. (3) One can control temperature and pressure so that we can compare the simulation results with those by experiments under the same conditions. The replica-exchange method21,22 is another powerful generalized-ensemble algorithm which is used as an alternative to the multicanonical algorithm. Similarly, the replica-exchange method in the isobaric-isothermal ensemble23–27 can be used as an alternative to the multibaric-multithermal algorithm. As far as the weight factors can be obtained, the multicanonical version is more efficient in sampling than the replica-exchange version,28 and this is why we used MUBATH in the present work. In a preceding article,20 we have applied the MUBATH MD algorithm to an alanine dipeptide in explicit water. It was the first work to apply the MUBATH MD algorithm to a biomolecular system. The AMBER parm96 force field29 was used for the peptide. However, discussion on the force-field dependence is necessary to investigate the temperature and pressure dependence of the peptide in more detail. In this article, both AMBER parm9930 and AMBER parm9629 are employed. The AMBER parm99 force field is known to tend to form R-helices, whereas the AMBER parm96 tends to make β-sheets in general.31 The alanine dipeptide has been studied theoretically extensively as a model biomolecule because of its simplicity (see, for example, refs 32–37). Experimentally, Takekiyo et al. have recently performed Raman spectroscopy experiments under several temperature and pressure conditions.38 They discussed partial molar enthalpy difference and partial molar volume difference among some states of alanine dipeptide. Estimation

10.1021/jp712109q CCC: $40.75  2008 American Chemical Society Published on Web 08/30/2008

Temperature/Pressure Dependence of Alanine Dipeptide

J. Phys. Chem. B, Vol. 112, No. 38, 2008 12039

of the partial molar enthalpy and partial molar volume is important in solution chemistry because these quantities control the population of each state when temperature and pressure are changed. There is no simulational study, except for our preceding article,20 that is related to these experiments. In this article, we discuss the temperature and pressure dependence of the peptide using the two force fields and compare the physical properties of the alanine dipeptide obtained by the Raman spectroscopy. To check the initial condition dependence, we used two sets of initial conditions in the present simulations, whereas only one set of initial conditions was employed in the preceding article.20 We show that the conventional isobaric-isothermal MD simulation depends on the initial conditions. On the other hand, the MUBATH MD simulation gives the results which are independent of initial conditions. In section 2 the algorithm of the MUBATH MD method is explained. Computations of the alanine dipeptide simulations in explicit water are described in section 3. Results are discussed in section 4. Section 5 is devoted to conclusions. 2. Methods In the isobaric-isothermal ensemble,5,6 the distribution PNPT(E,V) of potential energy E and volume V is given by

PNPT(E,V) ) n(E,V)e-β0H

(1)

where n(E,V) is the density of states as a function of E and V, β0 is the inverse of the product of the Boltzmann constant kB and absolute temperature T0, at which simulations are performed, and H is the “enthalpy” (without the kinetic energy contributions): H ) E + P0V. Here, P0 is the pressure at which simulations are performed. This ensemble has bellshaped distributions both in the potential-energy space and in the volume space if the density of states is well behaved or an increasing function of E and V and it is so in most of the biomolecular systems. To obtain the isobaric-isothermal ensemble, the combination of the Nose´ thermostat1,2 and the Andersen barostat5 is frequently employed. In the MUBATH ensemble,15–19 every state is sampled with a weight factor Wmbt(E,V) ≡ exp{-β0Hmbt(E,V)} so that a uniform distribution of both E and V may be obtained:

Pmbt(E,V) ) n(E,V) Wmbt(E,V) ) constant

(2)

Here, Wmbt(E,V) and Hmbt are referred to as the MUBATH weight factor and the MUBATH enthalpy, respectively. The difference between Hmbt and H is written as δH(E,V): Hmbt(E,V) ) H + δH(E,V). The difference δH(E,V) is therefore characteristic of the MUBATH simulation. The case of δH(E,V) ) 0 gives the regular isobaric-isothermal ensemble. Let us consider a system consisting of both spherical atoms and rigid-body molecules. The number of atoms in the flexible peptide molecule is N and that of the rigid-body water molecules is M. Atoms in the flexible model molecule are treated in the same way as monatomic molecules in the MD simulations. In the preceding article, the MUBATH algorithm was proposed on the basis of the Nose´ thermostat,1,2 the Andersen barostat,5 and the rigid-body-molecule Hamiltonian. In the present article, we describe a formalism based on the Nose´-Poincare´ thermostat39,40 instead of the Nose´ thermostat. The combination41 of the Nose´-Poincare´ thermostat,39,40 the Andersen barostat,5 and the symplectic quaternion scheme42 provides a symplectic MD integrator. The symplectic integrator has the advantage that the secular deviation of the value of the Hamiltonian is suppressed.43

The isobaric-isothermal Hamiltonian Hmbt for the present system is given by

Hmbt ) s[HNA(r˜ ,p ˜,q,π˜ ,V,PV,s,Ps) - H0]

(3)

where N+M

HNA )

∑ i)1

p˜ i2

M

∑

1 T + π˜ j S(qj)DjST(qj)π˜ j + 2miV2/3s2 j)1 8s2 PV2 Ps2 + P0V + + gkBT0 log s + 2W 2Q

E(V1/3r˜ ,q) +

δH(E(V1/3r˜,q),V) (4) Here, the superscripts T in π˜ jT and ST stand for transpose, the summation with respect to i is taken over both N atoms and M rigid-body molecules and that with respect to j over only M rigid-body molecules. The variable r˜i is the coordinate scaled by the length of the simulation box: r˜i ) V-1/3ri (ri is the real coordinate). We have introduced a simplified notation for the set of the scaled coordinates: r˜ ≡ {r˜1, r˜2,..., r˜N+M}. When i indicates a rigid-body molecule, ri stands for the coordinate of its center of mass. The variable p˜i is the conjugate momentum for r˜i. The real momentum pi is related to the virtual momentum p˜i by pi ) p˜i/V1/3s, where s is the additional degree of freedom for the Nose´ thermostat. The variables PV and Ps are the conjugate momenta for V and s, respectively. The constant mi is the mass of atom i or rigidbody molecule i. The constants W and Q are the artificial “mass” related to V and s, respectively. In the case of the system consisting of N atoms and M rigid-body molecules, g equals 3N + 6M. The variable q is a quaternion which specifies the orientation of the rigid-body molecule. Here, the quaternion q ) (q0, q1, q2, q3)T is related to the Euler angle (φ, θ, ψ) as follows:

( θ2 ) cos( φ +2 ψ ) θ φ-ψ q ) sin( ) cos( 2 2 ) θ φ-ψ q ) sin( ) sin( 2 2 ) θ φ+ψ q ) cos( ) sin( 2 2 )

q0 ) cos

(5)

1

(6)

2

(7)

3

(8)

(

)

The elements of the matrix S(q) are given by

q0 -q1 -q2 -q3 q1 q0 -q3 q2 S(q) ) q2 q3 q0 -q1 q3 -q2 q1 q0

(9)

The variable π˜ j is the conjugate momentum for qj. The real momentum of the quaternion πj is related to π˜ j by πj ) π˜ j/s. The matrix D is a 4 × 4 matrix consisting of the inverse of the principal moments of inertia I1, I2, and I3:

D)

(

I0-1

0

0

I1-1

0

0

0

0

I2-1

0

0

0

0

I3-1

0

0

)

(10)

where I0 is an artificial constant that gives the correct equations of motion in the limit of I0 f ∞.

12040 J. Phys. Chem. B, Vol. 112, No. 38, 2008

Okumura and Okamoto

The canonical equations of motion are obtained in virtual time from the Hamiltonian based on the Nose´ thermostat.20 On the other hand, the canonical equations of motion are obtained in real time from the Hamiltonian based on the Nose´-Poincare´ thermostat39,40 in eq 3. To write the equations of motion more elegantly, we may introduce the angular velocity

ω ) (ω1, ω2, ω3)T

(11)

and the four-dimensional angular velocity

ω ) (0, ω1, ω2, ω3) (4)

T

(12)

where ω1, ω2, and ω3 are the angular velocity along each of the corresponding principal axes. Here, ω and ω(4) are treated as a 3 × 1 matrix and a 4 × 1 matrix, respectively. The angular velocity ω(4) is related to π by

1 ωj(4) ) DjST(qj)πj 2

(13)

Transforming r˜ to r, p˜ to p, and π˜ to ω, the equations of motion are given by

˙ri )

˙ pi V + ri mi 3V

(14)

(

)

˙ V ∂δH s˙ Fi + p ∂E 3V s i 1 ˙j ) S(qj)ωj(4) q 2 ∂δH s˙ ˙j ) 1 + N - ωj × (Ijωj) - Ijωj Ijω ∂E j s PV ˙ V)s W

(

)

p ˙i ) 1 +

(

[ {∑

1 ˙ PV ) s 3V

N+M i)1

)

pi2 ∂δH + 1+ mi ∂E

(

(15) (16) (17)

N+M

˙ Ps )

∑ i)1

Pmbt(E,V) Wmbt-1(E,V)e-β(E+PV)

∫ dV ∫ dE Pmbt(E,V) Wmbt-1(E,V)e-β(E+PV)

(22)

The expectation value of a physical quantity A at T and P is then obtained from

〈A〉NPT ) )

∫ dV ∫ dE A(E,V) PNPT(E,V;T,P) 〈A(E,V) Wmbt-1(E,V)e-β(E+PV) 〉mbt 〈Wmbt-1(E,V)e-β(E+PV) 〉mbt

(23)

where 〈· · ·〉mbt is the MUBATH ensemble average. Because of the random walks both in the potential-energy space and in the volume space, we can calculate physical quantities in wide ranges of T and P. In the case of a large system, the isobaric-isothermal distribution at some temperature T and pressure P is narrow. Thus, the ratio of samples which are utilized to calculate the isobaric-isothermal distribution in eqs 22 and 23 to total samples decreases. We have to perform a long simulation as the system size becomes large or the uncertainties in ensemble averages calculated in eq 23 become large. This problem occurs not only in the multibaric-multithermal algorithm but also in other generalized-ensemble algorithms such as the multicanonical algorithm and the replica-exchange method.

) ∑ F · r } - (P + ∂δH ∂V )]

3. Computational Details

N+M

i

i

0

i)1

Ps Q

M pi2 + ωTI ω - gkBT0 mi j)1 j j j

∑

PNPT(E,V;T,P) )

(18)

(19) s˙ ) s

in the isobaric-isothermal ensemble at the desired temperature T and pressure P is given by

(20) (21)

where I is the 3 × 3 diagonal matrix. Its diagonal elements are the principal moments of inertia I1, I2, and I3. The vector Fi stands for the force acting on atom i if i indicates an atom in the flexible molecule or the total force on all the atoms in molecule i if i indicates a rigid-body molecule. The vector Nj is the torque on rigid molecule j. Performing the MD simulation by the equations of motion in eqs 14–21, the MUBATH distribution Pmbt(E,V) in eq 2 is obtained. In the case of δH(E,V) ) ∂δH(E,V)/∂E ) ∂δH(E,V)/∂V ) 0, the Hamiltonian in eq 3 and the equations of motion in eqs 14–21 become those for the regular isobaric-isothermal MD simulation of the Nose´-Poincare´-Andersen formulation. The MUBATH weight factor is, however, not a priori known and has to be determined, for example, by iterations of short simulations.44–46 After an optimal weight factor Wmbt(E,V) is determined, a long production run is performed for data collection. The reweighting techniques14 are used for the results of the production run to calculate the isobaric-isothermalensemble averages. The probability distribution PNPT(E,V; T, P)

We applied the MUBATH MD algorithm to a system consisting of one alanine dipeptide molecule ((S)-2-(acetylamino)-N-methylpropanamide) and 63 water molecules. We used enough water molecules so that the alanine dipeptide molecule was always held perfectly within the simulation box. We used both AMBER parm9930 and AMBER parm9629 force fields for the alanine dipeptide molecule and the TIP3P47 rigidbody model for the water molecules. We employed a cubic unit cell with periodic boundary conditions. The electrostatic potential was calculated by the Ewald method. We calculated the van der Waals interaction, which is given by the Lennard-Jones 12-6 term, for all pairs of the atoms within the minimum image convention instead of introducing the spherical potential cutoff. The time step was taken to be ∆t ) 0.5 fs. To obtain an optimal MUBATH weight factor Wmbt(E,V), two-dimensional versions of both the iterative method44–46 and the energy landscape paving method,48 which is a special case of the Wang-Landau techniques,49 for multicanonical weightfactor determination were employed. This is because this combination is more efficient to determine Wmbt(E,V) than employing only the iterative method. This technique updates Wmbt(E,V) at every MD step n as Wmbt(E,V, n) ≡ exp[-{E + P0V + δH(E,V,n)}/kBT0] to lower energy barriers. First, we updated δH(E,V,n) using only a histogram of E as in the isobaric-multithermal simulation.16 Thus, δH(E,V,n) is a function of only E and n: δH(E,n). The difference δH(i)(E,n) between the MUBATH enthalpy and the normal enthalpy at time step n in the ith iteration is given by48,49



δH(i)(E, n) )

{

δH(i-1)(E,V) )

δH(i-1)(E) + kBT0 log P(i)(E, n) for P(i)(E, n) g 1 for P(i)(E, n) ) 0 (24) δH(i-1)(E)

where P(i)(E,n) is the histogram of E accumulated until one time step before n, that is, until time step (n - 1). Here, the time step in the ith iteration takes the values n ) 1, 2,..., N(i) t . For the first iteration (i ) 1), the MUBATH enthalpy starts from the normal enthalpy with

δH(0)(E) ) 0

(25)

In the case of i g 2, δH(i-1)(E) is given by the histogram of potential energy P(i-1)(E) accumulated until the last time step in the (i-1)th iteration as follows:

{

δH(i-1)(E) ) δH(i-2)(E) + kBT0 ×

[logP(i-1)(E) - log P(i-1)(E0)]

for E e E0 and P(i-1)(E) g 1

δH(i-2)(E) - kBT0 log P(i-1)(E0) for E e E0 and P(i-1)(E) ) 0 for E > E0 δH(i-2)(E) ) 0 (26)

where log P(i-1)(E) was shifted by log P(i-1)(E0) for potential energy less than a threshold E0 (E e E0) to make the potentialenergy distribution wider toward the lower value of E. This is because it is easier to widen the potential-energy distribution at high temperature toward the lower values of E than to widen the distribution at low temperature toward the higher values of E. We set the threshold as E0 ) -400 kcal/mol, the temperature as T0 ) 600 K and the pressure as P0 ) 300 MPa. To prevent vaporization of the water, the pressure was set at such a high value. We iterated this procedure 18 times (i ) 1-18). In the first iteration the MD simulation was performed for 10 ps (N(1) t ) 20 000). At later iterations, we carried out the MD simulations for longer time, e.g., for 180 ps (N(18) ) 360 000) in the 18th t iteration. These numbers of time steps are the same for both AMBER parm99 simulation and AMBER parm96 simulation.

{

δH(i-2)(E,V) + kBT0 log P(i-1)(E,V) for P(i-1)(E,V) g 1 for P(i-1)(E,V) ) 0 δH(i-2)(E, V) (28)

The time steps for 400 ps in the 19th iteration and 500 ps in the 20th iteration were necessary in the AMBER parm96 simulation. On the other hand, 200 ps in the 19th iteration and 300 ps in the 20th iteration were enough in the AMBER parm99 simulation. An optimal weight factor, Wmbt(E,V), was obtained after these preliminary simulations in total for 2.21 ns in the AMBER parm99 simulation and for 2.61 ns in the AMBER parm96 simulation. Although we used an initial condition of φ ) ψ ) 180° for the weight-factor determination, two sets of the initial values of the alanine-dipeptide backbone dihedral angles were prepared for the data-production runs. Initial condition 1 was φ ) ψ ) 180° and initial condition 2 was φ ) 180° and ψ ) 0°. These initial conformations are shown in Figure 1. In a conformation of φ ) 0° and ψ ) 0°, the oxygen atom of the acetyl group is too close to the hydrogen which bonds to the nitrogen in the N-methyl group. In a conformation of φ ) 0° and ψ ) 180°, the oxygen atom of the acetyl group is again too close to the oxygen in the alanine. Due to these close distances, the potential energy was too high to be used as initial conditions. We performed a MUBATH MD simulation for an equilibration process for 10 ps before sampling from each initial condition with each force field. We then performed a long MUBATH MD simulation for data collection. The simulation time was 20 ns with each force field for each initial condition. Summing up these time steps for the two initial conditions, the total time steps were 40 ns with each force field. For convenience, we discretized E and V into bins and used histograms to calculate δH(E,V). We chose bin sizes of ∆E ) 20 kcal/mol and ∆V ) 0.1 nm3. To calculate the derivative of δH(E,V), such as ∂δH(E,V)/∂E and ∂δH(E,V)/ ∂V, it is necessary to interpolate the histogram. We interpolated the histogram by using the third-order polynomial in both E and V to make these derivatives continuous at the boundaries between two adjacent bins as follows:19

After we elongated δH(E) using only the potential energy distribution P(E), we now widen δH(E,V) by the histogram P(E,V) of both E and V as follows:

δH(i)(E,V,n) )

{

δH(i-1)(E,V) + kBT0 log P(i)(E,V,n) for P(i)(E,V,n) g 1 for P(i)(E,V,n) ) 0 δH(i-1)(E, V) (27)

We iterated the procedure in eq 27 twice, which corresponds to the 19th and 20th iterations. In the 19th iteration, we used δH(i-1)(E,V) calculated from log P(i-1)(E) in the 18th iteration shifted by log P(i-1)(E0) in eq 26. On the other hand, we did not shift log P(i-1)(E) in the 20th iteration and the subsequent MUBATH MD simulation for data collection because δH(i-1)(E,V) has been widened in both the E space and the V space as follows:

where f(x,y) is the value of δH(E,V):


f(x,y) ) δH(E,V)

Okumura and Okamoto

(30)

Here, x and y are defined by

x ) (E - nx∆E)/∆E

(31)

y ) (V - ny∆V)/∆V

(32)

where nx and ny are the largest integer less than E/∆E and that less than V/∆V, respectively. The ranges of x and y are 0 < x e 1 and 0 < y e 1, respectively. The subscripts x and y in fx and fy mean x-derivative of f and y-derivative of f, respectively. The values of f, fx, and fy at (x, y) ) (0,0), (1,0), (0,1), and (1,1) are given by the histogram for the MUBATH MD simulation. Deviation of the value of the Hamiltonian from its initial value is suppressed by this polynomial. For details, see ref. 19. To compare the MUBATH simulation with the conventional isobaric-isothermalone,wealsoperformedtheisobaric-isothermal MD simulations of the same system, which consists of an alanine dipeptide molecule in 63 water molecules. The initial conditions of these simulations were the same as those for the MUBATH simulation. The simulations were carried out for 2 ns each at (T0, P0) ) (240 K, 0.1 MPa), (298 K, 0.1 MPa), and (298 K, 300 MPa). 4. Results and Discussion We present time series of potential energy E, volume V, and backbone dihedral angles φ and ψ in Figures 2–6. These figures are obtained from the AMBER parm99 simulations. Those from the AMBER parm96 simulations were discussed in the preceding paper.20 The time series of E from initial condition 1 are shown in Figure 2. The potential energy fluctuated in a narrow range in the conventional isobaric-isothermal MD simulation. On the other hand, Figure 2d shows that the MUBATH MD simulation covered a wide energy range. The range was 4-5 times wider than that by the isobaric-isothermal MD simulations. Figure 3 shows the time series of V. The MUBATH MD simulation performed a random walk in a wide range of V, as shown in Figure 3d. It was 3-5 times wider than that by the isobaric-isothermal MD simulations. Such a large fluctuation of V can be realized by the MUBATH algorithm. It is not possible by the multicanonical algorithm because the volume is fixed. The isobaric-isothermal MD simulations covered only half a turn (a 180° range) of φ and did not make a complete rotation as shown in Figure 4a-c. This fact means that the simulations were trapped in local-minimum free-energy states and did not rotate beyond the energy barriers. In the MUBATH MD simulation, on the other hand, φ underwent a complete rotation (covering 360 °) four times (four turns) in the time range of 2 ns as shown in Figure 4d. It made 26 turns throughout the simulation of 20 ns. Here, we define one turn of rotation when φ undergoes all values between -180° and 180° and count another turn when φ undergoes all the values after that. Similarly to the simulations from initial condition 1, every isobaricisothermal MD simulation that started from initial condition 2 did not undergo all values of φ. On the other hand, the MUBATH MD simulation made a complete rotation 32 times during 20 ns. The numbers of rotations by the MUBATH MD simulations are listed in Table 1. The time series of E, V, and φ from initial condition 2 were essentially the same as those from initial condition 1. The difference between the two initial conditions lies in ψ. The time series of ψ that started from initial condition 1 and that from initial condition 2 are shown in Figure 5 and in Figure 6,

Figure 1. Initial conformations of alanine dipeptide for the MD simulations. The dihedral angles are (a) φ ) ψ ) 180° and (b) φ ) 180° and ψ ) 0°. The figures were created with RasMol.50

respectively. The isobaric-isothermal MD simulation at T0 ) 298 K and P0 ) 0.1 MPa showed relatively better sampling among the isobaric-isothermal simulations that we performed. It sampled states both around ψ ) 180° and around ψ ) 0°, as shown in Figures 5a and 6a. However, the dihedral angle ψ does not rotate completely at any temperature and pressure in the isobaric-isothermal MD simulation as shown in Figures 5a-c and 6a-c. This feature is similar to the results of φ. As the pressure increases or the temperature decreases, sampling efficiency gets worse. In the case of high pressure, ψ fluctuated around 0° soon after the simulation started (20 ps) from ψ ) 180° and never came back to a state around 180°, as shown in Figure 5b. When the simulation was started from initial condition 2 (Figure 6b), it sampled a state around 180° only once. In the case of low temperature, ψ takes values around 180° for 730 ps and then move to a state around 0°, as shown in Figure 5c. It also never came back to a state around 180°. When the simulation was started from 0°, it did not sample a state around 180° at all, as shown in Figure 6c. On the other hand, the MUBATH simulation showed high sampling efficiency for ψ. As shown in Figures 5d and 6d, there seems to be a potential energy barrier around ψ ) -110°, which cannot be overcome by the isobaric-isothermal simulation (this energy barrier corresponds to the boundaries between the localminimum states defined in Table 3 below). The MUBATH simulation overcame this potential energy barrier and realized complete rotations. The numbers of the rotations of ψ are also listed in Table 1. Such wide rotation ranges of φ and ψ show that the simulation in the MUBATH ensemble can escape from local-minimum free-energy states and samples the conformational space more widely than the conventional isobaric-isothermal ensemble. There are two reasons for the high sampling efficiency in the MUBATH simulation. One reason is that the sampling in a wide range of potential energy makes the simulation overcome the potential energy barrier between the states in the dihedral angle space. The other is that when V is large, the number density of water molecules is less than that when V is small. This means a decrease in the obstacles (water molecules) around the peptide and makes the rotation of φ and ψ easier. Although the multicanonical algorithm has high sampling efficiency due to the former effect, the latter arises from the MUBATH algorithm. We can see high sampling efficiency of the MUBATH algorithm also in the probability distributions P(φ,ψ) of φ and ψ. We discretized φ and ψ using the bin sizes of ∆φ ) ∆ψ ) 10° to calculate P(φ,ψ) as a histogram. Figure 7 shows the distributions P(φ,ψ) and the contour maps of log P(φ,ψ), for example, at T ) 240 K and P ) 0.1 MPa. Although the isobaric-isothermal distribution P(φ,ψ) from initial condition 1 (Figure 7a) has four peaks of the RR, RP, PII, and C5 states, the peaks of the PII and C5 states are too large when compared to the MUBATH distribution P(φ,ψ) as shown in Figure 7b. The distribution P(φ,ψ) from initial condition 2 has only two peaks of the RR and RP states, as shown in Figure 7c. On the other hand, the MUBATH MD simulation



Figure 2. Time series of potential energy E with the AMBER parm99 force field and initial condition 1 from (a) the isobaric-isothermal MD simulation at T0 ) 298 K and P0 ) 0.1 MPa, (b) the isobaric-isothermal MD simulation at T0 ) 298 K and P0 ) 300 MPa, (c) the isobaric-isothermal MD simulation at T0 ) 240 K and P0 ) 0.1 MPa, and (d) the multibaric-multithermal MD simulation.

Figure 3. Time series of volume V with the AMBER parm99 force field and initial condition 1. See the caption of Figure 2 for further details.

from each initial condition provides correct distributions P(φ,ψ) with all the five peaks of the RR, RP, PII, C5, and RL states, as shown in Figure 7b,d. These results indicate again that the MUBATH simulation can get over the free energy barriers and has a high sampling efficiency. Although we have shown that φ and ψ rotated better in the MUBATH simulation than in the conventional isobaricisothermal simulation, the other backbone dihedral angles ω around the amide bond between the acetyl group and the alanine and that between the alanine and the N-methyl group did not rotate freely even in the MUBATH simulation. They fluctuated

TABLE 1: Numbers of the Complete Rotations of the Backbone Dihedral Angles O and ψ by the MUBATH MD Simulations for 20 ns from the Two Initial Conditions (IC) with the Two Force Fieldsa AMBER parm99 dihedral angle φ ψ a

AMBER parm96

IC 1

IC 2

total

IC 1

IC 2

total

26 59

32 77

58 136

43 491

47 502

90 993

Total numbers of rotations from both initial conditions for 40 ns are also listed with each force field.


Okumura and Okamoto

Figure 4. Time series of the backbone dihedral angle φ with the AMBER parm99 force field and initial condition 1. See the caption of Figure 2 for further details.

Figure 5. Time series of the backbone dihedral angle ψ with the AMBER parm99 force field and initial condition 1. See the caption of Figure 2 for further details.

only around ω ) 180°. We did not observe the stereo inversion of the R-carbon, either. The alanine dipeptide was always in the L-form during the MUBATH simulation. Much higher energy is required to realize the cis-trans inversions around the amide bonds or the chiral inversion between the L-form and the D-form. We now compare the two force fields for the alanine dipeptide. Figure 8 shows P(φ,ψ) obtained from the MUBATH MD simulations by the reweighting techniques at T ) 298 K and P ) 0.1 MPa. These distributions are the average of the two initial-condition results. The peaks of the

RR and RP states are higher than those of the PII and C5 states with the AMBER parm99 force field. On the other hand, the AMBER parm96 has high peaks of the PII and C5 states and low peaks of the RR and RP states. In the case of longer peptides or proteins, the RR state corresponds to an R-helix structure and the PII and C5 states correspond to a β-strand structure. It is known that, in general, the AMBER parm99 force field tends to form an R-helix structure and the AMBER parm96 force field tends to form a β-strand structure.31 The distributions P(φ,ψ) in Figure 8 are consistent with this feature. Figure 8a′,b′ show the contour maps of log P(φ,ψ).



Figure 6. Time series of the backbone dihedral angle ψ with the AMBER parm99 force field and initial condition 2. See the caption of Figure 2 for further details.

TABLE 2: Peak Positions (O, ψ) of the Probability Distribution P(O,ψ) at T ) 298 K and P ) 0.1 MPaa state

AMBER parm99

AMBER parm96

PII C5 RR RP RL Cax 7

(-65°, 155°) (-145°, 165°) (-65°, -25°) (-145°, 5°) (55°, 35°)

(-65°, 145°) (-145°, 155°) (-65°, -45°) (-145°, -65°) (35°, 85°) (45°, -125°)

a The results were obtained by the reweighting techniques from the MUBATH MD simulation. The state labels are taken from ref 37.

The peak position of each state is slightly different depending on the force field. It is also shown in Table 2. The Cax 7 state was sampled only with the AMBER parm96 force field. Figure 8a′ shows that the free energy barriers with the AMBER parm99 force field around φ ) 0° between the PII and RL states and around ψ ) -110 to -120° between the PII and RR states and between the C5 and RP states are higher than the respective barriers with the AMBER parm96 force field. This feature is consistent with the fact that the dihedral angles φ and ψ rotate with the AMBER parm96 force field better than with the AMBER parm99 force field as listed in Table 1. The volume under the surface P(φ,ψ) around each peak corresponds to the population W of each state. To calculate W, the whole (φ, ψ) plane was divided into six states as shown in Figure 8a′,b′ and Table 3. For example, the population WPII of the PII state is calculated by the integral of P(φ,ψ) in the area in which φ and ψ take the PII conformation:

WPII )

∫(φ,ψ)∈P

II

dφ dψ P(φ,ψ)

(33)

where the integration range of (φ, ψ) means the range for the corresponding state in Table 3. The population of each state at T ) 298 K and P ) 0.1 MPa is listed in Table 4. The statistical uncertainties were estimated by using the jackknife method51 in which the production run was divided into ten segments.

Figure 9 shows the ratio of each state population against the PII state population as a function of the inverse of temperature 1/T at the constant pressure of P ) 0.1 MPa. As the temperature increases, the relative population WC5/WPII of the C5 state increases with both force fields, as shown in Figure 9a,a′. On the other hand, population ratios WRR/WPII and WRP/ WPII for the RR and RP states decrease with the AMBER parm99 force field (Figure 9b,c), whereas they increase with the AMBER parm96 force field (Figure 9b′,c′). Such a difference due to the force fields is interpreted as follows. The distribution P(φ,ψ) must be flat in the limit of high temperature (T f ∞). As the temperature decreases, the RR and RP states get populated with the AMBER parm99 force field, but they do not get populated so much as the PII state with the AMBER parm96 force field so that the distributions P(φ,ψ) in Figure 8 can be obtained. This is why the population ratios of the RR and RP states against the PII state is different between the two force fields. The population ratio WRL/WPII of the RL state increases as the temperature increases with the AMBER parm99 force field (Figure 9d), although the error in WRL/WPII is too large to discuss its temperature dependence with the AMBER parm96 force field (Figure 9d′). This is because the population of the RL state was originally low so that it may not be estimated accurately. The population ratio of the Cax 7 state increases as the temperature increases with the AMBER parm96 force field (Figure 9e′). From thermodynamics, the increase in temperature at a constant pressure causes, in general, the increase of enthalpy. The increases in the population ratios W/WPII of some state against the PII state by the temperature increase indicates that enthalpy for that state is higher than that of the PII state. The difference in partial molar enthalpy ∆H of the C5 state from that of the PII state, for example, is calculated from the derivative of log (WC5/WPII) with respect to 1/T by

[

∆H ) -R

∂ log(WC5/WPII) ∂(1/T)

]

P

(34)

where R is the gas constant: R ) 8.3145 J/(mol · K). The derivative of log(WC5/WPII) was calculated here by the least-


Okumura and Okamoto

Figure 7. Probability distributions P(φ,ψ) of the backbone dihedral angles φ and ψ and the contour maps of log P(φ,ψ) with the AMBER parm99 force field at T ) 240 K and P ) 0.1 MPa. (a) P(φ,ψ) obtained with initial condition 1 by the conventional isobaric-isothermal MD simulation. (b) P(φ,ψ) obtained with initial condition 1 by the reweighting techniques from the results of the multibaric-multithermal MD simulation. (c) P(φ,ψ) obtained with initial condition 2 by the conventional isobaric-isothermal MD simulation. (d) P(φ,ψ) obtained with initial condition 2 by the reweighting techniques from the results of the multibaric-multithermal MD simulation.

Figure 8. Probability distributions P(φ,ψ) of the backbone dihedral angles φ and ψ at T ) 298 K and P ) 0.1 MPa, which were obtained by the reweighting techniques from the results of the multibaric-multithermal MD simulation (a) with the AMBER parm99 force field and (b) with the AMBER parm96 force field. (a′) and (b′) are the contour map of (a) and that of (b) in logarithmic scale, respectively.

squares fitting. The error bars were estimated by the jackknife method51 again in which the production run was divided into ten segments. The differences ∆H between the other states and the PII state were also obtained in the same way. The ∆H values are listed in Table 5. Table 5 also lists the experimental data by Raman spectroscopy experiment for the C5 and RR states.38 The RP, RL, and Cax 7 states were not observed by the Raman spectroscopy experiment. The partial molar enthalpy differences ∆H between

TABLE 3: Dihedral Angle Ranges of (O,ψ) AMBER parm99

AMBER parm96

state

φ

ψ

φ

ψ


(-100°, 0°) (130°, -100°) (-100°, 0°) (130°, -100°) (0°, 130°) (0°, 130°)

(90°, -110°) (90°, -110°) (-110°, 90°) (-110°, 90°) (-90°, 120°) (120°, -90°)

(-100°, 0°) (120°, -100°) (-100°, 0°) (120°, -100°) (0°, 120°) (0°, 120°)

(30°, -120°) (30°, -120°) (-120°, 30°) (-120°, 30°) (0°, 140°) (140°, 0°)


Figure 9. Population ratios W/WPII against the PII state as functions of the inverse of temperature 1/T at constant pressure of P ) 0.1 MPa, which was obtained by the reweighting techniques from the results of the multibaric-multithermal MD simulation of (a) and (a′) the C5 state, (b) and (b′) the RR state, (c) and (c′) the RP state, (d) and (d′) the RL state, and (e′) the Cax 7 state. (a)-(d) Population ratios obtained with the AMBER parm99 force field and (a′)-(e′) those obtained with the AMBER parm96 force field.

TABLE 4: Population W of Each State at T ) 298 K and P ) 0.1 MPa state

AMBER parm99

AMBER parm96


0.017(3) 0.046(6) 0.281(11) 0.648(11) 0.007(2) 0.000(0)

0.421(13) 0.478(11) 0.055(6) 0.044(4) 0.0011(9) 0.0016(7)

TABLE 5: Differences ∆H/(kJ/mol) in Partial Molar Enthalpy of the C5, rR, rP, rL, and Cax 7 States from That of the PII State Calculated by the MUBATH MD Simulationsa state

AMBER parm99

AMBER parm96

Raman

C5 RR RP RL Cax 7

3.6 ( 2.2 -1.5 ( 2.0 -3.4 ( 1.9 5.3 ( 3.3

1.8 ( 0.5 4.4 ( 1.1 7.5 ( 1.4 6 ( 31 25 ( 7

2.5 ( 0.3 4.4 ( 1.5

a

Raman experimental data are taken from ref 38.

the C5 state and the PII state are consistent in the two force fields and the Raman spectroscopy within their error bars. The ∆H value for the RR state with the AMBER parm96 force field agrees very well with the Raman spectroscopy data, whereas that with the AMBER parm99 force field was estimated too low. Moreover, both ∆H values for the C5 and RR states by the Raman spectroscopy lie between those for the AMBER parm99 and AMBER parm96 force fields.


Figure 10. Population ratios W/WPII against the PII state as functions of pressure P at constant temperature of T ) 298 K, which was obtained bythereweightingtechniquesfromtheresultsofthemultibaric-multithermal MD simulation. See the caption of Figure 9 for further details.

The RISM theory52–54 is a powerful theory to calculate thermodynamical quantities in aqueous solution. Takekiyo et al. have tried to calculate ∆H by the RISM theory.38 In this case, however, they did not report ∆H because of significant uncertainties in enthalpy. The multicanonical algorithm is another simulation technique to sample a wide range of conformational space, which is also frequently applied to biomolecular simulations. Because the volume is fixed in the multicanonical simulation, we cannot study the temperature dependence of the peptide structures under the same pressure as in the experimental conditions. Accordingly, ∆H cannot be calculated. The conventional isobaric-isothermal MD simulation can specify both temperature and pressure. However, it gets trapped in local-minimum free-energy states as shown in Figure 7. For this reason, the sampling efficiency is not sufficient for biomolecules. The MUBATH method has the advantages of both multicanonical algorithm and isobaric-isothermal algorithm. It can escape from local-minimum free-energy states and specify temperature and pressure. The MUBATH method, therefore, enables us to calculate accurate ∆H among several states of a peptide. As far as we know, this is the first computational work to calculate ∆H of a peptide. Figure 10 shows the population ratio of each state and the PII state as a function of P at the constant temperature of T ) 298 K. As the pressure increases, the population ratio of the C5 state against the PII state decreases with the AMBER parm99 force field, whereas it does not show clear pressure dependence with the AMBER parm96 force field, as shown in Figure 10a,a′. Relative population on the RR state decreases with the AMBER


Okumura and Okamoto

TABLE 6: Differences ∆V/(cm3/mol) in Partial Molar Volume of the C5, rR, rP, rL, and Cax 7 States from That of the PII State Calculated by the MUBATH MD Simulationsa state

AMBER parm99

AMBER parm96

Raman

RISM-KB

C5 RR RP RL Cax 7

1.5 ( 0.9 1.8 ( 0.8 0.6 ( 0.8 -5.2 ( 1.0

-0.3 ( 0.7 -2.3 ( 1.3 0.0 ( 0.6 -9 ( 7 1 ( 10

0.1 ( 0.3 1.1 ( 0.2

-1.8 -1.2

a

Raman experimental and RISM-KB data are taken from ref 38.

parm99 force field (Figure 10b), on the other hand, it increases with the AMBER parm96 force field (Figure 10b′). Although there are differences in the pressure dependence on these two states, the population ratio of the RP state against the PII state does not show clear pressure dependence with both force fields, as shown in Figure 10c,c′. The population ratio of the RL state against the PII state increases with both force fields, too (Figure 10d,d′). The Cax 7 state also does not show pressure dependence with the AMBER parm96 force field, as shown in Figure 10e′. A pressure increase at constant temperature generally causes a decrease in the volume. The decreases in the population ratio of some state and the PII state mean that the volume of that state is larger than that of the PII state. The difference in partial molar volume ∆V of the C5 state from that of the PII state, for example, is calculated from the derivative of log (WC5/WPII) with respect to P by

[

∆V ) -RT

∂ log(WC5/WPII) ∂P

]

T

(35)

The difference between the partial molar volume of the other states and that of the PII state was also obtained in the same way. The values of ∆V are listed in Table 6. As we described above on the pressure dependence of the population ratio, ∆V on the C5 and RR states are different in both force fields, whereas ∆V on the RP and RL states are consistent in the two force fields within their error bars. The differences in the two force-field parameters are the bond angle, torsion angle, and improper torsion-angle terms. The other terms of the bond length, van der Waals, and electrostatic terms are common. One might think that ∆V should be consistent in the two force fields because the molecular volume is often calculated by the van der Waals radius and it is common with the both force fields. The volume here, however, means the variation of each state population induced by the pressure increase. The population at several pressures depends on not only the van der Waals parameter but also all the other parameters. It is, therefore, reasonable for ∆V to be different with the two force fields. Table 6 also lists that the experimental data for the C5 and RR states by the Raman spectroscopy experiment by Takekiyo et al.38 They also calculated these values by the RISMKirkwood-Buff (KB) theory.38 In RISM-KB theory they used the AMBER parm99 force field for the peptide and the SPC/E force field for water molecules. The ∆V values on the RP, RL, and Cax 7 states by the Raman spectroscopy and the RISM-KB theory were not reported. Both Raman spectroscopy values of ∆V for the C5 state and that for the RR state are between the respective values with the AMBER parm99 force field and that with the AMBER parm96 force field. The ∆V value for the C5 state by the Raman spectroscopy is closer to ∆V with the AMBER parm96 force field than that with the AMBER parm99 force field. It is consistent with the AMBER parm96 value. On the other hand, ∆V between the RR state and the PII state by the

Raman spectroscopy agrees with that with the AMBER parm99 force field. Although both ∆V values by the RISM-KB theory for the C5 and RR states are different from those by the Raman spectroscopy, ∆V for the RR state by the RISM-KB theory agrees with that with the AMBER parm96 force field. As discussed above, the multicanonical simulation can escape from local-minimum free-energy states but cannot change volume. The conventional isobaric-isothermal MD simulation can change volume but is trapped in a local free energy minimum. Therefore, the partial molar volume cannot be calculated by these simulation algorithms. The MUBATH method has the advantages of both algorithms. It is an effective simulation technique to calculate ∆V and is, to our knowledge, the first simulational work to calculate the partial molar volume difference for a biomolecular system. 5. Conclusions We applied the multibaric-multithermal (MUBATH) MD algorithm to an alanine dipeptide in explicit water and showed that this algorithm is more effective than the conventional isobaric-isothermal MD algorithm. The MUBATH algorithm has an advantage that the simulation performs a two-dimensional random walk both in the potential-energy space and in the volume space so that it allows the simulation to escape from local-minimum free-energy states. It sampled the PII, C5, RR, RP, and RL state with both AMBER parm99 and AMBER parm96 force fields. With the AMBER parm96 force field, the Cax 7 state was also sampled. On the other hand, the conventional isobaric-isothermal simulation was trapped in local-minimum free-energy states and sampled not all of them. The MUBATH algorithm can provide accurate dihedral angle distributions P(φ,ψ). Other advantages of the MUBATH MD algorithm are that one can control temperature and pressure as in real experimental conditions and that one can obtain various isobaric-isothermal ensemble averages at different temperatures and pressures from only one simulation run. Utilizing these advantages, we calculated the temperature and pressure dependences of the population for each state. From these temperature and pressure dependences, the differences in the partial molar enthalpy ∆H and the partial molar volume ∆V were calculated between the PII state and other states. We also discussed the force-field dependencies of ∆H and ∆V using the AMBER parm99 and AMBER parm96 force field and compared our MUBATH results with the Raman spectroscopy data. The Raman spectroscopy data of ∆H for the C5 state agreed with both MUBATH data with the AMBER parm99 and AMBER parm96 force fields. The partial molar enthalpy difference ∆H for the RR state and the partial molar volume difference ∆V for the C5 state by the Raman spectroscopy were consistent with the respective data with the AMBER parm96 force field. On the other hand, the Raman spectroscopy data of ∆V for the RR state agreed with the AMBER parm99 result. Moreover, all the experimental data lie in between the corresponding simulation results with the two force fields. A better force field thus may be obtained by taking some intermediate parameters between the AMBER parm99 and AMBER parm96 force fields. The MUBATH algorithm is therefore a powerful simulation technique to calculate ∆H and ∆V. It overcomes the difficulty of the multicanonical method, in which a volume is fixed, and that of the isobaric-isothermal method, in which the simulation gets trapped in states of local free-energy minima. To study the temperature and pressure dependences of the protein structures, one has to investigate its free energy surface at several

Temperature/Pressure Dependence of Alanine Dipeptide values of temperature and pressure. The MUBATH algorithm will thus be of great use for the protein folding problem under various conditions of temperature and pressure. Not only for the protein folding problem but also for general problems in condensed matter physics, which deals with a complicated free energy surface such as phase transition, this algorithm will be useful. Acknowledgment. This work was supported, in part, by the Grants-in-Aid for the Next Generation Super Computing Project, Nanoscience Program and for Scientific Research in Priority Areas, “Molecular Theory for Real Systems”, from the Ministry of Education, Culture, Sports, Science and Technology, Japan. References and Notes (1) Nose´, S. Mol. Phys. 1984, 52, 255. (2) Nose´, S. J. Chem. Phys. 1984, 81, 511. (3) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (4) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (5) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (6) McDonald, I. R. Mol. Phys. 1972, 23, 41. (7) Hansmann, U. H. E.; Okamoto, Y. In Annual ReView of Computational Physics VI; Stauffer, D., Ed.; World Scientific: Singapore, 1999; p 129. (8) Mitsutake, A.; Sugita, Y.; Okamoto, Y. Biopolymers (Pept. Sci.) 2001, 60, 96. (9) Berg, B. A. Comput. Phys. Commun. 2002, 147, 52. (10) Berg, B. A.; Neuhaus, T. Phys. Lett. B 1991, 267, 249. (11) Berg, B. A.; Neuhaus, T. Phys. ReV. Lett. 1992, 68, 9. (12) Hansmann, U. H. E.; Okamoto, Y.; Eisenmenger, F. Chem. Phys. Lett. 1996, 259, 321. (13) Nakajima, N.; Nakamura, H.; Kidera, A. J. Phys. Chem. B 1997, 101, 817. (14) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635; 1989, 63, 1658(E). (15) Okumura, H.; Okamoto, Y. Chem. Phys. Lett. 2004, 383, 391. (16) Okumura, H.; Okamoto, Y. Phys. ReV. E 2004, 70, 026702. (17) Okumura, H.; Okamoto, Y. J. Phys. Soc. Jpn. 2004, 73, 3304. (18) Okumura, H.; Okamoto, Y. Chem. Phys. Lett. 2004, 391, 248. (19) Okumura, H.; Okamoto, Y. J. Comput. Chem. 2006, 27, 379. (20) Okumura, H.; Okamoto, Y. Bull. Chem. Soc. Jpn. 2007, 80, 1114. (21) Hukushima, K.; Nemoto, K. J. Phys. Soc. Jpn. 1996, 65, 1604. (22) Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 1999, 314, 141. (23) Nishikawa, T.; Ohtsuka, H.; Sugita, Y.; Mikami, M.; Okamoto, Y. Prog. Theor. Phys. Suppl. 2000, 138, 270. (24) Okabe, T.; Kawata, M.; Okamoto, Y.; Mikami, M. Chem. Phys. Lett. 2001, 335, 435.

J. Phys. Chem. B, Vol. 112, No. 38, 2008 12049 (25) Sugita, Y.; Okamoto, Y. In Lecture Notes in Computational Science and Engineering; Schlick, T., Gan, H. H., Eds.; Springer-Verlag: Berlin, 2002; pp 303-331; also see e-print: cond-mat/0102296. (26) Paschek, D.; Garcia, A. E. Phys. ReV. Lett. 2004, 93, 238105. (27) Paschek, D.; Gnanakaran, S.; Garcia, A. E. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6765. (28) Mitsutake, A.; Sugita, Y.; Okamoto, Y. J. Chem. Phys. 2003, 118, 6676. (29) Kollman, P. A.; Dixon, R.; Cornell, W.; Fox, T.; Chipot, C.; Pohorille, A. In Computer Simulation of Biomolecular Systems 3; Wilkinson, A., Weiner, P., van Gunsteren, W. F., Eds.; Elsevier: Dordrecht, 1997; p 83. (30) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. Phys. 2000, 21, 1049. (31) Yoda, T.; Sugita, Y.; Okamoto, Y. Chem. Phys. Lett. 2004, 386, 460. (32) Gould, I. R.; Kollman, P. A. J. Phys. Chem. 1992, 96, 9255. (33) Apostolakis, J.; Ferrara, P.; Caflisch, A. J. Chem. Phys. 1999, 110, 2099. (34) Smith, P. E. J. Chem. Phys. 1999, 111, 5568. (35) Mackerell, A. D., Jr.; Feig, M.; Brooks, C. L. J. Comput. Chem. 2004, 25, 1400. (36) Chekmarev, D. S.; Ishida, T.; Levy, R. M. J. Phys. Chem. B 2004, 108, 19487. (37) Chodera, J. D.; Swope, W. C.; Pitera, J. W.; Dill, K. A. Multiscale Model. Simul. 2006, 5, 1214. (38) Takekiyo, T.; Imai, T.; Kato, M.; Taniguchi, Y. Biopolymers 2004, 73, 283. (39) Bond, S. D.; Leimkuhler, B. J.; Laird, B. B. J. Comput. Phys. 1999, 151, 114. (40) Nose´, S. J. Phys. Soc. Jpn. 2001, 70, 75. (41) Okumura, H.; Itoh, S. G.; Okamoto, Y. J. Chem. Phys. 2007, 126, 084103. (42) Miller, T. F.; Eleftheriou, M.; Pattnaik, P.; Ndirango, A.; Newns, D.; Martyna, G. J. J. Chem. Phys. 2002, 116, 8649. (43) Yoshida, H. Phys. Lett. A 1990, 150, 262. (44) Berg, B. A.; Celik, T. Phys. ReV. Lett. 1992, 69, 2292. (45) Lee, J. Phys. ReV. Lett. 1993, 71, 211 and 2353(E). (46) Okamoto, Y.; Hansmann, U. H. E. J. Phys. Chem. 1995, 99, 11276. (47) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (48) Hansmann, U. H. E.; Wille, L. T. Phys. ReV. Lett. 2002, 88, 068105. (49) Wang, F.; Landau, D. P. Phys. ReV. Lett. 2001, 86, 2050. (50) Sayle, R. A.; Milner-White, E. J. Trends. Biochem. Sci. 1995, 20, 374. (51) Berg, B. A. Introduction to Monte Carlo Simulations and Their Statistical Analysis; World Scientific: Singapore, 2004. (52) Chandler, D.; Andersen, H. C. J. Chem. Phys. 1972, 57, 1930. (53) Hirata, F.; Rossky, P. J. Chem. Phys. Lett. 1981, 83, 329. (54) Perkyns, J. S.; Pettitt, B. M. Chem. Phys. Lett. 1992, 190, 626.

JP712109Q

Temperature and Pressure Dependence of Alanine Dipeptide Studied

Recommend Documents