Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2747743

Formats

Article sections

- Abstract
- I. Introduction
- II. Theory and methodology
- III. Results and discussion
- IV. Summary and conclusions
- References

Authors

Related links

J Phys Chem B. Author manuscript; available in PMC 2010 June 4.

Published in final edited form as:

PMCID: PMC2747743

NIHMSID: NIHMS117155

Department of Computational Biology, University of Pittsburgh School of Medicine, 3059 BST3, Pittsburgh, PA 15260

The publisher's final edited version of this article is available at J Phys Chem B

See other articles in PMC that cite the published article.

The loop 287–290 (Ile, Phe, Arg, and Phe) of the protein AcetylCholineEsterase (AChE) changes its structure upon interaction of AChE with dii*so*propylphosphorofluor*i*date (DFP). Reversible dissociation measurements suggest that the free energy (*F*) penalty for the loop displacement is Δ*F=F*_{free} − *F*_{bound} ~ −4 kcal/mol. Therefore, this loop has been the target of two studies by Olson’s group for testing the efficiency of procedures for calculating *F*. In this paper we test for the first time the performance of our “hypothetical scanning molecular dynamics” (HSMD) method and the validity of the related modeling for a loop with bulky sidechains in explicit water. Thus, we consider only atoms of the protein that are the closest to the loop (they constitute the “template”) where the rest of the atoms are ignored. The template’s atoms are fixed in the x-ray coordinates of the free protein, and the loop is capped with a sphere of TIP3P water molecules; also, the x-ray structure of the bound loop is attached to the free template. We carry out two *separate* MD simulations starting from the free and bound x-ray structures, where only the atoms of the loop and waters are allowed to move while the template-water and template-loop (AMBER) interactions are considered. The *absolute F*_{free} and *F*_{bound} (of the loop and water) are calculated from the corresponding trajectories. A main objective of this paper is to assess the reliability of this model and for that several template sizes are studied capped with 80–220 water molecules. We find that consistent results for the free energy (which also agree with the experimental data above) require a template larger than a minimal size, and a number of water molecules, which lead approximately to the experimental density of bulk water. For example, we obtain Δ*F*_{total} = Δ*F*_{water} +Δ*F*_{loop} = −3.1 ± 2.5 and −3.6 ± 4 kcal/mol for a template consisting of 944 atoms and a sphere containing 160 and 180 waters, respectively. Our calculations demonstrate the important contribution of water to the total free energy. Namely, for water densities close to the experimental value Δ*F*_{water} is always negative leading thereby to negative Δ*F*_{total} (while Δ*F*_{loop} is always positive). Also, the contribution of the water entropy *T*Δ*S*_{water} to Δ*F*_{total} is significant. Various aspects related to the efficiency of HSMD are tested and improved, and plans for future studies are discussed.

