Generalized-ensemble molecular dynamics and Monte Carlo algorithms beyond the limit of the multicanonical algorithm

Biomolecules such as proteins have complicated free energy surfaces with many local minima. Conventional molecular dynamics (MD) and Monte Carlo (MC) simulations in physical ensembles, such as the canonical [1–4] and isobaric–isothermal [5, 6] ensemble, tend to get trapped in these local-minimum states. One of the powerful techniques to avoid this difficulty is generalized-ensemble algorithms [7–11], such as the multicanonical algorithm [12–15]. In the multicanonical ensemble, a free one-dimensional random walk is realized in the potential-energy space and a simulation does not get trapped in free-energy-minimum states.

However, because the multicanonical simulation is performed in a fixed volume, neither the pressure dependence nor the temperature dependence at certain pressure can be investigated as in experiments. In order to overcome this difficulty, I proposed recently the multibaric–multithermal algorithm [16–22]. In this ensemble, two-dimensional random walks are realized both in the potential-energy space and in the volume space, so that the temperature and pressure dependence can be discussed.

Another problem of the multicanonical algorithm is that the determination of the multicanonical weight factor to give a flat distribution becomes difficult in a large system. In order to alleviate this difficulty, I proposed the partial multicanonical algorithm [23, 24]. The partial multicanonical simulation samples a wide range of a part of the potential energy terms, which is necessary to sample the conformational space widely. Thus, one can concentrate the effort to determine a multicanonical weight factor only on the important energy terms, so that the partial multicanonical algorithm has a higher sampling efficiency than the multicanonical and canonical algorithms.

In this review article, I explain the multibaric–multithermal algorithm and partial multicanonical algorithm. Applications of these algorithms to an alanine dipeptide in explicit water are also presented.

The multicanonical algorithm is briefly reviewed in section 2. In section 3, the multibaric–multithermal algorithm is explained. The partial multicanonical algorithm is described in section 4. Section 5 is devoted to conclusions.

In the canonical ensemble, the distribution of total potential energy E is given by

where n(E) is the density of states, β ₀ is the inverse of the product of the Boltzmann constant k ^_B and absolute temperature T₀, at which simulations are performed. This ensemble has a bell-shaped distribution in the potential-energy space.

In the multicanonical ensemble [12–15], every state is sampled with a weight factor W ^_muca (E) ≡ ^exp{−β ₀ E ^_muca (E)} so that a uniform distribution of E may be obtained,

where E ^_muca (E) and W ^_muca (E) are referred to as the multicanonical energy and the multicanonical weight factor, respectively. The difference between E ^_muca and E is written as δE(E):E ^_muca (E) ≡E+δE(E). The difference δE(E) is characteristic of the multicanonical simulation. The case of δE(E)=0 gives the regular canonical ensemble. The multicanonical simulation covers a wide range of E and escapes from local-minimum free-energy states. The expectation value of physical quantity A at temperature T is obtained by reweighting techniques [25] from

where 〈 ... 〉 ^_muca is the multicanonical ensemble average. Because of the random walks in the potential-energy space, we can calculate physical quantities in a wide range of T. The multicanonical weight factor W ^_muca (E) is not a priori known and has to be determined, for example, by the iterations of short simulations [26–28] and the Wang–Landau techniques [29].

In the isobaric–isothermal ensemble [5, 6], the distribution of potential energy E and volume V is given by

where n(E,V) is the density of states as a function of E and V and H is the 'enthalpy' (without the kinetic energy contributions): H=E+P₀ V. Here, P₀ is the pressure at which simulations are performed. This ensemble has bell-shaped distributions both in the potential-energy space and in the volume space, as shown in figure 1(a). In order to obtain the isobaric–isothermal ensemble, the combination of the Nosé thermostat [1, 2] and the Andersen barostat [5] is frequently employed.

Zoom In Zoom Out Reset image size

