Hamiltonian REM

In this program we adopt a variant of the replica exchange called Hamiltonian REM, that is far more flexible than the standard temperature REM technique illustrated above. In the Hamiltonian REM, each replica is characterized by a different potential energy rather than by a temperature. In its simplest implementation, the potential energies of the replicas differ by a scaling factor $c_{m}$, with $c_{1}=1$ for the target replica. Clearly, as long as the exchanged states differ only in the coordinates (i.e. momenta are not exchanged), the scaling of the potential energy of a canonical system (NVT) is equivalent to an inverse temperature scaling. Thus, Hamiltonian REM with full potential energy scaling and temperature REM are perfectly equivalent in an extended Monte Carlo simulation. When the replica simulations are done by numerically, integrating the Nosé-Hoover equations of motion at constant volume, Hamiltonian REM with full potential energy scaling and temperature REM are clearly no longer equivalent. Since momenta are not exchanged, Eq. 5.6 is valid for both full Hamiltonian REM and temperature REM, but in the latter technique both the kinetic and the potential energy are scaled, while in the former as implemented in ORAC one scales only the potential energy.

The advantage of using the Hamiltonian REM is two-fold: i) as all the replica have the same operating temperature, one does not have, like in temperature REM, to reinitialize the velocities after one successful configuration exchange and ii) since the mean atomic velocities are the same throughout the extended system, one does not have to adapt the time step size for preserving the quality of r-RESPA integrator as it should be done in temperature REM.

Hamiltonian REM can also be applied to a specific part of the potential, thereby localizing the effect of the configurational exchanges to specific part of the systems. Given a potential made up of a sum of various $i=1,...k$ contributions (e.g. stretching, bending, torsional, solute-solvent solute-solute solvent-solvent non bonded etc.), then one can define in a general way the $m$-th replica of the extended system as

$$V_m(X) = \sum_{i=1}^k c_i^{(m)} v_i(X)$$ (5.13)

where $c_i^{(m)}$ is the scaling factor for th $i$-th contribution, $v_i(X)$ of the potential of the $m$-the replica. So each replica is characterized by a $k$-dimensional scaling vector ${\bf c}_m = (c_{1}^{(m)},..c_{i}^{(m)},..c_{k}^{(m)})$ whose component are the scaling factors of $k$ contributions of the interaction potential for that replica. The target replica, replica 1, is such that $V_1(x)=V(x)$, the unscaled potential corresponding to the target system for which ${\bf c}_m = (1,1,1,1..)$. In vector notation we may compactly write Eq. 5.13 as

$$V_m(X)= {\bf c}_m \cdot {\bf v}(X)$$ (5.14)

Using this formalism, the probability of a configuration $X$ in the $m$-the replica may be written as

$$P_m(X)=\frac{1}{Z_m}e^{-\beta {\bf c}_m\cdot {\bf v}(X))}.$$ (5.15)

with $Z_m = \int e^{-\beta {\bf c}_m\cdot {\bf v}(X))} dX$. The detailed balance condition for the exchange of configurations between replica $m$ (characterized by the scaling vector ${\bf c}_m$) and replica $n$ (characterized by the scaling vector ${\bf c}_n$) is then given by