Calculation of the entropy, *S* and (Helmholtz) free energy, *F* constitutes a central problem in computer simulation, in spite of the significant progress achieved in the last 50 years.^{1}^{–}^{10} This problem is in particular severe in structural biology due to the flexibility and strong long-range interactions characterizing bio-macromolecules such as proteins. Thus, the potential energy surface of a protein, *E*(**x**) is rugged (**x** is the 3*N*-dimensional vector of the Cartesian coordinates of the molecule’s *N* atoms), i.e., this surface is “decorated” by a tremendous number of localized wells and “wider” ones, defined over regions, Ω_{m} (called microstates) – each consisting of many localized wells (an example for a microstate is the α-helical region of a peptide (see further discussion in our Refs. ^{11} and ^{12}). A microstate Ω_{m}, which typically constitutes only a tiny part of the entire conformational space Ω, can be represented by a sample (trajectory) generated by a *local* molecular dynamics (MD)^{13}^{,}^{14} simulation starting from a structure belonging to Ω_{m}. MD studies have shown that a molecule will visit a localized well only for a very short time [several femtoseconds (fs)] while staying for a much longer time within a microstate,^{15}^{,}^{16} meaning that the microstates are of a greater physical significance than the localized wells. Typically one is interested in calculating the free energy of the most stable microstates, rather than calculating the total free energy (of Ω). Identifying these microstates, in particular that with the global free energy minimum, is the extremely daunting task of protein folding. Notice that the microstates discussed above are meta-stable states; however, non-stable microstates, such as a transition state, might also be of interest.

Differences Δ*F* (Δ*S*) are commonly calculated by thermodynamic integration (TI) over physical quantities such as the energy, temperature, and the specific heat, (“calorimetric TI”) as well as non-thermodynamic parameters [free energy perturbation (FEP) is also included in this category].^{1}^{–}^{10} This is a robust and highly versatile approach, which enables one calculating a *small* difference in the binding *F* of two ligands **a** to **b** in the active site of a *large* enzyme solvated by water. However, while the mutation process leading from **a** to **b** is well controlled by TI, conformational changes in the entire protein (i.e.,”jumps” of side chains among rotamers) occur constantly and therefore the results might converge only after extremely long simulation times. Also, it is sometimes difficult to control the size and shape of the active site after mutation and the correct position of **b** in it. In many cases one is interested in calculating Δ*F _{mn}* between two microstates Ω

The *absolute S* and *F* can be calculated by harmonic^{17}^{–}^{19} and quasi-harmonic^{20}^{–}^{22}^{,}^{10} approximations and more general methods that are not limited to harmonic conditions, such as the local states (LS)^{23}^{–}^{25} and hypothetical scanning (HS)^{26}^{–}^{28} methods of Meirovitch, and other techniques.^{29}^{,}^{30} However, all these methods are not applicable as yet to diffusive systems such as explicit water. Notice that the absolute *F* can also be obtained with TI provided that a reference state *R* with known *F* is available and an efficient integration path *R*→*m* can be defined. A classic example is the calculation of *F* of liquid argon or water by integrating the free energy from an ideal gas reference state, where for such homogeneous systems TI constitutes a very efficient method. However, for non-homogeneous systems such integration might not be trivial, and in models of peptides and proteins defining adequate reference states and integration paths is a standing problem.^{7}

In recent years we have developed a new method for calculating the absolute *S* and *F* from a single sample called the hypothetical scanning Monte Carlo (HSMC) (or HSMD, where MD is used).^{11}^{,}^{12}^{,}^{31}^{–}^{36} HSMC(D) is based on ideas of the LS and HS methods mentioned above. Namely, in all these methods each conformation *i* of a sample (generated by MC MD or any other technique) is *reconstructed* step-by-step (from nothing, see II.4) using transition probabilities (TPs). The product of these TPs leads to an approximation for the correct Boltzmann probability *P _{i}*

HSMC(D) has been developed systematically as applied to liquid argon, TIP3P water,^{31}^{,}^{32} self-avoiding walks on a square lattice,^{33} and peptides,^{34}^{–}^{36} where for the first three models HSMC(D) results have been found to agree within error bars to TI results obtained by extensive MC or MD simulations. Also, for polyglycine molecules differences Δ*F _{mn}* and Δ

HSMD has also been applied to a mobile loop of the protein α-amylase,^{11}^{–}^{12} where the system is modeled by the AMBER96 force field^{37} alone and by AMBER96 with the GB/SA implicit solvent of Still and coworkers.^{38} In this work only protein atoms close to the loop (the template) were considered; however, they were kept fixed in their x-ray crystallographic positions, while only the loop’s atoms were moved by MD. A further development step of HSMD has been achieved recently, where the implicit solvent was replaced by explicit solvent, i.e., the α-amylase loop was capped with 70 TIP3P^{39} water molecules, which (together with the loop’s atoms) were moved in the MD simulation. Because the application of HSMD to water has not been optimized yet, the contribution of water to the free energy was calculated by a TI procedure incorporated within the framework of HSMD – this procedure is called HSMD-TI; in this TI process the water-loop interactions are gradually decreased to zero. Notice that the related statistical errors are relatively small because the loop structure is held fixed during the integration and only the water molecules are moved by MD. Previous studies of ${N}_{\text{water}}^{\text{min}}$, the minimal number of water molecules required for obtaining stable structures of proteins^{40} and loops^{41} (based on a root mean square deviation criterion) suggested that in some cases ${N}_{\text{water}}^{\text{min}}$ might be small (even 5 or 10 for certain loops). However, for loops, no systematic study has been carried out to determine the smallest template size and the minimal number of water molecules that would lead to consistent free energy results, i.e. results that are unchanged for larger template/water systems.

This study is carried out here, where an additional goal is to optimize HSMD(TI) further. Thus, HSMD(TI) is applied to the four-residue loop 287–290 ( Ile, Phe, Arg, and Phe) of the protein acetylcholinesterase (AChE) from *Torpedo California*. AChE degrades the neurotransmitter acetylcholine (ACh), producing choline and an acetate group. It is mainly found at neuromuscular junctions and cholinergic synapses in the central nervous system, where its activity serves to terminate synaptic transmission. AChE has a very high catalytic activity - each molecule degrades about 5000 ACh molecules per second. Reduction in the activity of the cholinergic neurons is a well-known feature of Alzheimer's disease. Thus, AChE inhibitors are employed to reduce the rate at which ACh is broken down, thereby increasing the concentration of ACh in the brain and combating the loss of ACh caused by the death of cholinergic neurons. AChE is also the target of many nerve gases, particularly organophosphates inhibitors (e.g. sarin), which block the function of AChE and thus cause excessive ACh to accumulate in the synaptic clefts. The excess ACh causes neuromuscular paralysis leading to death.^{42}

Of interest is the reaction of the inhibitor dii*so*propylphosphorofluor*i*date (DFP) with AChE^{43}^{–}^{48} which leads to a displacement of the loop’s backbone roughly by 4 Å. Moreover, comparison of experimental reversible dissociation constants measured for a variety of inhibitors of differing molecular size suggests that the free energy penalty for the loop displacement is on the order of 4 kcal/mol (i.e., *F*_{free} − *F*_{bound} ~ −4 kcal/mol).^{43}^{,}^{45}^{,}^{49}

The fact that the crystal structures of the free^{50} and bound^{43} enzyme are available, makes this loop a convenient target for testing free energy procedures, and indeed such tests were already carried out by Olson and collaborators.^{51}^{,}^{52} In the present study we shall also extend the scope of HSMD for an internal loop (the AchE loop is “hidden” in a cleft) with large sidechains; the earlier investigation of amylase has been applied to an external loop consisting of small sidechains.

Before describing the system and our procedures in detail, it should be stressed that our study relies heavily on the notion of a microstate introduced earlier. However, determining the *exact* limits of a microstate in conformational space is practically impossible and therefore it is commonly defined by an MC or MD sample initiated from a microstate’s structure. Thus, the microstate’s size typically increases with simulation time, *t* and estimation of *E*, *S*, and *F* depend on *t* as well. Therefore, estimation of Δ*E _{mn}*, Δ

As was pointed out above we study the 4-residue loop, Ile, Phe, Arg, Phe, (287–290) of AChE in two microstates related to the free and bound loop structures. The starting point is the available crystal structures of AChE from the protein data bank (PDB), where the unbound (free) structure 2ace^{50} is based on residues 4 to 535 and 2dfp, the structure of AChE with DFP consists of residues 2 to 535 (the bound structure).^{43} To be able to compare these structures 2dfp was trimmed so that both structures consist of the 531 residues, 4–535 and crystal water molecules were retained.

A close analysis shows that these conformations differ overall by (all heavy atoms) root mean square deviation (RMSD) of 1.22 Å, while the loop conformations differ by RMSD=1.38, 2.91, and 2.71 Å for the backbone, sidechains, and all loop heavy atoms, respectively, meaning that most of the deviation is due to sidechains movement. The RMSD of the templates, i.e. all atoms excluding the loop’s atoms, is 1.10 Å.

We used the program TINKER,^{53} and the AMBER force field,^{37} where Lys, Arg, and His are positively charged and Asp and Glu have a negative charge. Hydrogens were added to the free structure (including the crystal water) and the potential energy was minimized, with all heavy atoms restrained to their crystal positions by harmonic forces (with a force constant of 100 kcal/Å^{2}). Afterwards the loop and (TIP3P) water atoms were allowed to relax in the presence of a fixed template. The minimization eliminates bad atomic overlaps and strains in the original structures, while keeping the atoms reasonably close to the PDB coordinates.

Because the free and bound protein structures are similar, we adopt here the same strategy used in our previous studies,^{11}^{–}^{12} i.e., the structure of the bound loop (with hydrogens added) is attached to the free template. It would be impossible to compare the free energies of the free and bound microstates keeping the DFP attached to the bound template, because the two systems will be different. Thus, we calculate the free energy required to move the loop from the free to the bound microstate in the presence of the free template. In principle, one could carry out a similar study based on the bound template; however, many coordinates of the bound structure, 2dfp appear with large B-factors above 40, while the free structure (2ace) is much better resolved (B-factor values lower than 30). Therefore, the calculations were carried out with the free template. Also, taking into account the whole protein (of 8284 atoms) would be computationally prohibitive; therefore (as in Refs ^{11} and ^{12}), the template size is reduced to the *N*_{tmpl} atoms closest to the loop, where the rest of the atoms of the protein are ignored. More specifically (see Figure 1), the center of mass of the backbone atoms of the free loop is calculated as a (3D) reference point denoted **x**_{cmb} and a distance (*R*_{tmpl}) is chosen. If the distance of any atom of a residue from **x**_{cmb} is less than *R*_{tmpl}, the entire residue is included in the template; otherwise, the residue is eliminated. To assess the minimal template size needed we have studied *R*_{tmpl} =11, 12, and 13 Å corresponding to *N*_{tmpl}=800, 944, and 1100 protein atoms and 20, 30 and 40 crystal water molecules, respectively.

A two dimensional diagram of the template and the spherical water restraining region. The loop is represented as the heavy black curve, where the symbol, , denotes **x**_{cmb} - the center of mass of the loop backbone. **x**_{cmb} is the center of the dashed **...**

After defining the template, the water molecules and the orientations of the polar hydrogens of the free and bound loops (on the same template) are subjected to an optimization procedure. This procedure (where all the loop’s heavy atoms are kept fixed), is carried out in several rounds, each consisting of a 0.1 ns MD run (1 fs time step) at high temperature (1250 K) followed by energy minimization; the process is stopped after a fixed number of unsuccessful rounds (typically 100), namely rounds which do not reduce the energy further. Although the heavy atoms are fixed, considerable gain in potential energy is achieved. Next, keeping the template atoms fixed, the system energy is minimized where the loop and water molecules are allowed to move *freely*. The final “free” and “bound” structures constitute starting points for a further analysis.

While our considerations thus far are based on the crystal structures, we seek to simulate the loop in solution. Therefore, it is not clear whether the positions of the crystal waters are relevant for the solution environment. In particular, water molecules that are caged within the crystal structure are expected to stay there during the simulation, and thus can be considered as part of the template. Therefore, the number and arrangement of these waters should be globally optimized, practically by energy minimization, which, however, is a non-trivial problem (see below).

Our strategy is to add more water molecules to the already existing crystallographic waters. Thus, we define a sphere centered at **x**_{cmb} with a radius, *R*_{water} (*R*_{water}=*R*_{tmpl}+1 Å) where waters are added at random to the hemisphere oriented towards the exterior of the template. To hold these waters around the loop they are restrained with a flat-welled half-harmonic potential (with a force constant of 10 kcal mol^{−1}Å^{−2}) based on their distance from **x**_{cmb}. That is, if the distance of a water oxygen from **x**_{cmb} is greater than *R*_{water} a harmonic restoring force is applied, otherwise the restraining force is zero. In a previous paper^{41} the minimal number of waters, ${N}_{\text{water}}^{\text{min}}$ required to cap a loop has been studied with respect to the loop stability, i.e., the change in the loop’s RMSD from the initial (PDB) structure during 4 ns MD simulations. The objective of the present paper is similar but based on a free energy criterion; thus, we seek to find ${N}_{\text{water}}^{\text{min}}$ that for ${N}_{\text{water}}\ge {N}_{\text{water}}^{\text{min}},\phantom{\rule{thinmathspace}{0ex}}{F}_{\text{free}}-\phantom{\rule{thinmathspace}{0ex}}{F}_{\text{bound}}$ shows stability. To find ${N}_{\text{water}}^{\text{min}}$ we study 80 ≤ *N*_{water} ≤ 220 for different template sizes.

The fully hydrated system is then treated by a sequence of MD/minimization procedure similar to that outlined above, applied only to the hydrogen atoms of the loop and *all* water molecules restrained within the full sphere of radius *R*_{water}. Again, at the end of this optimization the energy is minimized where the loop’s atoms and the water molecules are free to move. The conformation of the free loop changed minimally from its crystal structure (RMSD=0.22, 0.63, and 0.57 Å for backbone, side chains and all loop atoms, respectively), while the “bound” conformation changed more (RMSD=1.16, 2.46, and 2.32 Å for backbone, side chains and all loop atoms, respectively; see discussion of the results in Table 1).

Minimum and maximum values of dihedral angles, α_{k}(min) and α_{k}(max) and their differences Δα_{k} (in degrees) for the free and bound samples^{a}

For each of the optimized “free” and “bound” structures (with several template sizes, *N*_{tmpl} and several levels of hydration defined by *N*_{water}) an MD run at 300 K is performed, where only the loop and water atoms are moved, while the template atoms are kept fixed. Thus, 1000 loop/water configurations are collected by retaining a configuration every 0.5 ps along the 0.5 ns MD trajectory. An equilibration run of 0.5 ns is carried out prior to the production run. The total potential energy *E*_{total} is the sum of partial energies related to the loop and water (the template-template energy is constant and thus is ignored),

$${E}_{\text{total}}=[{E}_{\text{loop-loop}}+{E}_{\text{loop-tmpl}}]+[{E}_{\text{water-water}}+{E}_{\text{water-tmpl}}+{E}_{\text{water-loop}}]={E}_{\text{loop}}+{E}_{\text{water}}$$

(1)

where *E*_{loop-loop} is the intra loop energy, *E*_{loop-tmpl} is the energy due to loop-template interactions; these energies define the total loop energy *E*_{loop}, and the interactions related to water are defined in a similar way, where their total is defined as *E*_{water}. The reconstruction of the loop structure is carried out in internal coordinates; therefore, the loop conformations simulated by MD are transferred from Cartesians to the dihedral angles (_{i}, ψ_{i}, and ω_{i} (*i*=1,*N*=4), the bond angles θ_{i,l} (*i*=1,*N, l*=1,3), the side chain angles χ, and the corresponding bond angles. For convenience, all these angles (ordered along the backbone) are denoted by α* _{k}*,

The partition function of the loop/water system is

$${Z}_{m}={\displaystyle \underset{m}{\int}\text{exp}-[E({\mathbf{x}}_{\text{loop}},{\mathbf{x}}^{N})/{k}_{B}T]d{\mathbf{x}}_{\text{loop}}d{\mathbf{x}}^{N}}$$

(2)

where *E*(**x**_{loop},**x**^{N})=*E*_{tot} defined in eq 1; **x**_{loop} is the Cartesian coordinates of the loop in microstate *m* (for brevity we omit the letter *m* in most equations). **x**^{N} is the 9*N* Cartesian coordinates of the water molecules, where for brevity we denote in the theoretical section *N*_{water} by *N* (*N*=*N*_{water}). However, it is convenient to change the variables of integration from **x**_{loop} to internal coordinates, which makes the integral dependent also on a Jacobian.^{17}^{,}^{18}^{,}^{20} This transformation is applied under the assumption that the potentials of the bond lengths (“the hard variables”) are strong and therefore their average values can be assigned to their corresponding Jacobian, *J*, which to a good approximation can be taken out of the integral (however, see a later discussion in this section). For the same reason one can carry out the integration over the bond lengths (assuming that they are not correlated with the α* _{k}*) and the remaining integral becomes a function of the

$$Z\prime =DZ=D{\displaystyle \underset{m}{\int}\text{exp}([-{E}_{\text{loop}}([{\mathrm{\alpha}}_{k}])-{E}_{\text{water}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})]/{k}_{B}T)d[{\mathrm{\alpha}}_{k}]d{\mathbf{x}}^{N}},$$

(3)

where [α_{k}] = [α_{1},…α_{K}] and *d*[α_{k}] = *d*α_{1}…*d*α_{K}. *D*(*T*) is a product of the integral over the bond lengths and their Jacobian, *J*. Due to the strong bond stretching potentials, *D* is *assumed* to be the same (i.e., constant) for different microstates of the same loop and therefore ln*D* cancels and can be ignored in calculations of free energy and entropy *differences*. The Jacobian of the bond angles, which should appear under the integral, is omitted because we have shown^{11} that it cancels out in entropy and free energy differences (our main interest). The Boltzmann probability density corresponding to *Z* (eq 3) is

$${\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})=\text{exp}\left\{-E([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})/{k}_{B}T\right\}/Z,$$