In the multibaric–multithermal ensemble [16–20], every state is sampled with a weight factor W ^_mbt (E,V) ≡ ^exp {−β ₀ H ^_mbt (E,V)} so that a uniform distribution of both E and V, as shown in figure 1(b), may be obtained,

Here, W ^_mbt (E,V) and H ^_mbt are referred to as the multibaric–multithermal weight factor and the multibaric–multithermal enthalpy, respectively. The difference between H ^_mbt and H is written as δH(E,V):H ^_mbt (E,V)=H+δH(E,V). The difference δH(E,V) is therefore characteristic of the multibaric–multithermal simulation. The case of δH(E,V)=0 gives the regular isobaric–isothermal ensemble.

The isobaric–isothermal Hamiltonian for a system consisting of N atoms based on the Nosé thermostat [1, 2] and the Andersen barostat [5] is given by

The variable is the coordinate scaled by the length of the simulation box: ( r _i is the real coordinate). We have introduced a simplified notation for the set of the scaled coordinates, . The variable is the conjugate momentum for . The real momentum p _i is related to the virtual momentum by , where s is the additional degree of freedom for the Nosé thermostat. The variable P_V and P_s are the conjugate momenta for V and s, respectively. The real-time interval Δt is associated with the virtual time interval Δt' by Δt=Δt'/s. The constant m_i is the mass of atom i. The constants W and Q are the artificial 'mass' related to V and s, respectively. In the case of a system consisting of Natoms, g equals 3N if the time development is performed in real time t, or g equals 3N+1 if the time development is performed in virtual time t'.

The equations of motion in real time are given from the Hamiltonian in equation (6) by

where the dot ( ⋅) stands for the real-time derivative ^d/^d t and F _i stands for the force acting on atom i. Performing the MD simulation by the equations of motion in equations (7)–(12), the multibaric–multithermal distribution in equation (5) is obtained.

In the case of δH(E,V)=∂ δH(E,V)/∂E=∂ δH(E,V)/∂V=0, the Hamiltonian in equation (6) and the equations of motion in equations (7)–(12) become those for the regular isobaric–isothermal MD simulation of the Nosé–Andersen formulation.

In order to perform a multibaric–multithermal MC simulation, we employ the Metropolis criterion on and V. The trial moves from to and is generated by uniform random numbers. The multibaric–multithermal enthalpy is consequently changed from to by the trial moves. The trial moves are now accepted with the probability

The multibaric–multithermal probability distribution is obtained by this scheme.

After an optimal weight factor W ^_mbt (E^, V) is determined—for example, by the iterations of short simulations [26–28] and by the Wang–Landau techniques [29]—a long production run is performed for data collection. The reweighting techniques [25] are used for the results of the production run to calculate the isobaric–isothermal-ensemble averages. The expectation value of a physical quantity A at T and P is then obtained by

where 〈 ... 〉 ^_mbt is the multibaric–multithermal 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.

I introduce here an application of the multibaric–multithermal MD algorithm to a system consisting of one alanine dipeptide molecule ((S)-2-(acetylamino)-N-methylpropanamide) and 63 water molecules [22]. The AMBER parm96 force field [30] was used for the alanine dipeptide molecule and the TIP3P rigid-body [31] model was employed for the water molecules. The electrostatic potential was calculated by the Ewald method. The time step was taken to be Δt=0.5 ^fs. The temperature was set as T₀=600 ^K and the pressure was set as P₀=300 ^MPa. The MD simulations were performed with the symplectic combination [32] of the Nosé–Poincaré thermostat [33, 34], the Andersen barostat [5] and the symplectic quaternion scheme [35]. In order to obtain an optimal multibaric–multithermal weight factor W ^_mbt (E,V), two-dimensional versions of both the iterative method [26–28] and the Wang–Landau techniques [29] for multicanonical weight-factor determination were employed. An optimal weight factor was obtained after these preliminary simulations in total for 2.61 ns in the AMBER parm96 simulation. Two sets of the initial values of the alanine-dipeptide backbone dihedral angles were prepared for the data-production runs. The initial condition 1 was ϕ=ψ=180° and the initial condition 2 was ϕ=180° and ψ=0°. The initial conformation 1 is shown in figure 2. We then performed a long multibaric–multithermal MD simulation for data collection. The simulation time was 20 ns for each initial condition. The isobaric–isothermal MD simulations of the same system were also performed for comparisons.