$$\frac{W(X,{\bf c}_m,X^{\prime},{\bf c}_n, \beta) }{W(X^{\prime},{... ...ta)}=e^{-\beta ({\bf c} _m-{\bf c} _n)\cdot [{\bf v} (X^{\prime})-{\bf v} (X)]}$$ (5.16)

Again, the detailed balance is implemented through Eq. 5.6 with

$$\Delta = \beta ({\bf c} _m-{\bf c} _n)\cdot [ ({\bf v} (X^{\prime})-{\bf v} (X)]$$

$$= \beta \sum_{i=1}^k (c^{(m)}_i- c^{(n)}_i)[V_i(X^{\prime}) - V_i(X)]$$ (5.17)

There is considerable freedom in the splitting of the potential (Eq. 5.13) and in the selection of the corresponding scaling factors. These factors are always positive and can be either smaller or greater than one, meaning that the corresponding potential contributions, for $m>1$, imply a heating and a cooling, respectively, of the involved degrees of freedom. For example we could use $c<1$ for torsions and and $c>1$ for bending, so that, with increasing $m$, torsional degrees of freedom are heated up while bending are frozen down.

Global scaling: In the present implementation of ORAC , one can do a global subdivision (i.e. ignoring the distinction between solvent and solute) of the overall atomistic interaction potential for biomolecular system according to the following:

$$V_m(X) = c_{\rm ba}^{m} (V_{\rm Bonds} + V_{\rm Angle} + V_{\rm i-tors}) + c_{t}^{(m)} ( V_{\rm tors} + V_{\rm 14}) +$$ (5.18)

$$+ c_{\rm nb}^{(m)} (V_{\rm vdw} + V_{qr} + V_{\rm qd})$$ (5.19)

where the meaning of the subscripts is given in Sec. 4.1. Typically one then sets $c_{\rm ba}^{(m)}=1~~\forall m$, as there is little advantage for conformational sampling in exchanging configurations involving stiff degrees of freedom such as bending, stretching and improper torsion. On the other hand, conformational transitions in proteins are mainly driven by torsional and intraprotein and protein-solvent non bonded interactions. It is thus convenient to heat up these degrees of freedom by scaling with $c_t^{(m)} < 1$ and $c_{\rm nb}^{(m)} < 1$ the corresponding potential functions. With this choice the quantity $\Delta$ in Eq. 5.17 is given by

$$\Delta = \beta (c_t^{(m)} - c_t^{(n)}) [ V_{\rm tors} (X^{\prime}) + V_{\rm 14}(X^{\prime}) - V_{\rm tors} (X)- V_{\rm 14}(X)] +$$

$$+ \beta (c_{\rm nb}^{(m)} - c_{\rm nb}^{(n)}) [ V_{\rm vdw}(X^{\pri... ...prime}) + V_{\rm qd}(X^{\prime}) - (V_{\rm vdw}(X) + V_{qr}(X) + V_{\rm qd}(X))$$ (5.20)

Local scaling: Hamiltonian REM in ORAC can work also by tempering only a user defined ``solute''. Unlike standard implementation of the solute tempering techniques[106], the ``solute'' in the present version can be any portion of the system including solvent molecules. Once the solute has been defined, the complementary ``non solute'' portion of the system is by definition the ``solvent''. In this manner, the scaling (i.e. the heating or freezing) can be localized in a specific part of the system with the remainder (the ``solvent'') of the system behaving normally (i.e. with the target interaction potential). In order to clarify how local scaling work, we illustrate the technique with a working general example. Suppose to choose a subset of atoms $n$ in the system that define the ``solute''. This subset can be chosen arbitrarily and may include disconnected portions of the protein, as well as selected solvent molecules. The solvent is then made up of the remaining $N-n$ atoms. According to this subdivision, the global potential of the system may be written as

$$V(X) = V^{\rm (Slt)}(X_n) + V^{\rm (Slt-Slv)} (X_n,X_{N-n}) + V^{\rm (Slv-Slv)} (X_{N-n})$$ (5.21)

where

$$V^{\rm (Slt)}(X_n) = V_{\rm tors}(X_n) + V_{\rm 14}(X_n) + V_{\rm vdw}(X_n) + V_{qd}(X_n)$$

$$V^{\rm (Slt-Slv)}(X_n,X_{N-n}) = V_{\rm vdw}(X_n,X_{N-n}) + V_{qd}(X_n,X_{N-n})$$

$$V^{\rm (Slv-Slv)} (X_{N-n}) = V_{\rm vdw}(X_n,X_{N-n}) + V_{qr}(X_n,X_{N-n}) + V_{\rm qr} + V_{\rm bd}(X_N)$$ (5.22)

The solute potential $V^{\rm (Slt)}(X_n)$ includes all the proper torsions and the 14 non bonded interactions involving the $n$ atoms of the solute.^5.2 The solute-solvent interaction involve all non-bonded interactions between the $N-n$ solvent atoms and the $n$ solute atoms. The solvent-solvent interaction involve all non-bonded interactions among the $N-n$ solvent atoms. As one can see, the global fast bonded potential $V_{\rm bd}(X_N) = (V_{\rm Bonds} + V_{\rm Angle} + V_{\rm i-tors})$ is assimilated to a solvent-solvent contribution. It should also be remarked that the reciprocal lattice contribution $V_{\rm qr}$, i.e. the long range electrostatics, is in any case assigned to the solvent-solvent term even if it includes all kinds of non bonded interaction (solute-solute, solvent-solute and solvent-solvent). The reason why $V_{\rm qr}$ is not split in the solute-solute, solute-solvent and solvent-solvent components is both physical and practical. Firstly, the long-range potential associated to each of three component of this term is expected, in general, to be rather insensitive along arbitrary reaction coordinates, such that a scaling of $V_{\rm qr}$ do not correspondingly produce a significant heating of any conformational coordinate. Secondly, in the Particle Mesh Approach approach the solute-solute, solvent-solute and solvent-solvent contribution to $V_{\rm qr}$ can no longer be easily separated, and this term must be thus arbitrarily assigned to one of the three components. Given the subdivision Eq. 5.22, the local scaling for replica $m$ in ORAC is implemented as

$$V_m(X) = c^{(m)}_{\rm Slt} V^{\rm (Slt)}(X_n) + c_{\rm (Slt-Slv)}^{(m)} V^{\rm (Slt-Slv)} (X_n,X_{N-n}) + V^{\rm (Slv-Slv)} (X_{N-n})$$ (5.23)

The solvent-solvent interactions, including the global bonded potential and the long-range electrostatic interactions, are not scaled in the local approach.

Solute-solute interactions and solute-solvent interactions as defined in Eq. 5.22 are scaled independently, thereby generalizing the so-called solute-tempering approach recently proposed by Liu et al.[106] This generality allows a complete freedom in the choice of the scaling protocol. For example, one can choose to set $c_{\rm (Slt-Slv)}^{(m)} > 1$, i.e. to progressively ``freeze'' the solute-solvent interaction as the replica index $m$ grows, while at the same time setting $c_{\rm Slt}^{(m)}=1$ for all replicas, thereby favoring, at large $c_{\rm (Slt-Slv)}^{(m)}$ the ``solvation'' of the solute, i.e., for example, favoring the unfolding.

The global REM algorithm (i.e. uniform scaling of the full interaction potential) as implemented in ORAC works also for constant pressure simulation (see ISOSTRESS directive). In that case, the selected external pressure pressure refers to that of the target replica ($m=1$). Since the $PV$ is a configurational term and is not scaled in the current implementation, the non target replicas sample coordinate configurations according to a higher external pressure, i.e. $P_m = P_1/c_m$ where $c_m$ is the scaling factor of replica $m$. This choice is done in order to avoid, through an increase of the external pressure, a catastrophic expansion of the simulation box for low scaling factors (or high temperatures).