(4)

and the exact entropy *S* and exact free energy *F* (defined up to an additive constant) are

$$S=-{k}_{B}{\displaystyle \underset{m}{\int}{\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})\phantom{\rule{thinmathspace}{0ex}}\text{ln}{\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})d[{\mathrm{\alpha}}_{k}]d{\mathbf{x}}^{N}}$$

(5)

and

$$F={\displaystyle \underset{m}{\int}{\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})\{E([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})+{k}_{\mathrm{B}}T\phantom{\rule{thinmathspace}{0ex}}\text{ln}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})\}d[{\mathrm{\alpha}}_{k}]d{\mathbf{x}}^{N}}$$

(6)

It should be pointed out that the fluctuation of the *exact F* is zero,^{54}^{,}^{55} because by substituting ρ^{B}([α_{k}]) (eq 4) inside the curly brackets of eq 6 one obtains, *E*([α_{k}]) + *k*_{B}*T* ln ρ^{B}([α_{k}]) = −*kT* ln *Z* = *F*, i.e. the expression in the curly brackets is constant and equal to *F* for any set [α_{k}] within *m*. This means that the free energy can be obtained from *any single* conformation if its Boltzmann probability density is known. However, the fluctuation of an approximate free energy (i.e., which is based on an approximate probability density) is finite and it is expected to decrease as the approximation improves.^{54}^{,}^{55}^{,}^{30}^{–}^{33} Because HSMC(D) provides an approximation for ρ^{B}([α_{k}],**x**^{N}), it enables one, in principle, to estimate the free energy of the system from *any single* structure^{31}^{,}^{33} [Notice, however, that calculation of ρ^{B}([α_{k}],**x**^{N}) for a single conformation depends on the entire microstate as is also evident from the HSMC(D) procedure discussed later].

With MD the bond stretching energy is taken into account in eq 6 (and in free energy functionals defined later) while the corresponding entropy is ignored. The contribution of this energy to the free energy becomes an additive constant if one accepts the assumptions about the stretching energy and the corresponding Jacobian made prior to eq 6. This is a very good approximation; however, if the bond stretching entropy should be considered, we have argued in Ref. ^{11}, II.5 that it can be estimated approximately within the framework of HSMD.

HSMC(D) is based on the ideas of the *exact* scanning method where a system is constructed (from nothing) step-by-step using transition probabilities (TPs). The product of these TPs is equal to the Boltzmann probability (eq 4) from which the entropy and free energy can be calculated. Practically, a loop/water configuration is generated by initially building a loop structure followed by the construction of a configuration of the surrounding water molecules. In this way a sample of statistically independent system configurations can be obtained.

For simplicity this construction is described for a loop consisting of *M* Gly residues (with dihedral and bond angles denoted α_{k}, 1≤ α_{k} ≤ 6*M*=*K*)in microstate *m*; the loop is surrounded by *N*_{water} water molecules moving within the volume defined by the sphere of radius *R*_{water}, the template, and the loop. Starting from nothing, a conformation of the loop is built first by defining the angles α_{k} step-by-step using transition probabilities (TPs) and adding the related atoms.^{53} Thus, at step *k*, *k*−1 angles α_{1}, … ,α_{k−1}; these angles and the related structure (the past) are kept constant, and α_{k} is defined with the *exact* TP density ρ(α_{k}|α_{k−1} … α_{1}),

$$\rho {(\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1})={Z}_{\text{future}}({\mathrm{\alpha}}_{k},\cdots ,{\mathrm{\alpha}}_{1})/[{Z}_{\text{future}}({\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1})]$$

(7)

where *Z*_{future} (α_{k}, …,α_{1}) is a future partition function. The term “future” indicates that the integration defining *Z*_{future} is carried out over the variables α_{k+1},…,α_{K} and the 9*N* coordinates **x**^{N} of the water molecules which will be determined in future steps of the build-up process. In this integration the atoms treated in the past are held fixed in their coordinates (which are determined by α_{1} … α_{k}), while α_{k+1},…,α_{K} are varied in a restrictive way where the corresponding conformations of the “future” part of the loop remain in microstate *m*. Thus

$${Z}_{\text{future}}({\mathrm{\alpha}}_{k},\cdots ,{\mathrm{\alpha}}_{1})={\displaystyle \underset{m}{\int}\text{exp}-[(E({\mathrm{\alpha}}_{K},\cdots ,{\mathrm{\alpha}}_{1},{\mathbf{x}}^{N})/{k}_{\mathrm{B}}T]d{\mathrm{\alpha}}_{k+1}\cdots d{\mathrm{\alpha}}_{K}d{\mathbf{x}}^{N}}$$

(8)

where *E* (eq 1) is the total potential energy of the loop/template/water system, which also imposes the loop closure condition. The product of the TPs (eq 7) leads to the (Boltzmann) probability density of the entire loop conformation,

$${\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}({\mathrm{\alpha}}_{K},\cdots ,{\mathrm{\alpha}}_{1})={\displaystyle \prod _{k=1}^{K}\mathrm{\rho}({\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1}).}$$

(9)

After the loop structure has been constructed a configuration of water molecules is generated step-by-step, where the TP density ρ_{water}(**x**_{k}|α_{K},…,α_{1},**x**^{k−1}) for placing water molecule *k* at **x**_{k} is defined in a similar way to eq 7 based on *Z*_{future} ([α_{k}],**x**^{k}) and *Z*_{future}([α_{k}],**x**^{k−1}) where the loop conformation is kept constant and the *k*−1 water molecules that have already been treated are fixed at their coordinates, **x**^{k−1} and the summation in *Z*_{future}(**x**^{k}) is over the as yet undecided *N*−*k*+1 water molecules. (Notice that **x**_{k} denotes the 9 Cartesian coordinates of water molecule *k*, while **x**^{k} denotes the set of Cartesian coordinates of the *k* molecules 1,2,….,*k*). The Boltzmann probability density of the water is

$${\mathrm{\rho}}_{\text{water}}^{\mathrm{B}}({\mathrm{\alpha}}_{K},\cdots ,{\mathrm{\alpha}}_{1},{\mathbf{x}}^{N})={\displaystyle \prod _{k=1}^{N}{\mathrm{\rho}}_{\text{water}}({\mathbf{x}}_{k}|{\mathrm{\alpha}}_{K},\cdots ,{\mathrm{\alpha}}_{1},{\mathbf{x}}^{k-1})}$$

(10)

and the probability density ρ^{B}([α_{k}],**x**^{N}) of the loop/water configuration is the product of ${\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}])\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})$. One can define for *m* the loop entropy, *S*_{loop},

$${S}_{\text{loop}}=-{k}_{\mathrm{B}}{\displaystyle \underset{m}{\int}{\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}])\phantom{\rule{thinmathspace}{0ex}}\text{ln}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}])d[{\mathrm{\alpha}}_{k}]}$$

(11)

*S*_{loop} is defined up to an additive constant. Extending the exact scanning procedure to side chains is straightforward.

This construction procedure (which is not feasible for a large loop/water system) provides the theoretical basis for HSMC(D). Thus, the *exact* scanning method is equivalent to any other exact simulation technique (in particular, to Metropolis MC or MD) in the sense that large samples generated by such methods lead to the same averages and fluctuations. Therefore, one can assume that a given MC or MD sample has rather been generated by the exact scanning method, which enables one to reconstruct each conformation *i* by calculating the TP densities that *hypothetically* were used to create it step-by-step; this is the basis for HSMC(D) (as well as the HS and LS methods).

The theory of HSMD is again described as applied to a loop consisting of *M* Gly residues. One starts by generating an MD sample of microstate *m* with water molecules; the conformations are then represented in terms of the dihedral and bond angles α_{k},1≤α_{k}≤ 6*M*=*K*, and the variability range Δα_{k} is calculated,

$$\mathrm{\Delta}{\mathrm{\alpha}}_{k}={\mathrm{\alpha}}_{k}(\text{max})-{\mathrm{\alpha}}_{k}(\text{min}),$$

(12)

where α_{k}(max) and α_{k}(min) are the maximum and minimum values of α_{k} found in the sample, respectively. Δα_{k}, α_{k}(max), and α_{k}(min) enable one to verify that the sample has not “escaped” from microstate *m*.

System configuration ([α_{k}],**x**^{N}) (denoted *i* for brevity) is reconstructed in two stages, where the loop structure is reconstructed first followed by the reconstruction of the water configuration. Thus, at step *k* of stage 1, *k*−1 angles α_{k−1} … α_{1} have already been reconstructed and the TP density of α_{k}, ρ(α_{k}|α_{k−1},…,α_{1}) is calculated from an MD sample of *n _{f}* conformations (generated in Cartesian coordinates), where the