Zoom In Zoom Out Reset image size

Figure 3 shows the contour maps of the distributions  (ϕ, ψ) at T=240 ^K and P=0.1 ^MPa by the MD simulations started from the initial condition 1. The distribution  (ϕ, ψ) by the multibaric–multithermal MD were obtained using the reweighting techniques. The isobaric–isothermal distribution  (ϕ, ψ) has only four peaks of the α ^_R ,α ^_P ,P ^_II and C₅ states, as shown in figure 3(a). On the other hand, the multibaric–multithermal MD simulation provides correct distributions  (ϕ, ψ) with all the six peaks of the α ^_R ,α ^_P ,P ^_II ,C₅,α ^_L , and C₇ ^ax states, as shown in figure 3(b). These results mean that the multibaric–multithermal simulation can get over the free energy barriers and has a high sampling efficiency.

Zoom In Zoom Out Reset image size

The volume under the surface  (ϕ, ψ) around each peak corresponds to the population W of each state. The ratio of each state population W against the ^{P _II} state population was calculated as a function of the inverse of temperature 1/T at the constant pressure of P=0.1 ^MPa, as shown in figure 4(a) and as a function of the inverse of pressure P at the constant temperature of T=298 ^K, as shown in figure 4(b). The difference in partial molar enthalpy ΔH and the difference in partial molar volume ΔV can be obtained from the gradient in figure 4(a) and that in figure 4(b), respectively. Partial molar enthalpy and partial molar volume are important quantities in solution chemistry, which describes the change in the ratio of each state population due to temperature and pressure changes. As the temperature increases, the relative population W _{^C 5} /W ^{_{P II}} of the ^C ₅ state increases, as shown in figure 4(a). This indicates that the enthalpy for the ^C ₅ state is higher than that of the ^{P _II} state. The difference in partial molar enthalpy ΔH and the difference in partial molar volume ΔV of the ^C ₅ state from that of the ^{P _II} state, for example, are calculated from the derivative of W _{^C 5} /W ^{_{P II}} by

where R is the gas constant: R=8.3145 ^J (^mol K)⁻¹. The values of ΔH and ΔV are listed in tables 1 and 2, respectively. The differences ΔH and ΔV of the ^C ₅ state and ΔH of the α _R state are consistent in the multibaric–multithermal MD and the Raman spectroscopy experimental data [37]. The value of ΔV of the α _R state, however, is smaller than that of the Raman spectroscopy.

Zoom In Zoom Out Reset image size

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, ΔH and ΔV cannot be calculated by these simulation algorithms. The multibaric–multithermal algorithm has the advantages of both algorithms. This is the first simulational work to calculate ΔH and ΔV of biomolecules.

Although the multicanonical algorithm is a powerful technique, the determination of the multicanonical weight factor is time consuming for a large system. In order to alleviate the effort to determine the weight factor, the total potential energy is divided into two terms in the partial multicanonical algorithm [23, 24],

where E₁ is a part of the potential energy which is important to sample a wide range of the conformational space. The partial potential energy E₁ can be taken in any manner, in principle. The electrostatic, Lennard Jones and torsion-angle energy terms between solute atoms may be included in a biomolecular system. The second term E₂ is the rest of the potential energy. The partial multicanonical algorithm covers a wide range of E₁. The partial multicanonical simulation is performed with a weight factor W ^_pmuca (E₁,E₂) ≡ ^exp {−β ₀ E ^_pmuca (E₁,E₂)} so that a uniform distribution of E₁ may be obtained,

where E ^_pmuca (E₁,E₂) and W ^_pmuca (E₁,E₂) are referred to as the partial multicanonical energy and the partial multicanonical weight factor, respectively. The difference between E ^_pmuca and E depends only on E₁ and is written as δE₁(E₁),