$${\mathrm{\rho}}_{\text{loop}}({\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1})\approx {\mathrm{\rho}}^{\text{HS}}({\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1})={n}_{\text{visit}}/[{n}_{f}\delta {\mathrm{\alpha}}_{k}]$$

(13)

where ρ^{HS}(α_{k}|α_{k−1}, ,α_{1}) becomes exact for very large *n _{f}*(

$${\mathrm{\rho}}^{\text{HS}}({\mathrm{\alpha}}_{k+n-1},\dots ,{\mathrm{\alpha}}_{k+1},{\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1})={n}_{\text{visit}}/[{n}_{f}{\displaystyle {\prod}_{j=k}^{j=k+n-1}\delta {\mathrm{\alpha}}_{j}],}$$

(14)

where in Ref. ^{11} we have shown that δα_{k} and δα_{k+1} can be optimized. Notice that with HSMD the future loop conformations generated by MD at each step *k* remain in general within the limits of *m*, which is represented by the analyzed MD sample. The corresponding probability density is

$${\mathrm{\rho}}^{\text{HS}}({\mathrm{\alpha}}_{K},\cdots ,{\mathrm{\alpha}}_{1})={\displaystyle \prod _{k=1}^{K-n+1}{\mathrm{\rho}}^{\text{HS}}({\mathrm{\alpha}}_{k+n-1},\dots ,{\mathrm{\alpha}}_{k+1},{\mathrm{\alpha}}_{k}|{\mathrm{\alpha}}_{k-1},\cdots ,{\mathrm{\alpha}}_{1}).}$$

(15)

ρ^{HS}([α_{k}]) defines an approximate entropy functional, denoted ${S}_{\text{loop}}^{\mathrm{A}}$, which can be shown (using Jensen’s inequality) to constitute a *rigorous* upper bound for *S*_{loop} (eq 11),^{26}^{,}^{30}

$${S}_{\text{loop}}^{\mathrm{A}}=-{k}_{\mathrm{B}}{\displaystyle \underset{m}{\int}{\mathrm{\rho}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}])\phantom{\rule{thinmathspace}{0ex}}\text{ln}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}^{\text{HS}}([{\mathrm{\alpha}}_{k}])d[{\mathrm{\alpha}}_{K}].}$$

(16)

${\mathrm{\rho}}_{\text{loop}}^{\mathrm{B}}$ (eq 9) is the Boltzmann probability density of [α_{K}] in *m*. Thus, for microstate *m*, ${S}_{\text{loop}}^{\mathrm{A}}$ can be estimated from a Boltzmann sample (of size *n*_{s}) generated by MD using the arithmetic average,

$${\overline{S}}_{\text{loop}}^{\mathrm{A}}(m)=-\frac{{k}_{\mathrm{B}}}{{n}_{\mathrm{s}}}{\displaystyle \sum _{t=1}^{{n}_{\mathrm{s}}}\text{ln}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}^{\text{HS}}(t,m)}$$

(17)

where ρ^{HS}(*t*,*m*) is the value of ρ^{HS}([α_{k}]) obtained for configuration *t* of the sample of *m*. ${\overline{S}}_{\text{loop}}^{\mathrm{A}}$ (with the bar) is an estimation of the ensemble average ${S}_{\text{loop}}^{\mathrm{A}}$ (eq 16); correspondingly, the ensemble averages of the energy are estimated from a sample of size *n*_{s} and should appear with a bar as well. However, from now on only estimations will be considered and for simplicity, all of them will appear without the bar, like the energies defined in eq 1. ${S}_{\text{loop}}^{\mathrm{A}}$ (eq 16 and eq 17) constitutes a measure for the loop flexibility of a pure geometrical character, i.e. with no *direct* dependence on the interaction energy. When the *converged* value of ${S}_{\text{loop}}^{\mathrm{A}}$ is considered, it will be denoted by *S*_{loop}, which is expected to be the exact value within the statistical error. In the same way, the difference in the loop entropies between two microstates obtained for a specific set of parameters is denoted by $\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}$ while the converged difference is denoted by Δ*S*_{loop}, thus

$$\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}={\overline{S}}_{\text{loop}}^{\mathrm{A}}(m)-{\overline{S}}_{\text{loop}}^{\mathrm{A}}(n)$$

(18a)

where

$$\mathrm{\Delta}{S}_{\text{loop}}={\overline{S}}_{\text{loop}}^{\mathrm{A}}(m)-{\overline{S}}_{\text{loop}}^{\mathrm{A}}(n)\phantom{\rule{thinmathspace}{0ex}}(\text{converged}).$$

(18b)

The difference in the loop energy between two microstates is (see eq 1),

$$\mathrm{\Delta}{E}_{\text{loop}}={E}_{\text{loop}}(m)-{E}_{\text{loop}}(n).$$

(19)

One can also define a free energy difference for the loop, Δ*F*_{loop},

$$\mathrm{\Delta}{F}_{\text{loop}}=\mathrm{\Delta}{E}_{\text{loop}}-T\mathrm{\Delta}{S}_{\text{loop}}$$

(20)

To reconstruct the water configuration one can use in principle the HSMC(D) procedure for fluids developed previously, which would lead to ${\mathrm{\rho}}_{\text{water}}^{\text{HS}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})$ [as an approximation for ${\mathrm{\rho}}_{\text{water}}^{\mathrm{B}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})$ (eq 10)] and then to the contribution of the water configuration to the free energy

$${F}_{\text{water}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})={E}_{\text{water}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N})+{k}_{\mathrm{B}}T\phantom{\rule{thinmathspace}{0ex}}\text{ln}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{\rho}}_{\text{water}}^{\text{HS}}([{\mathrm{\alpha}}_{k}],{\mathbf{x}}^{N}).$$

However, this procedure for fluids has not been optimized yet and it is relatively time consuming.

Alternatively, as in Ref. ^{12}, one can obtain *F*_{water} ([α_{k}],**x**^{N}) by a thermodynamic integration (TI) procedure based on the same reference state for the free and bound structures. In this state the water-water and water-template interactions are preserved but the (fixed) loop structure [α_{k} ] does not “see” the surrounding waters, i.e. the loop-water interactions (electrostatic and Lennard Jones) are switched off. These interactions are gradually increased (from zero) during an MD simulation of water (while the loop structure remains fixed at [α_{k}]). For [α_{k} ] of microstate *m* one obtains from the integration ${F}_{\text{water}}^{\text{TI}}([{\mathrm{\alpha}}_{k}],m)$ which is then averaged over the *n _{s}* sample configurations (as in eq 17). As described in the Appendix (V.1), the integration is carried out in two stages but in an

$${F}_{\text{water}}^{\text{TI}}(m)={F}_{\text{water}}^{\text{TI}}(\text{ch})+{F}_{\text{water}}^{\text{TI}}(\text{LJ})=\frac{1}{{n}_{\mathrm{s}}}{\displaystyle \sum _{t=1}^{{n}_{\mathrm{s}}}{F}_{\text{water}}^{\text{TI}}(\text{ch},t)+{F}_{\text{water}}^{\text{TI}}(\text{LJ},t)}$$

(21)

The difference in the free energy of water between *m* and *n* denoted Δ*F*_{water} is

$$\mathrm{\Delta}{F}_{\text{water}}={F}_{\text{water}}^{\text{TI}}(m)-{F}_{\text{water}}^{\text{TI}}(n)$$

(22)

Δ*E*_{water} (see eq 1) is

$$\mathrm{\Delta}{E}_{\text{water}}={E}_{\text{water}}(m)-{E}_{\text{water}}(n)$$

(23)

and

$$\mathrm{\Delta}{E}_{\text{total}}=\mathrm{\Delta}{E}_{\text{water}}+\mathrm{\Delta}{E}_{\text{loop}}$$

(24)

The total free energy is

$${F}_{\text{total}}(m)={F}_{\text{water}}^{\text{TI}}(m)+{E}_{\text{loop}}(m)-T{S}_{\text{loop}}(m),$$

(25)

and the difference in the total free energy between microstates *m* and *n* is

$$\mathrm{\Delta}{F}_{\text{total}}=\mathrm{\Delta}{E}_{\text{loop}}-T\mathrm{\Delta}{S}_{\text{loop}}+\mathrm{\Delta}{F}_{\text{water}}$$

(26)

The difference in the water entropy between *m* and *n* is

$$T\mathrm{\Delta}{S}_{\text{water}}=T[{S}_{\text{water}}^{\text{TI}}(m)-{S}_{\text{water}}^{\text{TI}}(n)]=\mathrm{\Delta}{E}_{\text{water}}-\mathrm{\Delta}{F}_{\text{water}},$$

(27)

where the corresponding difference in the total entropy is

$$T\mathrm{\Delta}{S}_{\text{total}}=T\mathrm{\Delta}{S}_{\text{water}}+T\mathrm{\Delta}{S}_{\text{loop}}.$$

(28)

Notice that all entropies and free energies are defined up to an additive constant.

The HSMD reconstruction procedure needs further discussion. Thus, the MD simulation of the future chain at step *k* starts from the reconstructed conformation *i*, and every *g*=10 fs the current conformation is retained, where the *n*_{init} initial retained conformations are discarded for equilibration. The next *n _{f}* (retained) future conformations are represented in internal coordinates and their contribution to

In previous work^{11}^{,}^{12}^{,}^{34}^{–}^{36} we developed procedures for keeping the loop in its original microstate and measures for estimating the extent of coverage of the microstate by the reconstructed samples (an important test for an adequate coverage is verifying that entropy differences are stable as *n _{f}* is increased.) In this paper these procedures have not been applied because the maximal

As has already been pointed out, our main interest is in the difference Δ*F _{mn}* (Δ

As stated above and earlier, an essential aim of this work is to find the minimal template size and minimal number of water molecules required for a reliable modeling of the loop behavior (i.e., that free energy differences for the minimal and larger models are approximately the same). We have tested three template sizes, defined by radii, *R*_{tmpl} =11,12, and 13 Å consisting of 783, 944, and 1141 atoms, respectively, where each template was studied with an increasing number of water molecules, *N*_{water}, ranging from 80 to 220.

For each pair (*R*_{tmpl},*N*_{water}) the solvated free and bound structures were initially optimized as described in sections II.1 and II.2 leading to two optimized structures. Then, starting from each optimized structure, an 1 ns MD trajectory was generated at 300 K, where the initial 0.5 ns part was used for equilibration and a sample of 1000 structures was generated from the last (half) part of the trajectory by retaining a structure to a sample every 0.5 ps. From these samples (which represent the free and bound microstates) the free energy and entropy are calculated.

These simulations and the reconstruction simulations (for generating the future samples) were carried out with the velocity-Verlet algorithm^{57} based on a time step of 2 fs, where bonds involving hydrogens (including those of water) were frozen to their ideal values by the RATTLE algorithm;^{57} the Berendsen^{57} heat bath controlled the temperature. Cut-offs on long-range interactions were not imposed, and in the reconstruction process a structure was added to the sample every *g*=10 fs, where the *n*_{init}=250 initial structures (2.5 ps) were discarded for equilibration. The future samples were generated for several bin sizes, where results are presented for δα_{k}=Δα_{k}/*l, l*=5, 10, 20, 30, 40, and 50, centered at α_{k} (i.e., α_{k} ±δα_{k}/2) (eqs 12–14). If the counts of the smallest bin are smaller than 50 the bin size is increased to the next size, and if necessary to the next one, etc. In the case of zero counts, *n*_{visit} is taken to be 1; however, an event of zero counts is very rare.

In Table 1 we present results for α_{k}(min), α_{k}(max) and Δα_{k} (eq 12) for the backbone dihedral angles and ψ and the sidechains, χ. These values are based on the samples of 1000 conformations generated for the free and bound microstates for the template of *R*_{tmpl}=12 Å and *N*_{water}=160. The table shows that the Δα_{k} values are relatively small (in most cases smaller than 80°) and they are comparable to those obtained by other samples [based on different (*R*_{tmpl},*N*_{water}) pairs]. This suggests that |Δ*S*_{loop}| (eq 18b) would be small as well (probably not larger than 2 kcal/mol; compare with Δα_{k} in Table 1 of Ref. ^{12}). These small Δα_{k} values stem from the constraints imposed by the template on the inner loop.

For comparison we also provide in Table 1 the dihedral values of the crystal structures, α_{k}(crystal). These angles enable one to determine whether the samples have escaped from their original microstates. While exact definition of a microstate is practically unfeasible (see discussion in section I.3 of Ref. ^{11}), we have accepted an “escape” criterion for a dihedral angle when α_{k}(crystal)+60° is smaller than α_{k}(max) or α_{k}(crystal)−60° is larger than α_{k}(min), i.e., if some angle values fall beyond the range α_{k}(crystal)±60°; these angles are bold-faced in the table. The table reveals that three and two backbone angles of the free and bound microstates, respectively have been escaped but the deviations are small; the corresponding sidechain (escaped) angles are five and three, among them are χ^{1} – χ^{4} of Arg(free) (where the deviation of χ^{1} is small) and the slightly deviating χ^{1} of Phe(bound). Thus, the original microstates are not completely retained (mainly for the sidechains) but this might be expected since we model the loops in solution rather in the crystal environment.

In Table 2 results are presented for *E*_{loop-temp}, *E*_{loop-loop}, *E*_{loop} and *E*_{total} (eq 1), and for ${F}_{\text{water}}^{\text{TI}}(\text{ch}),{F}_{\text{water}}^{\text{TI}}(\text{LJ}),\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{F}_{\text{water}}^{\text{TI}}$ (eq 21), calculated as described in the Appendix, V.1; we also provide results for ${F}_{\text{sum}}={E}_{\text{loop}}+{F}_{\text{water}}^{\text{TI}}$, which is the total free energy without the contribution of the loop entropy, *S*_{loop} (eqs 16 and 17). For each quantity, we have calculated the difference, Δ between the free and bound values. These results were obtained for *R*_{tmpl}= 12 and *R*_{water}= 13 Å for 80 ≤ *N*_{water}≤ 180. Similar results are presented in Table 3 for *R*_{tmpl}= 13 and *R*_{water}= 14 Å and 80≤*N*_{water}≤ 220.

Results in kcal/mol for energy and free energy components and their difference (Δ) between the free and bound microstates for *R*_{tmpl}=12 and *R*_{water}=13Å^{#}

Results in kcal/mol for energy and free energy components and their difference (Δ) between the free and bound microstates for *R*_{tmpl}=13 and *R*_{water}=14Å^{#}

To check the stability of these results, and assess their statistical errors they were calculated for an increasing sample size of *n*_{s}=10, 20 (not shown), *n*_{s}=40, and 80 for the higher *N*_{water} values. In some calculations the thermodynamic integration is doubled. The tables show that the results presented (in particular those for the larger *N*_{water}) are very stable. Statistical errors, *s*/(*n*_{s})^{1/2}, where *s* is the standard deviation, are smaller than 2.5 and 1.5 kcal/mol for *E*_{total} and the other quantities, respectively [because for *n*_{s}=40 and 80 the conformations are selected from the trajectory every 12.5 and 6.25 ps, respectively, the energy correlations are expected to be low. Notice that the correlations of the correct free energy are zero, because every conformation leads to the correct result (see II.3); for an approximate free energy the correlations will be smaller than the energy correlations.] The errors in Δ are differences between results obtained for the largest and smaller *n*_{s} values. The Δ results (in particular those for the larger *N*_{water}) are again very stable, i.e., the errors in most cases are relatively small, within 2 kcal/mol, due to cancelation of the individual errors in the difference; for example, ${F}_{\text{water}}^{\text{TI}}(\text{ch}),\phantom{\rule{thinmathspace}{0ex}}{F}_{\text{water}}^{\text{TI}}(\text{LJ})$, (and thus ${F}_{\text{water}}^{\text{TI}}$) (eq 21) change significantly in going from a single to a double integration, however, this change is comparable in the free and bound calculations and get cancelled in Δ. It should be pointed out that to verify our error estimation, for some of the larger values of *N*_{water} (*R*_{tmpl}=12 and 13 Å) several TI runs (for calculating ${F}_{\text{water}}^{\text{TI}}$) were carried out starting from different sets of velocities.

To estimate the minimum values of *R*_{tmpl} and *N*_{water} it is sufficient to study Δ*F*_{sum} - the difference in the sums of all contributions to the free energies excluding that of *S*_{loop} (eqs 16 and 17), where Δ*S*_{loop} is expected to be small. The results of Δ*F*_{sum} in Table 2 are positive, 31, 30 and 14 kcal/mol for the low density water, *N*_{water}=80, 100, and 120, respectively, meaning that the bound microstate is more stable (has lower free energy) than the free microstate (again, neglecting *S*_{loop}). On the other hand, as *N*_{water} is increased to 140, 160, and 180, Δ*F*_{sum} becomes negative, −8 −1 and −5 kcal/mol, respectively, i.e., the free loop becomes the most stable. We also provide a measure for the density of water, ρ_{water}=*N*_{water}/(hemisphere volume), i.e., ρ_{water}=*N*_{water}/[2π(*R*_{water})^{3}/3], which increases to 0.0304, 0.0348 and 0.0391 Å^{−3} for *N*_{water} 140, 160, and 180, respectively. These values are comparable to the experimental density of water, ρ_{water}=0.0350 Å^{−3}, which corresponds to 154 waters in the hemisphere; however, these densities are somewhat approximate since the hemisphere is not totally empty but contains part of the template. Thus, the density is lower for waters arranged initially in crystal water positions, or those which are put randomly inside the template. Also, bulk water can move during the simulation to crevices inside the template, and some might “seep” to the back of the template; however, because *R*_{water} - *R*_{tmpl},=1 Å (see figure 1) this “escape of waters is avoided to a large extent, as can be learnt from computer graphics.

The same behavior is shown in Table 3, where Δ*F*_{sum} is positive, 33, 21, 48, 14, and 14 kcal/mol for *N*_{water} = 80, 100, 120,140, and 160, respectively, becoming negative, Δ*F*_{sum}= −20, −18, and −3 kcal/mol for *N*_{water}=180, 200, and 220, i.e., for ρ_{water}= 0.0313, 0.0348, and 0.0383 Å^{−3} , respectively. Thus, as in Table 2, only ρ_{water} values close to the experimental density, 0.0350 Å^{−3} (corresponding to 193 waters in the hemisphere) lead to negative Δ*F*_{sum}, where the increase in *N*_{water} (as compared to *R*
_{tmpl}=12 Å) is due to the increase in *R*_{water} (from 13 to 14 Å).

We have also carried out calculations for *R*_{tmpl}= 11 and *R*_{water}=12 Å but have found Δ*F*_{sum} to be positive (~30 kcal/mol) for all values of *N*_{water}, even for ρ_{water} ~ 0.0350 Å^{−3} . This suggests that the template defined by *R*_{tmpl}= 11 Å (783 atoms) is too small. Indeed, computer graphics has shown that this template does not provide the required cover for the loop and as a result some dihedral angles, especially for Arg deviate significantly from their corresponding values in Table 1. On the other hand, the fact that for water densities close to the experimental value Δ*F*_{sum} becomes negative for both *R*_{tmpl}= 12 and 13 Å suggests that a template defined by *R*_{tmpl}≥ 12 Å (944 atoms) would be adequate. The strongly fluctuating (negative) values Δ*F*_{sum}= −20, −8, and −3 kcal/mol obtained for *R*tmpl= 13 Å and *N*_{water}≥ 180 probably reflect the difficulty to adequately optimize starting structures for MD simulations for these relatively large systems (see discussion in II.1 and II.2). Therefore, we calculate *S*_{loop} only for the smaller systems based on *R*_{tmpl}= 12 Å, where the (negative) Δ*F*_{sum} results are closer to each other.

Sometimes the energy rather than the free energy has been used in the literature as a criterion of stability; therefore, it is of interest to compare Δ*F*_{sum} to the corresponding values of Δ*E*_{total}. Table 2 and Table 3 show that these quantities (while not equal) are correlated, both decrease as *N*_{water} increases, with *R*^{2} = 0.74 and 0.79, respectively; see Figure, 2 and Figure, 3.

It should be pointed out that changes in *N*_{water} affect more the energy values of the free microstate than those of the bound one. Thus, in Table 2, *E*_{loop-tmpl} is changed within the ranges δ_{free}=102−52=50 and δ_{bound}=128−120=8 kcal/mol, where the corresponding ranges for *E*_{loop-loop} also differ significantly, δ_{free}=46−10=36 versus δ_{bound}= −10+2=12 kcal/mol. Similar relations (but somewhat less pronounced) are observed in Table 3, where the ranges for *E*_{loop-tmpl} are δ_{free}=91−68=23 and δ_{bound} =146−131=15, and for *E*_{loop-loop} are δ_{free}=49−20=29 and δ_{bound}=10−2=8 kcal/mol. The results for ${F}_{\text{water}}^{\text{TI}}(\text{ch})$ show the same tendency, where δ_{free}=84−63=21 and δ_{bound} =53−48=5 (Table 2), and δ_{free}=105−41=64 and δ_{bound} =62−47=15 kcal/mol (Table 3). On the other hand, ${F}_{\text{water}}^{\text{TI}}(\text{LJ})$ increases in most cases with *N*_{water} for both microstates.

It is of interest to determine which energy components lead to the change of Δ*F*_{sum} from positive to negative at *N*_{water} =140 (Table 2) and *N*_{water} =180 (Table 3). Thus, for each energy (and free energy) component in Table 2 we calculate two averages of Δ, one for the three lower values of *N*_{water} (80, 100, and 120) and the second for *N*_{water}=140, 160 and 180 (these two averages are denoted by Δ_{1}). The calculations show that two components contribute to the decrease of Δ*F*_{sum}, while one component contributes slightly to its increase. More specifically, Δ_{1}(*E*_{loop-tmpl}) decreases by 16 (from 60 to 44), ${\mathrm{\Delta}}_{1}({F}_{\text{water}}^{\text{TI}})$ decreases by 18 (from −6 to −24), while Δ_{1}(*E*_{loop-loop}) increases (slightly) by 4 (from −29 to −25 kcal/mol). In Table 3 the first group consists of *N*_{water}=80, 100, 120, 140 and 160 and the second group of *N*_{water}=180, 200 and 220. Unlike in Table 2, Δ_{1}(*E*_{loop-tmpl}) is practically the same for both groups, Δ_{1}(*E*_{loop-loop}) increases slightly (as in Table 2) by 3 from −31 to −28, while ${\mathrm{\Delta}}_{1}({F}_{\text{water}}^{\text{TI}})$ decreases significantly by 41 (from ~0 to −41 kcal/mol) and thus provides the sole contribution for the decrease of Δ*F*_{sum}. Thus, a consistent effect on the change of Δ*F*_{sum} is provided by the water component, ${F}_{\text{water}}^{\text{TI}}$. From a structural point of view, the results of Table 2 suggest that on average the loop moves (but not necessarily much) to decrease Δ_{1}(*E*_{loop-tmpl}) [where ${\mathrm{\Delta}}_{1}({F}_{\text{water}}^{\text{TI}})$ is decreased as well]. In table 3 the loop moves less, its energy is unchanged, but ${\mathrm{\Delta}}_{1}({F}_{\text{water}}^{\text{TI}})$, which consists of the loop-water interactions decreases significantly. It has been found difficult to verify this picture by a structural analysis. The consistent behavior of our model for the two templates (as a function of *N*_{water}) is reflected by the similar behavior of Δ*F*_{sum} for both templates (the Helmholtz free energy is the characteristic thermodynamic potential in the canonical ensemble and *F*_{sum} is very close the total Helmholtz free energy). The fact that some components of Δ*F*_{sum} behave differently is not unexpected, since the template size, the number of buried crystal waters in the template, and the optimization of the water-template system is different for *R*_{tmpl}=12 and 13 Å.

The discussion in the last two paragraphs demonstrates that *N*_{water} (which defines the water pressure) affects significantly the energy and free energy components of the free and bound microstates (in particular it leads to the change in the sign of Δ*F*_{sum}). Hence, one would expect that these energetic changes will correlate with the local water density around the loop. Thus, we calculated (for *R*_{tmpl}=12 Å) the average number of water molecules within spheres of radii 3, 4, 5, and 6 Å around the loops’ residues. However, for each *N*_{water} studied (*N*_{water}=120, 160 and 180) these numbers were found to be comparable for the free and bound microstates. On the other hand (as expected), the loop structures (in terms of dihedral angles) have been affected by *N*_{water} but not in a systematic way and therefore these changes are not discussed here.

Results for the loop entropy, ${S}_{\text{loop}}^{\mathrm{A}}$ (eq 17) appear in Table 4 for *R*_{tmpl}=12 and *R*_{water}=13 Å. Two sets of results are presented for *N*_{water}= 160 and 180 for the free and bound microstates and for their difference $T[{S}_{\text{loop}}^{\mathrm{A}}(\text{free})-{S}_{\text{loop}}^{\mathrm{A}}(\text{bound})]=T\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}$ (see the discussion preceding eq 18a). These results were obtained by reconstructing *n*_{s}=80 loop structures, distributed homogeneously along the entire sample of 1000 system configurations. The simulated future consists of the future part of the loop including all the surrounding water molecules. The results are presented for several values of *n _{f}* - the sample size of the future chains (eqs 13 and 14), where