The difference δE₁(E₁) is characteristic of the partial multicanonical simulation. If we define here a distribution of E₁ and the density of states n(E₁) of E₁ as

equation (17) is written as

The expectation value of a physical quantity A is obtained by the reweighting techniques [25] from

where 〈 ... 〉 ^_pmuca is the partial multicanonical ensemble average. Note that 〈A〉 ^_NVT can be calculated correctly at a temperature T only in the vicinity of T₀, although 〈A〉 ^_NVT can be calculated in a wide range of T in the multicanonical algorithm. This is because a limited range of E₂ is sampled with the normal Boltzmann weight factor ^exp(−β ₀ E₂) in partial multicanonical algorithm.

The scheme of the partial multicanonical algorithm includes the canonical and the multicanonical algorithms. In the case of E₁=0 and E₂=E, this algorithm becomes the regular canonical algorithm. In the case of E₁=E and E₂=0, the multicanonical algorithm is obtained.

In order to perform a partial multicanonical MD simulation, a constant-temperature MD method is employed. The partial multicanonical Hamiltonian based on the Nosé thermostat [1, 2] is given by

where p' _i is the conjugate momentum for r _i . The real momentum p _i is related to the virtual momentum p' _i by p _i = p' _i /s. In the Nosé–Hoover formalism [1–3], a variable is introduced and the equations of motion are given from the Hamiltonian in equation (23) by

Equation (25) is obtained from regular canonical-ensemble equations of motion by modifying force from

to