HSMD results (in kcal/mol) for the loop entropy, $T{S}_{\text{loop}}^{\mathrm{A}}$ and for entropy differences $T{S}_{\text{loop}}^{\mathrm{A}}$ between the free and bound microstates at *T*=300^{a}

Being an upper bound, $T{S}_{\text{loop}}^{\mathrm{A}}$ (18a) is expected to decrease as the approximation improves, i.e., with decreasing the bin size – an expectation that is fully satisfied. For example, for *N*_{water}=160 and *n _{f}* =2000, the $T{S}_{\text{loop}}^{\mathrm{A}}(\text{free})$ values are 64.3, 61.2, 60.3, 60.1, 60.1, and 60.1 (kcal/mol + constant), i.e. they decrease for

This problem can be cured by dividing a (long) trajectory of size *n _{f}* into

The HSMD results for the entropy are compared to those obtained with the quasi-hermonic (QH) method from larger MD samples of 10,000 loop-water configurations, which are needed for achieving a reasonable precision. To avoid the “escape” of a sample from the original microstate, it consists of 10 separate samples of 1000 configurations (0.5 ns), each started from the same structure with a different set of initial velocities, where the initial trajectory of 0.5 ns is used for equilibration and is thus discarded. The central values of $T{S}_{\text{loop}}^{\text{QH}}$ (eq A3 in the Appendix) exceed the HSMD results for *TS*_{loop} (for *n*_{s}=2000) by ~13–16 kcal/mol. These elevated results are in accord with ${S}_{\text{loop}}^{\text{QH}}$ being an upper bound and are comparable to the overestimation of ${S}_{\text{loop}}^{\text{QH}}$ values found in our previous studies.^{11}^{,}^{12}^{,}^{34}^{–}^{36}

As stated above, we are mostly interested in the results for the difference in entropy between the free and bound microstates, *T*Δ*S*_{loop} (eq 18b). Table 4 shows that for *N*_{water}=160, the *converged* value is *T*Δ*S*_{loop}(*n* = 80) =1.2±0.1 kcal/mol which “covers” the $T\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}$ results (eq 18a) for Δα_{k}/*l, l* ≥30 for all *n _{f}* values, even for

For *N*_{water}=180, *T*Δ*S*_{loop}(*n*_{s} = 80) = −0.7±0.3 kcal/mol, representing the $T\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}$ results for Δα_{k}/*l, l* ≥30 for all *n _{f}* values, i.e. the bound microstate has the higher entropy. Here, the results for ${S}_{\text{loop}}^{\mathrm{A}}({n}_{\mathrm{s}}=40)$ for the free microstate are systematically higher than for ${S}_{\text{loop}}^{\mathrm{A}}({n}_{\mathrm{s}}=80)$ (as for

The computer time required to reconstruct a loop structure capped with *N*_{water}=160 and 180 is 8.1 and 9.2 hours CPU on a 2.1 GHz Athlon processor, meaning that the entire reconstructions required 1296 and 1472 h CPU. However, we have shown that considering only 10% (*n _{f}*=200) of the maximal reconstruction samples and using smaller samples of

In summary, the fact that *T*Δ*S*_{loop} changes sign in going from *N*_{water}=160 to 180 is not surprising since significant changes are also observed in other components of the energy in Table 2 (e.g., *E*_{loop-loop} = −10 and −37 kcal/mol for *N*_{water}=160 and 180, respectively). Also (like in previous studies), it is demonstrated that a limited future sampling in the reconstruction process (e.g., *n _{f}*=200) is sufficient for obtaining the correct

In Table 5 we summarize the contribution of the loop and water to the free energy averaged over samples of of *n*_{s}=80 configurations. We provide *E*_{water}, (eq 1) which includes the water-water, water-template and water-loop interactions, and the contribution of water to the free energy, ${F}_{\text{water}}^{\text{TI}}$ (eq 21); *TS*_{water} (which is not provided) can be obtained from *E*_{water} and ${F}_{\text{water}}^{\text{TI}}$. *E*_{loop}, which contains the loop-loop and loop-template interactions (eq 1), leads together with *TS*_{loop} (taken from Table 4) to *F*_{loop}. The entropy and free energy are defined up to additive constants, which are cancelled in the differences - our main interest.

Contributions (in kcal/mol) of the loop and water to the energy, entropy, and free energy and the differences of these values between the free and bound microstates^{#}

The table shows that for *N*_{water}=160, the *absolute* values of Δ*F*_{water} and Δ*F*_{loop} are comparable; however, the contribution of *T*Δ*S*_{water} to Δ*F*_{water} is significant, being ~33% of Δ*E*_{water} (in absolute values), while the contribution of *T*Δ*S*_{loop} to Δ*F*_{loop} is small where *T*Δ*S*_{loop} constitutes only ~3.6% of Δ*E*_{loop}. For *N*_{water}=180 the situation is even more extreme, where the contribution of *T*Δ*S*_{water} to Δ*F*_{water} (in absolute values) is larger (by 155%) than that of Δ*E*_{water} while the corresponding contribution of *T*Δ*S*_{loop} is again small (6.5%). In other words, for *N*_{water}=160, *E*_{water}(free) < *E*_{water}(bound) significantly and correspondingly also *TS*_{water}(free) < *TS*_{water}(bound) significantly. For *N*_{water}=180, *E*_{water}(free) is only slightly smaller that *E*_{water}(bound) while *TS*
_{water}(free) is *larger* than *TS*_{water}(bound). In any case, the effect of the entropy of water is significant. Also, the results of Table 5 and results in Table 2 and Table 3 demonstrate the important contribution of water to the total free energy, where for water densities close to the experimental value $\mathrm{\Delta}{F}_{\text{water}}^{\text{TI}}$ is always negative leading thereby to negative Δ*F*_{total} (while Δ*F*_{loop} or Δ*E*_{loop} are always positive; see also Table 6)

Total energy, entropy, and free energy (in kcal/mol) at *T*=300 K and their differences between the free and bound microstates

The total contributions of *E*_{total} (=*E*_{water}+*E*_{loop}) and *TS*_{total} (=*TS*_{water}+*TS*_{loop}) to *F*_{total} (=*F*_{water}+*F*_{loop}) are summarized in Table 6 (again, for the entropy only the results for *T*Δ*S*_{total} are given). For *N*_{water}=160, Δ*E*_{total} and *T*Δ*S*_{total} are comparable with the same sign (as for water above), and thus leading to a small negative Δ*F*_{total}= −3.1±2.5 kcal/mol. For *N*_{water}=180, Δ*F*_{total}= −3.6±4 kcal/mol is equal to the value obtained for *N*_{water}=160 within a larger error; however, Δ*F*_{total}(180) is based on *positive* Δ*E*_{total} and *T*Δ*S*_{total} values (i.e., both the energy and entropy of the free microstate are larger than their bound counterparts). The fact that the entropic effects are significant means that (as the table also demonstrates) Δ*E*_{total} by itself does not constitute an adequate criterion of stability. We also provide in the table the results for Δ*F*_{sum} from Table 2, which are very close to those of *F*_{total} due to the small contribution of *T*Δ*S*_{loop}, meaning that in these cases *F*_{sum} serves as a reliable measure of stability.

The above results for Δ*F*_{total} are equal to the experimental value, ~ −4 kcal/mol within the error bars. Furthermore, one would expect |*T*Δ*S*_{loop}(*N*_{water}=140)| to be small, similar to the results obtained for *N*_{water}=160 and 180 (this is based on values for Δα_{k} as those presented in Table 1). Therefore, Δ*F*_{total}(*N*_{water}=140) is approximately represented by Δ*F*_{sum}(*N*_{water}=140) = −8±2 kcal/mol (see Table 2), which is close to the experimental value; the same is expected for the calculations based on *N*_{water}=220 (for *R*_{tmpl}=13 and *R*_{water}=14 Å), where Δ*F*_{sum}(*N*_{water}=220) = −3±2 kcal/mol (Table 3). This agreement of our calculations with the experiment should be accepted with some caution because for *R*_{tmpl}=13 and *R*_{water}=14 Å the Δ*F*_{sum} values for *N*_{water}=180 and 200 (−20±2 and −18±2 kcal/mol, respectively; see Table 4) are too low to lead to Δ*F*_{total} values close to the experimental result. However, the significant difference between the Δ*F*_{sum} results for *N*_{water}=180, 200, and 220, suggests that the initial water configurations in the larger systems of *R*_{tmpl}=13 and *R*_{water}=14 Å have not sufficiently optimized (see II.2). We find this global energy optimization to be the most uncertain, while the accuracy of the simulation results (including TI) is adequate.

In the present paper HSMD-TI has been applied to the mobile loop 287–290 (Ile, Phe, Arg, and Phe) of the protein Acetylcholineesterase (AChE), where the difference in free energy between the free and bound structures of the loop, *F*_{free}−*F*_{bound} has been estimated experimentally to be ~ − 4 kcal/mol. In view of this result, the main objectives of the present paper have been related to modeling issues: are a fixed template and capped water constitute an adequate model? and if they are, what is the minimal template size and minimal number of water molecules required to obtain reliable results? Another objective has been to further improve the efficiency of various components of HSMD-TI, in particular due to the fact that HSMD is applied for the first time to an internal loop consisting of residues with large sidechains. To achieve these aims, we have carried out a systematic study consisting of several template sizes which are capped with 80–220 water molecules.

We have emphasized the difficulty to determine the optimal distribution of water. Thus, the hemisphere contains bulk water (which can be simulated adequately by MD) as-well- as *internal* water molecules that reside in crevices on the surface and inside the template; because some of these waters are practically not mobile (by MD), they can be considered as part of the template, and their number, spatial positions, and orientations affect significantly the system energy. Using crystal water as internal ones might not always be adequate as our objective is to model the system in solution rather than in the crystal. Alternatively, one can seek to optimize (prior to the production simulations) the distribution of internal waters by energy minimization which, however, is an extremely difficult task. We have adopted a strategy where the initial configuration of water consists of ~30 crystal water positions and water molecules distributed randomly in the hemisphere; this configuration is optimized by a procedure based on a series of high temperature MD runs followed by energy minimizations. Clearly, finding the *global* energy minimum is practically unfeasible, but this is a general modeling problem not specific to HSMD-TI.

We have found that minimizing the energy of water before each integration step (during the TI process for calculating the contribution of water to the free energy) improves the convergence of the TI procedure significantly. Also, we have shown that satisfactory accuracy can be obtained by reconstructing a relatively small number of system configurations (80 or even 40), provided that they are selected homogeneously from the entire sample. This leads to a reduction in computer time by a factor of 7 as compared to our previous studies. In the reconstruction of the sidechain of Arg, three successive χ angles were treated successfully in each step, suggesting that 4 successive backbone angles (i.e., two pairs of a dihedral and the following bond angle) could be treated as well, which would decrease computer time further.

The limited number of atoms in the loop and waters, and the (relatively small) *constant* template lead to relatively small statistical errors (which for an *N*-atom system increase as *N*^{1/2}). As in our previous studies (and in accord with theoretical arguments discussed in Ref. ^{11}), we find that systematic errors in ${S}_{\text{loop}}^{\mathrm{A}}(m)$ are cancelled to a large extent in differences, $\mathrm{\Delta}{S}_{\text{loop}}^{\mathrm{A}}$ (eq 18a), and a similar cancelation occurs for the free energy of water, Δ*F*_{water} and other energy components. It should be pointed out (that in agreement with our previous studies and Ref. ^{58}) the quasi-harmonic approximation has been found to overestimate the entropy significantly, which might reflect strong long-range correlations and an-harmonic effects within the loop due to the loop-template, loop-loop and loop-water interactions; also, the results for $\left|T\mathrm{\Delta}{S}_{\text{loop}}^{\text{QH}}\right|$ overestimate the HSMD values for |Δ*S*_{loop}|. Finally, notice that the calculations of the transition probabilities of different steps are completely independent and they are also independent of the integration of water. Therefore, the reconstruction steps and TI of water can be fully parallelized.

The main conclusions from the present study (besides the above points which are mostly of a technical character) is that our approximate model is reliable, at least for the loop studied. This model is based on the same constant template for the free and bound microstates, where the loop is capped with a sphere of water molecules. By studying several template sizes and an increasing number of water molecules, we have found that to obtain consistent results for the free energy, Δ*F*_{total}=*F*_{free}−*F*_{bound} the template should be larger than a minimal size, and the number of water molecules (in the hemisphere of *R*_{water} -*R*_{tmpl}=1 Å) should lead approximately to the experimental density of bulk water. Also, our results for Δ*F*_{total} agree with the experimental data, ~ −4 kcal/mol. Our results demonstrate the important contribution of water to the total free energy, where for water densities close to the experimental value Δ*F*_{water} is always negative leading thereby to negative Δ*F*_{total} (while Δ*F*_{loop} is always positive). Also, the contribution of the water entropy, *T*Δ*S*_{water} to Δ*F*_{total} is significant. The next step would be to apply HSMD-TI to the free and bound loops attached to their own templates (defined by the PDB structures) rather than to a common template as was done here.

This work was supported by NIH grant 2-R01 GM066090-4 A2.

As described earlier, in the TI process the interaction energy [electrostatic and Lennard Jones (LJ)] between a *fixed* loop structure and the (moving) water molecules is decreased gradually to zero (rather than increased from zero) at constant *T* and *V*, where the water-water and water-template potential energy is unchanged. For the (LJ) potential we have used the shifted scaling potential, introduced by Zacharias et *al*.,^{59}

$$\mathrm{\varphi}({r}_{ij},\mathrm{\lambda})=\mathrm{\lambda}4\mathrm{\epsilon}\left[\frac{{\mathrm{\sigma}}^{12}}{{\left({r}_{ij}^{2}+\text{\delta}(1-\mathrm{\lambda})\right)}^{6}}-\frac{{\mathrm{\sigma}}^{6}}{{\left({r}_{ij}^{2}+\delta (1-\mathrm{\lambda})\right)}^{3}}\right]\phantom{\rule{thinmathspace}{0ex}},$$

(A1)

where the shift parameter, δ=3 Å^{2}, prevents the divergence of the potential (and its derivative) at small pair separations; a similar scaling function is used for the Coulomb interactions. The free energy derivatives with respect to λ, *F*/λ is