For a partial multicanonical MC simulation, we perform Metropolis sampling on the coordinates r . The trial moves from r _i to r' _i are generated by uniform random numbers. The partial multicanonical energy is consequently changed from E ^_pmuca ( r ) to E' ^_pmuca ≡E ^_pmuca ( r' ) by these trial moves. The trial moves are accepted by the Metropolis criterion [4] with the probability

The partial multicanonical probability distribution is obtained by this scheme.

The partial multicanonical MD algorithm was applied to a system consisting of one alanine dipeptide molecule and 67 water molecules [23]. The total potential energy is the sum of the potential energy between solute atoms E ^slt–slt , that of the solute–solvent interaction E ^slt–svn , and that between solvent atoms E ^svn–svn , as follows:

Each term for the AMBER force field [30] is given by

where the subscripts bond, angle, torsion, imp, elec and LJ stand for the bond-length, bond-angle, torsion-angle, improper torsion-angle, electrostatic and Lennard Jones potential energy, respectively. Although these terms may vary for different force fields, the partial multicanonical algorithm can be applied in essentially the same way. Even if we sample wide ranges of bond-length, bond-angle and improper torsion-angle energy in E ^slt–slt , the structure of the solute molecule is not expected to be changed very much. Sampling a wide range of the solvent energy E ^svn–svn would not change the solute structure very much, either. Sampling a wide range of the solute–solvent potential energy E ^slt–svn may sample a wide range of the interface area between the solute molecule and solvent molecules. The solute molecule of the alanine dipeptide in this simulation, however, is small and this effect was not taken into account in this simulation. Finally, the partial potential energy E₁ was here taken to be the sum of the torsion-angle, electrostatic and Lennard Johnes energy between the solute molecules as follows:

and E₂ was the rest of the potential energy,

Note that this choice of E₁ is not always the best. The best choice of E₁ should depend on the system and the force field.

In order to obtain an optimal partial multicanonical weight factor W ^_pmuca (E₁,E₂), both the iterative method [26–28] and the Wang-Landau techniques [29] were employed. An optimal weight factor was obtained after these preliminary simulations for 2 ns in total. I then performed a long partial multicanonical MD simulation for data collection for 6 ns at T₀=300 ^K.

The probability distributions of E and E₁ are shown in figure 5. The distribution of E in the multicanonical simulation was much wider than those in the other two simulations. The partial multicanonical algorithm sampled a much wider range of E₁. These generalized-ensemble MD simulations actually sampled a wide range of the potential energy, which and only which are predesignated to be sampled widely.

Zoom In Zoom Out Reset image size

The time series of ϕ is shown in figure 6. The canonical MD simulation did not make a complete rotation, as shown in figure 6(a). It covered only 65.6% of the whole range of ϕ. There seem to be potential energy barriers around ϕ=0° and ϕ=120°, which cannot be overcome by the canonical simulation (these energy barriers correspond to the boundaries between the local-minimum states). This means that the canonical simulation was trapped in local-minimum free-energy states. The multicanonical MD simulation also did not make a complete rotation and sampled 82.8% of the whole range of ϕ, as shown in figure 6(b). It overcame the potential energy barrier around ϕ=0° but did not overcome that around ϕ=120°. The partial multicanonical MD simulation, on the other hand, overcame both energy barriers and covered the whole range of ϕ, as shown in figure 6(c). The backbone dihedral angle ϕ underwent a complete rotation (covering 360°) six times (six turns) in this MD simulation for 6 ns.

Zoom In Zoom Out Reset image size

Note that these results do not mean that the multicanonical algorithm cannot sample the whole range of ϕ. If a better weight factor W ^_muca (E) were obtained with the longer preliminary simulations, the whole range of ϕ could be sampled also in the multicanonical ensemble. These results mean that the partial multicanonical algorithm shows higher sampling efficiency than the multicanonical algorithm if the weight factor is determined with the same effort.

The probability distributions (ϕ, ψ) of ϕ and ψ are shown in figure 7. In order to obtain the distribution (ϕ, ψ) in the canonical ensemble from the multicanonical simulation and that from the partial multicanonical simulation, the reweighting techniques in equation (3) and those in equation (22) were used, respectively. The probability distributions (ϕ, ψ) by the canonical simulation show four peaks of the ^{P _II} , ^C ₅, α ^_R , and α ^_P states, as shown in figure 7(a). However, they do not have any peaks in the range of 0° < ϕ < 120°. The distribution (ϕ, ψ) by the multicanonical MD simulation has five peaks of the ^{P _II} , ^C ₅, α ^_R , α ^_P and α ^_L states, as shown in figure 7(b). The partial multicanonical MD simulation also sampled the C₇ ^ax states in addition to the five peaks that were sampled in the multicanonical simulation, as shown in figure 7(c). In the present simulations, only the partial multicanonical simulation sampled all of the six states of the alanine dipeptide.

Zoom In Zoom Out Reset image size

I have reviewed two new generalized-ensemble MD and MC algorithms: the multibaric–multithermal and partial multicanonical algorithms. These generalized-ensemble algorithms allow the simulation to escape from local-minimum free-energy states and sample many conformations.

The multibaric–multithermal simulation performs a two-dimensional random walk both in the potential-energy space and in the volume space so that one can obtain various isobaric–isothermal ensemble averages at different temperatures and pressures from only one simulation run. From the temperature and pressure dependences of an alanine dipeptide in explicit water, the differences in the partial molar enthalpy ΔH and the partial molar volume ΔV were calculated between the ^{P _II} state and other states. This is the first work to calculate ΔH and ΔV by molecular dynamics simulations. The multibaric–multithermal algorithm will thus be a powerful simulation technique to study the temperature and pressure dependences of biomolecules like proteins.

In the partial multicanonical algorithm, the simulation covers a wide range of only important potential energy, which is necessary to investigate a wide range of the conformational space. Thus, one can concentrate the effort for the weight factor determination on the important energy terms. The MD simulation of an alanine dipeptide showed that the multicanonical simulation showed higher sampling efficiency than the canonical simulation. The partial multicanonical algorithm will thus be of great use for the protein folding problem. This idea, in which a simulation covers a wide range of only important potential energy terms, can be utilized not only for the multicanonical algorithm but also for multidimensional extensions of the multicanonical algorithm, such as the multibaric–multithermal [16–22] and multicanonical–multioverlap [38] algorithms.