$$\frac{\partial F}{\partial \mathrm{\lambda}}={\langle \frac{\partial E({\mathrm{x}}^{N},\mathrm{\lambda})}{\partial \mathrm{\lambda}}\rangle}_{\mathrm{\lambda}},$$

(A2)

where the derivative of the energy is calculated analytically. The integration with respect to λ is carried out by dividing the range [1,0] into 20 equal integration bins Δλ_{i}. The (λ=1→λ=0) integration of the electrostatic interactions (i.e. charge elimination) is carried out first (in the presence of intact LJ interactions) followed by a ε=λ =1→0 integration of the LJ interactions. Thus, the entire two-stage process is based on 40 *F*/λ_{i} integration steps.

The MD simulation consists of a 2 fs integration step. Each (Δλ_{i}) step starts with energy minimization (based on λ_{i}) of the last structure obtained in the simulation of the *i*-1 step, followed by 5 ps MD simulation for equilibration, which is discarded; the next 20 ps of MD simulation a configuration is retained every 0.02 ps, i.e., altogether 1000 configurations are used for evaluating <*F*/λ_{i}>. It should be pointed out that the energy minimization (which was not performed in paper 20) has contributed to a nice convergence of the integration results.

For a single loop structure the free energy integration requires approximately the same time for each procedure (electrostatic or LJ) i.e., ~2×9.2 =18.4 h CPU for *N*_{water}=160 and ~10.5×2=21 h CPU for *N*_{water}=180 on a 2.1 GHz Athlon processor; these times refer to the *double* TI, where the MD simulation at each step is doubled, i.e. it is based on 40 ps for each λ_{i} (see Table 2 and Table 3). These simulation lengths are adequate because the loop’s conformation is kept fixed (as well as the template) and only the water molecules are moved by MD during integration.

With the QH method introduced by Karplus and Kushick,^{20} the Boltzmann probability density of structures defining a microstate is approximated by a multivariate Gaussian. Thus,

$${S}_{\text{loop}}^{\text{QH}}(m)=({k}_{\mathrm{B}}/2)\{N+\text{ln}[{(2\mathrm{\pi})}^{N}\text{Det}(\mathbf{\sigma})]\}$$

(A3)

where the covariance matrix, **σ**, is obtained from a local MD (MC) sample and *N* is (usually) the number of internal coordinates. Clearly, *S*^{QH} constitutes an upper bound for *S* since correlations higher than quadratic are neglected; also, an-harmonic contributions are ignored, and QH is not suitable for diffusive systems such as water. While QH has been used extensively during the years, a systematic study of its performance has been carried out only recently by Gilson’s group^{55} who have found that the performance of QH deteriorates significantly in Cartesian coordinates and when applied to more than one microstate.^{7}

1. Beveridge DL, DiCapua FM. Annu. Rev. Biophys. Biophys. Chem. 1989;18:431. [PubMed]

2. Kollman PA. Chem. Rev. 1993;93:2395.

3. Jorgensen WL. Acc Chem. Res. 1989;22:184.

4. Meirovitch H. In: Reviews in Computational Chemistry. Lipkowitz KB, Boyd DB, editors. Vol. 12. New York: Wiley-VCH; 1998. p. 1.

5. Gilson MK, Given JA, Bush BL, McCammon JA. Biophys. J. 1997;72:1047. [PubMed]

6. Boresch S, Tettinger F, Leitgeb M, Karplus M. J. Phys. Chem. B. 2003;107:9535.

7. Meirovitch H. Curr. Opin. Struct. Biol. 2007;17:181. [PubMed]

8. Gilson MK, Zhou H-X. Ann. Rev. Biophys. Biomol. Struct. 2007;36:21. [PubMed]

9. Foloppe N, Hubbard R. Curr. Med. Chemistry. 2007;13:3583. [PubMed]

10. van Gunsteren WF, Bakowies D, Baron R, Chandrasekhar I, Christen M, Daura X, Gee PJ, Geerke DP, Glättli A, Hünenberger PH, Kastenholz MA, Oostenbrink C, Schenk M, Trzesniak D, van der Vegt NFA, Yu HB. Angew. Chem. Int. Ed. 2006;45:4064. [PubMed]

11. Cheluvaraja S, Meirovitch H. J. Chem. Theory Comput. 2008;4:192.

12. Cheluvaraja S, Mihailescu M, Meirovitch H. J. Phys, Chem. B. 2008;112:9512. [PMC free article] [PubMed]

13. Alder BJ, Wainwright TE. J. Chem. Phys. 1959;31:459.

14. McCammon JA, Gelin BR, Karplus M. Nature. 1977;267:585. [PubMed]

15. Elber R, Karplus M. Science. 1987;235:318. [PubMed]

16. Stillinger FH, Weber TA. Scienc. 1984;225:983. [PubMed]

17. Gō N, Scheraga HA. J. Chem. Phys. 1969;51:4751.

18. Gō N, Scheraga HA. Macromolecules. 1976;9:535.

19. Hagler AT, Stern PS, Sharon R, Becker JM, Naider F. J. Am. Chem. Soc. 1979;101:6842.

20. Karplus M, Kushick JN. Macromolecules. 1981;14:325.

21. Schlitter J. Chem. Phys. Lett. 1993;215:617.

22. Andricioaei I, Kaplus M. J. Chem. Phys. 2001;115:6289.

23. Meirovitch H. Chem. Phys. Lett. 1977;45:389.

24. Meirovitch H, Vásquez M, Scheraga HA. Biopolymers. 1987;26:651. [PubMed]

25. Meirovitch H, Koerber SC, Rivier J, Hagler AT. Biopolymers. 1994;34:815. [PubMed]

26. Meirovitch H. Phys. Rev. A. 1985;32:3709. [PubMed]

27. Meirovitch H, Scheraga HA. J. Chem. Phys. 1986;84:6369.

28. Meirovitch H. J. Chem. Phys. 2001;114:3859.

29. Hnizdo V, Darian E, Fedorowicz A, Demchuk E, Li S, Singh H. J. Comput. Chem. 2007;28:655. [PubMed]

30. Killian BJ, Kravitz JY, Gilson MK. J. Chem. Phys. 2007;127:024107. [PMC free article] [PubMed]

31. White RP, Meirovitch H. J. Chem. Phys. 2004;121:10889. [PubMed]

32. White RP, Meirovitch H. J. Chem. Phys. 2006;124:204108. [PMC free article] [PubMed]

33. White RP, Meirovitch H. J. Chem. Phys. 2005;123:214908. [PMC free article] [PubMed]

34. Cheluvaraja S, Meirovitch H. J. Chem. Phys. 2005;122:054903. [PubMed]

35. Cheluvaraja S, Meirovitch H. J. Phys. Chem. B. 2005;109:21963. [PMC free article] [PubMed]

36. Cheluvaraja S, Meirovitch H. J. Chem. Phys. 2006;125:024905. [PubMed]

37. Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz KM, Jr, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA. J. Am. Chem. Soc. 1995;117:5179.

38. Qiu D, Shenkin PS, Hollinger FP, Still WC. J. Phys. Chem. 1997;101:3005.

39. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J. Chem. Phys. 1983;79:926.

40. Steinbach PJ, Brooks BR. Proc. Natl. Acad. Sci. U.S.A. 1993;90:9135. [PubMed]

41. White RP, Meirovitch H. J. Chem. Theory Comput. 2006;2:1135. [PMC free article] [PubMed]

42. Rosenberry TL. Adv. Enzymol. Relat. Areas Mol. Biol. 1975;43:103. [PubMed]

43. Millard CB, Kryger G, Ordentlich A, Greenblatt HM, Harel M, Raves ML, Segall Y, Barak D, Shafferman A, Silman I, Sussman JL. Biochem. 1999;38:7032. [PubMed]

44. Millard CB, Koellner G, Ordentlich A, Shafferman A, Silman I, Sussman JL. J. Am. Chem. Soc. 1999;121:9883.

45. Chiu YC, Main AR, Dauterman WC. Biochem. Pharmacol. 1969;18:2171. [PubMed]

46. Hart GJ, O’Brien RD. Biochem. 1973;12:2940. [PubMed]

47. Forsberg A, Puu G. Eur. J. Biochem. 1984;140:153. [PubMed]

48. Ordentlich A, Barak D, Kronman C, Flashner Y, Leitner M, Segall Y, Ariel N, Cohen S, Velan B, Shafferman A. J. Biol. Chem. 1996;268:17083. [PubMed]

49. Ordentlich A, Kronman C, Barak D, Stein D, Ariel N, Marcus D, Velan B, Shafferman A. FEBS Lett. 1993;334:215. [PubMed]

50. Raves ML, Harel M, Pang Y-P, Silman I, Kozikowski AP, Sussman JL. Nat. Struct. Biol. 1997;4:57. [PubMed]

51. Carlacci L, Millard CB, Olson MA. Biophys. Chem. 2004;111:143. [PubMed]

52. Olson MA. Proteins. 2004;57:645. [PubMed]

53. Ponder JW. TINKER-software tools for molecular design, version 3.9. St. Louis, Mo: Department of Biochemistry and Molecular Biophysics, Washington University School of Medicine; 2001.

54. Meirovitch H, Alexandrowicz Z. J. Stat. Phys. 1976;15:123.

55. Meirovitch H. J. Chem. Phys. 1999;111:7215.

56. Meirovitch H, Vásquez M, Scheraga HA. Biopolymers. 1988;27:1189. [PubMed]

57. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Oxford: Clarenden Press; 1987.

58. Chang CE, Chen W, Gilson MK. J. Chem. Theory. Comput. 2005;1:1017. [PubMed]

59. Zacharias M, Straatsma TP, McCammon JA. J. Chem. Phys. 1994;100:9025.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |