INTRODUCTION

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

Structure prediction of the transmembrane portion of G-protein coupled receptors (GPCRs) has become accessible to computational techniques with the development of an -carbon template of the helices (Baldwin, 1997) based on the electron cryomicroscopy structure of rhodopsin (Unger et al., 1997). Construction of the transmembrane domain of GPCRs is today a common procedure that follows a reasonably reliable representation of the three-dimensional (3D) disposition of the helices, guided by a set of biophysical rules (Ballesteros and Weinstein, 1995). The topology of GPCR indicates that the helices are connected by three loops on the extra and intracellular sides of the membrane, respectively. However, no structure of the loops in GPCRs is available and due to the variability of their sequences, predictions of loop structure based on computational techniques remains largely unexplored. Nevertheless, several attempts have been made in constructing the loops with considerable success in predicting specific binding sites on the receptor surfaces (Colson et al., 1998; Perlman et al., 1997; Pogozheva et al., 1998; Moro et al., 1999). Clearly, the loops determine ligand selectivity on the extracellular side and initiate signal transduction through interaction with the G-proteins on the intracellular side. A reliable method of predicting their structure and dynamic properties is therefore highly desirable. In the absence of a relationship between sequence and structure, a method based on an energetic approach could be extremely useful. The main obstacles to such a method are the same as those in computational prediction of protein folding, namely the limited sampling in simulation approaches and the representation of the solvent. The sampling problem can be addressed to some extent through the locally enhanced sampling (LES) approximation (Elber and Karplus, 1990) but an effective representation of the solvent at a microscopic level requires further elaboration.

In this paper we address the question of locating the lowest energy minimum of a solvated peptide in an effective way that can also be applied to the determination of the 3D structure of specific parts of large flexible biopolymers. For example, in the case of loops on GPCRs, extensive sampling is required and the influence of the solvent should not be overlooked. To treat these systems we have developed a computational model that combines the LES (Elber and Karplus, 1990) and the SA methods with the primary hydration shell (PHS) approach (Beglov and Roux, 1995).

The LES approximation (Elber and Karplus, 1990; Roitberg and Elber, 1991) is an optimization based on mean-field theory (Gerber et al., 1982) that provides an efficient sampling of parts of peptides and proteins. It can also be coupled with a realistic representation of the aqueous environment. The method has been applied to various problems, such as the determination of diffusion pathways of a ligand through proteins (Elber and Karplus, 1990; Czerminski and Elber, 1991), studies of ligand recombination dynamics (Gibson et al., 1992), free energy calculations (Verkhivker et al., 1992), and modeling of side chains in peptides and proteins (Roitberg and Elber, 1991). In the applications of LES to structure prediction, it has been combined with the simulated annealing (SA) method (Otten and van Ginneken, 1989; Wilson et al., 1991; Morales et al., 1991) to optimize the energy of the systems and in the search for low energy minima of peptides and proteins. The combined LES/SA approach was also applied to structure prediction of solvated peptides with either periodic boundary conditions (PBC) (Simmerling and Elber, 1994, 1995) or with the spherical solvent boundary potential (SSBP) approach (Mohanty et al., 1997).

The recognition that the solvent plays an integral part in structure prediction led to several approaches to include its effect. In some methods, dominant effects of solvation are taken into account only implicitly, whereas in other methods a large number of solvent molecules are included explicitly. When the solvent is treated implicitly, as, for example, in continuum electrostatic approximations (Sharp, 1991; Sklenar et al., 1990; Zauhar and Morgan, 1985), the direct involvement of solvent molecules in influencing the structure of the solute is treated incorrectly (Tirado-Rives and Jorgensen, 1991). An explicit inclusion of solvent molecules is more realistic, but computationally expensive. The two most common approaches in explicit solvent simulations are the PBC (Brooks et al., 1988) and the SSBP (Beglov and Roux, 1994). In the PBC approach, a cell that contains the biomolecule and solvent molecules is surrounded by images of itself imitating an infinite periodic system. In the SSBP approach, the biomolecule is surrounded by a sphere of explicit water molecules, which are restrained by external forces that approximate the mean potential at the boundary of the sphere. The complete energy expression of the system also includes a contribution of the solvent outside the spherical boundary through a continuum approximation. The main problem with including a large number of solvent molecules in LES/SA simulations is the freezing of the torsional degrees of freedom around 500 K due to the fast transfer of kinetic energy from solute to solvent coupled with the slow distribution of kinetic energy throughout the solvent.

The major advantage of the SSBP method is its ability to reduce the number of water molecules required to solvate the solute. However, because the solvent in the SSBP approach is constrained to a fixed spherical shape, it cannot adjust dynamically according to conformational fluctuations of a flexible biopolymer (Beglov and Roux, 1995). In particular, when only a limited solvation is constructed to reduce the number of solvent molecules, the spherical shell can be too restrictive when the biopolymer may undergo significant conformational changes in SA.

Similarly to the SSBP method, the PHS (Beglov and Roux, 1995) method is a compromise between the two extreme treatments of solvation effects. It has, however, the advantage that it includes a flexible nonspherical solvent restraining potential, which eliminates the restrictions of the fixed spherical shape. The PHS, in distinction from the SSBP, has to sacrifice the bulk representation of the solvent and replace it with a restraining force. The restraining force in the PHS model balances the instantaneous pressure inside the primary solvent shell. Because this pressure is temperature-dependent, a generalization of the constant temperature PHS model to a temperature-dependent model is necessary to combine it with SA. Furthermore, combining this method with LES requires careful treatment of the restraining force in view of the existence of multiple copies of various parts of the solute, because the energy expression changes as a function of the number of copies in LES.

Further elucidation of these points is included in the next section, in which the basic ideas of the LES, SA, and PHS methods are summarized, and in the subsequent section, in which analysis of the issues that arise from combining the three methods together is presented. To demonstrate how the combined method, PHS/LES/SA, which we abbreviate as PLS, can be applied we have determined the structure of thyrotropin-releasing hormone (TRH) in water applying the PLS approach. The computational details are included in Computational Protocol, and the results obtained for TRH are presented in Results and Discussion. Finally, the results are summarized in the last section.

	THE LES, SA, AND PHS METHODS: THE BASIC IDEAS

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

The locally enhanced sampling approximation

As noted in the Introduction, the LES method (Elber and Karplus, 1990; Roitberg and Elber, 1991) is based on a mean-field theory (Gerber et al., 1982), in which oversampling is obtained by "multiplying" parts of the peptide or the protein whose structure one wishes to determine. The copies of each multiplied segment do not interact with each other, and each copy of a particular segment is subjected to a force due to its interaction with a copy of another segment that is weighted by the inverse product of the number of copies of both segments. Although the oversampling leads to a new potential surface, the global energy minimum for the new system is identical to the original one, implying that application of the LES approximation may lead to the correct global energy minimum. As a result of oversampling, the barriers separating local minima in the LES calculation are lower, which is a significant advantage when searching for the global minimum of a system. Lowering the barriers reduces the probability that the system would be trapped in a specific local minimum. Moreover, the computational savings that originate from the parallel computations on several system copies are significant. The number of copies chosen is made as a compromise between sampling efficiency and extra computational cost due to the increase in the number of copies. In view of experience described in LES papers, 3-5 copies seem to be an appropriate choice in most cases (e.g., Simmerling and Elber, 1994; Mohanty et al., 1997) although in some cases 10 copies have been applied (Roitberg and Elber, 1994).

Simulated annealing

In general, the potential surface of a molecule with many vibrational degrees of freedom has multiple minima. Indeed, most energy minimization procedures to locate local minima strongly depend on the starting point of the calculation. A grid search (Lipton and Still, 1988), which is a method that systematically scans torsional coordinates, can locate all the local minima, including the global minimum of the system, but is extremely expensive computationally for systems with more than 10 torsional angles.

The SA procedure (Kirkpatrick et al., 1983) is designed to locate the global minimum of the system efficiently, overcoming the problems encountered in other methods and making it applicable to conformational searches and protein folding. It is based on an analogy between the computational optimization of a multivariable function and the experimental annealing procedure for optimization of a crystal, in which an imperfect crystal is heated until it melts and then slowly cooled to form a perfect crystal. Thus, similarly to the experimental annealing method, in SA molecular simulations (Otten and van Ginneken, 1997; Wilson et al., 1991; Morales et al., 1991), the temperature is increased to allow an extended search of the conformation space, and then decreased slowly to anneal the system to its global minimum applying either the Metropolis algorithm (Metropolis et al., 1953) as the sampling method in a Monte Carlo simulation, or the continuity of the equations of motion in a molecular dynamics simulation of a canonical ensemble.

The primary hydration shell potential

The PHS method (Beglov and Roux, 1995) has been developed to provide an economical treatment of solvent effects on arbitrarily shaped peptides. In this method, only a limited number of solvent molecules in the primary hydration shell is included and maintained close to the peptide using a flexible nonspherical half-harmonic restraining potential. The potential has been designed to act between each solvent molecule that escapes the shell and the solute atom with the nearest van der Waals surface in accordance with the following equation:

(1)

In this equation, K_shell is the force constant of the restraining half-harmonic potential, r_ij is the distance between the ith solvent molecule and the jth solute atom, R_vdw^(j) is the van der Waals radius of the jth solute atom, and R_shell is the shell radius. The shell radius is adjusted in response to the instantaneous pressure inside the shell in a manner similar to the changes in the periodic box in constant pressure simulations (for more details see below). In the presence of this potential, the solute experiences a pseudo-bulk solvent environment that mimics conditions of constant pressure when the system is maintained at a constant temperature. As explained below, the shell adapts dynamically to conformational changes of the solute.

	THE TEMPERATURE-DEPENDENT PLS APPROACH

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

Combining SA with the PHS model

As noted above, the primary hydration shell potential has been developed and applied to constant temperature simulations. Because it was intended to mimic an environment of a bulk liquid at constant pressure, adjustments in the shell radius were made to maintain an average constant pressure. Technically, this was done by the following first-order dynamical equation:

(2)

where R_shell is the shell radius, C_relax is the solvent relaxation constant, is the instantaneous boundary force averaged over the solvent molecules related to the pressure inside the primary shell, and F_ref is a restraining boundary force.

According to Eq. 2, when is somewhat smaller than F_ref, which implies that the pressure in the shell is too low, its radius would decrease in accordance with the difference between and F_ref, F, and with the solvent relaxation constant, C_relax. Likewise, when is larger than F_ref, which implies that the pressure in the shell is too high, its radius would increase in accordance with F. Thus, eventually, the water may either escape (if F_ref is too small and there is no other constraint) or collapse to an extremely small volume (if F_ref is too large). Only an adequate value of F_ref will maintain a stable system; that is, the shell radius would fluctuate around an average value even though significant changes in the shell radius could still occur as a result of conformational changes in the solute. However, because simulations are finite, the precise determination of F_ref is unnecessary as long as the change in the shell radius is insignificant during the simulation time.

Clearly, the value of F_ref determines the shell radius, which in turn reflects the density of the solvent. The most sensitive parameter in these simulations was found to be the reference boundary force, and therefore its correct determination is crucial to the calculations. It was found that at 300 K a value of 0.9 kcal/molÅ maintains the appropriate density of bulk water, i.e., 0.0334/ų, using a spherical system of pure water that consisted of 25 water molecules (Beglov and Roux, 1995). However, when the temperature is increased, the water density decreases. Nonetheless, the temperature is increased in SA molecular simulations merely to provide additional kinetic energy to allow the biomolecule conformation to change. In the periodic boundary method with constant volume simulation, this problem is overcome by allowing an image water molecule to replace a water molecule that has escaped. In the present method, images do not exist, and therefore the solvent escape has to be prevented by the restraining potential. Because the pressures that accumulate in a restrained system increase with temperature, the restraining force F_ref should be increased with the temperature to prevent escape of the solvent. However, a slight increase in the equilibrium value of the shell radius is expected to account for a pseudo-linear decrease in the solvent density at higher temperatures. Quantitatively, a linear extrapolation of the changes in the water density below the boiling point to a temperature of 1800 K would amount to ~6% increase in the shell radius, which is ~20% change in volume.

To determine the temperature-dependence of the restraining force, we refer to the definition of the temperature in terms of the mean kinetic energy in molecular dynamics simulations (Brooks et al., 1988):

(3)

where (3N  n) is the total number of unconstrained degrees of freedom in the system, m_i is the mass of atom i, v_i(t) is its velocity at time t, and k_B is the Boltzmann constant.

The instantaneous boundary force is averaged over the solvent molecules that have escaped the hydration shell. It depends on the individual distance that each solvent molecule reached beyond the shell radius. Because the distance to which a solvent molecule would reach depends on its velocity, it follows that the dependence of the instantaneous boundary force on temperature would derive from the dependence of the average velocity of the solvent molecules on temperature.

From Eq. 3 it follows that the temperature-dependence of the average velocity is

(4)

The distance to which a solvent molecule reaches is obtained by integrating the velocity over time

(5)

which results in the following relation between the instantaneous shell force and the velocity:

(6)

Accordingly, the general temperature-dependent form of the instantaneous force, which will be further elaborated in Results, is

(7)

where the constants c₁ = (K_shellt)²c and c₂ = K_shellx₀ reflect the properties of the system including the interactions between the solvent and the solute that determine the actual value of the instantaneous boundary force.

Combining LES with the PHS model

In the PHS approach is calculated by averaging over the forces that act only on the solvent molecules that escape the shell. These forces depend on the distance between each solvent molecule and the closest solute atom. Thus, in an LES simulation, in which some or even all atoms are replicated, this distance is always the shortest and only one atom of the replicas is selected. It therefore follows that the required F_ref will always be smaller than in non-LES simulations, especially at lower temperatures when the potential energy due to the interaction with the LES copies becomes dominant. Consequently, the constants in the general temperature-dependent form of the instantaneous force, Eq. 7, are expected to depend on the number of copies, as will be demonstrated in Results. It should be noted in this context that although only one of the LES copies determines the force that acts on a specific water molecule at a particular time step, in long simulations all equivalent atoms that belong to different copies are treated equivalently. This affects the average distance that determines the force exerted on the solvent molecules, but does not affect the structure determination of the solute when appropriate scaling of the restraining force is made.

Another factor that may affect the restraining force when combining LES with the PHS model is the lack of mass scaling. In the LES approximation the energy and forces are scaled by using the BLOCK commands in CHARMM. The mass can also be scaled using the SCALAR commands. However, unless only equilibrium properties are desired, it has been demonstrated that in the original framework of the LES, mapping a mass-scaled replicated system back to the physical one-copy system at thermal equilibrium results in an overestimated temperature by the factor of the number of copies in the simplest case of a uniform weighting of all the copies (CHARMM c26bl user manual). Therefore, mass scaling is not always desirable and we leave the masses unscaled. This, however, is expected to affect the restraining force required in a replicated system as compared with an unreplicated one through c₁. Further discussion of the thermal properties of LES can be found in the papers by Straub and Karplus (1991) and Ulitsky and Elber (1993).

The effect of the hydration shell thickness in the PHS model

The thickness of the hydration shell is an important parameter of the PHS model because the final results of a simulation using this approach should be independent of the thickness of the shell. As the thickness of the shell increases, the outer solvent molecules on which the restraining force is applied become more distant from the solute. Therefore, the non-bonding interactions between these external solvent molecules and the solute decrease, an effect that should be compensated by increasing the restraining force, as will be demonstrated in Results as well. This increase should reach a limiting value when the external solvent molecules experience the rest of the system as a pure solvent. The thickness-dependence of the restraining force should be expressed as the dependence of c₂ on the shell radius in the expression for the instantaneous force, Eq. 7.

	COMPUTATIONAL PROTOCOL

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

The molecular dynamics simulations were performed using CHARMM (Brooks et al., 1983) version 26. Small changes were made in the PHS-related sections of the CHARMM program to obtain a comprehensive temperature- and/or time-dependent analysis of the instantaneous shell force and of the shell radius, as well as to control the restraining force adjustments in temperature-dependent SA simulations.

The tests of the combined PLS approach were performed on the TRH, whose initial structure for optimization was based on its x-ray structure (Kamiya et al., 1980; see Fig. 1 a). The initial structure of TRH was placed at the center of a 20 Å sphere of 1024 TIP3P (Jorgensen et al., 1983) water molecules. Water molecules that were within a distance of <2.6 Å from the TRH (17 waters) were deleted and the combined TRH-water system was then minimized in four stages. Initially, the TRH was kept frozen while its surrounding water was minimized using the steepest descent (SD) algorithm for 4000 steps followed by the adopted-basis Newton Raphson (ABNR) algorithm for an additional 4000 steps. This was followed by a minimization of the entire system by applying the SD and the ABNR procedures consecutively for 4000 steps each after releasing the TRH. An additional 933 water molecules that were at a distance larger than 5 Å from any TRH heavy atom were deleted, leaving the system with 74 water molecules that constitute the primary hydration shell (see Fig. 1 b).

View larger version (41K): [in this window] [in a new window]   FIGURE 1   Thyrotropin releasing hormone (TRH). (a) Initial conformation based on the x-ray structure; (b) same as (a) surrounded by a 5 Å water shell; (c) 5-LES-copy structure b (see text) at 1800 K; (d) final 5-LES-copy structure b obtained after cooling from 1800 K to 300 K.

The TRH was then replicated in three parts according to the residues that it comprises: pyroglutamic acid, histidine, and prolineamide. The boundaries between the parts were at the peptide bond that connects the residues. Each residue was multiplied in most simulations 5 times, but to test the dependence of the restraining force on the number of copies, TRH residues were multiplied 1, 2, 3, 4, and 10 times.

The temperature of the entire system was raised from 0 K to 1800 K using the leap-frog Verlet integrator with Langevin dynamics used in the presence of the miscellaneous mean-field potential (MMFP) in CHARMM, keeping the center of mass of TRH fixed and allowing the application of the PHS model. A time step of 0.001 ps was used in all simulations, the hydrogen atoms were subjected to the SHAKE constraint, and a friction constant of 25 ps was applied to the water oxygen atoms. The same friction constant was used in the original implementation of the PHS model by Beglov and Roux (1995). However, to prevent overdamping the motions of the TRH, we have applied the Langevin force only to the oxygen atoms of the solvent. K_shell was 3 kcal/molŲ and C_relax was 0.5 Å²/(pskcal/mol). The Langevin heat bath temperature was increased every 0.4 ps by 10° along with the restraining force of the shell potential, which was increased gradually from 0.12 to 3 kcal/molÅ.

The system was equilibrated at 1800 K for 450 ps, after which it was gradually cooled to 300 K at a rate of 10 K/3 ps from several different structures. Examination of the energy-dependence on temperature shows that it is monotonic and approximately linear. As the cooling process progressed, the restraining force was adjusted. However, because the SA process is much slower than the heating stage, the adjustments of the restraining force require a higher accuracy. Therefore, the details of the adjustments in the restraining force vary for the different conditions and will be specified for each case individually in the following section.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

Structure of TRH

To establish the applicability of the combined PLS method, we have conducted two independent SA of TRH in the PHS, starting with two different structures extracted from the trajectory at 1800 K. One structure was selected after 258.7 ps (structure a) and the other after 425.8 ps (structure b, see Fig. 1 c). To evaluate the enhanced sampling provided by LES we have calculated the all-atom rms deviations among all the possible conformations of these structures at 1800 K. Because there are 3 residues, each multiplied in 5 copies, the total number of conformations for each structure is 125 (5³). The (a,a) and (b,b) diagonal corners in Fig. 2 a show the rms deviations among the 250 conformations within the a and the b structures, respectively. The (a,b) off-diagonal corner shows the rms deviations among the 125 conformations of the a and b structures. At 1800 K the maximal rms deviations between the conformations of each structure is ~1.6 Å in a and b and reaches ~2.3 Å between the a and b structures. This illustrates the extent of the enhanced sampling at the high temperature. For each of these structures, the system was cooled at a rate of 10 K/3 ps to 300 K, decreasing the restraining force along the process in a manner that will be specified below. The all-atom rms deviation of the structures decreased gradually with the decrease in temperature, until at 300 K the structures nearly converged and could no longer be regarded as different structures, as can be seen from the rms deviation shown in Fig. 2 b. After the SA, the conformations within each structure converge to an rms deviation that does not exceed 0.28 Å and the difference between the structures does not exceed 0.39 Å. An approximate analysis that can provide insight into the trend of the rms deviations with temperature can be obtained by analyzing only 5 structures (instead of 125) obtained from the ith copies of the multiplied segments. For the present case, this would imply PGL(i)-HIS(i)-PRO(i), i = 1, ... , 5. The rms deviations at 1800 K between each of the copies of structures a and b range from 1.49 Å to 2.13 Å and within the structures they reach 1.4 Å. At 300 K these deviations reduce to 0.25 Å between copies of the same structure and 0.37 Å between structures a and b. Also, the overall rms deviation among the 5 copies of each structure decreased from 2.40 Å to 0.61 Å as the temperature was reduced from 1800 K to 300 K. Additional information regarding the convergence of the structures at lower temperature can be obtained from a cluster analysis of the conformations of structures a and b, which is presented in Fig. 2 c for 1800 K and in Fig. 2 d for 300 K. It can be seen that for the clustering value of 0.3937 Å, at which level the SA can be considered converged and which can thus be regarded as only one family at 300 K, there are six families at 1800 K. The resultant structure b at 300 K is drawn in Fig. 1 d. In this context it should be noted that previous LES/SA described by Simmerling and Elber (1994) found correct structures in simulations as short as 250 ps. We have performed simulations at a cooling rate of 10 K/2 ps for a total of 300 ps, but found them to be too short. The SA starting from two different structures (e.g., a and b) did not converge to the same final structure with an rms deviation at 300 K exceeding 1 Å).

View larger version (122K): [in this window] [in a new window]   FIGURE 2   RMS deviations between the 125 conformations of structure a and the 125 conformations of structure b of TRH in a five-copy LES simulation (a) extracted from the equilibration at 1800 K, (b) after cooling to 300 K, (c) same as (a) but in generic conformation order. The number of clustered families at the level of 0.3937 Å is shown on the bar below the matrix. (d) Same as (b) but in generic conformation order. Note the different ranges in the color scales at the bottom of the RMS matrices.

The details of the annealing process were further examined by following the variations in seven torsional angles of the TRH₁, ₁, ₂, ₂, ₃, ₁ and ₂with temperature. The results for 10 copies, which originated from structures a and b, are presented in Fig. 3, a-g. From Fig. 3 it is evident that the angle distributions, which ranged from 50° to 150° at 1800 K, converge at 300 K to the same value for each angle. The convergence of torsional angles can be presented by plotting the circular variance of these angles as a function of temperature, as shown in Fig. 3 h. The circular variance (Mardia, 1972; Allen and Johnson, 1991; MacArthur and Thornton, 1993) is defined as

(8)

and can serve as a continuous measure of the variability in the torsion angle space. This analysis takes proper account of the periodic nature of the angular domain instead of the regular linear statistics that cannot be applied here because the range of the angles is 0°    360°. The maximal variance for n data points is 1, and it approaches 0 when the variance is diminished. Fig. 3 h clearly demonstrates that all the angles converge at the end of the SA.

View larger version (56K): [in this window] [in a new window]   FIGURE 3   Annealing history of TRH represented as distribution of torsional angles of the 5-LES-copies of structures a and b. (a) 1 (b) 1 (c) 2 (d) 2 (e) 1 (f) 2 (g) 3. (h) Circular variance analysis of the torsional angle history presented in (a)-(g). The lines were added to clarify the general trend of the changes in variance with temperature. They do not imply a functional relationship.

To examine the reliability of the PLS method we have performed a 10-ns simulation starting from copy 1 of structure b at 300 K, using the PHS method. Due to the small size of the system and to the minimal number of water molecules included, this simulation was accomplished in approximately two days. The restraining force used was 0.87 kcal/molÅ and an average shell radius of 4.5 ± 0.5 Å was maintained during the simulation. The distribution and relative populations of the seven torsional angles listed above during the PHS simulation are given in Table 1. Graphical presentation of the evolution of three of these angles, ₁, ₂, and ₁, with time is given in Fig. 4, a-f, along with the angle histogram and its integration. The evolution of these angles represents the general characteristics of the system, demonstrating the high flexibility of some angles and the relative rigidity of others. Monitoring the variations in the distribution of the water molecules around the TRH during the simulation, we observed that the shell provides an efficient effective solvation shell around the TRH adjusting dynamically to changes in the conformations of TRH, and rearranging continuously. To determine the most stable conformation of TRH, we have analyzed the possible combinations of the different angular ranges of the most flexible torsional angles, ₁, ₃, ₁, and ₂. The sample space included 10,000 evenly spaced conformations recorded every 1 ps during the 10-ns simulation. The combination that included the highest number of conformations was 270°  ₁ < 330°, 140°  ₃ < 180°, 270°  ₁ < 330°, and 60°  ₂ < 160°. It included ~14% of the conformations. The next two most abundant conformation ranges were 150°  ₁ < 220°, 140°  ₃ < 180°, 160°  ₁ < 210°, 60°  ₂ < 160°, and 150°  ₁ < 220°, 140°  ₃ < 180°, 160°  ₁ < 210°, and 200°  ₂ < 290°. They included ~9% and 7% of the conformations, respectively. Six other conformation ranges were represented by 2%-4% and many others by <2%. This analysis leads to the conclusion that the energy barriers between different conformations is rather low, allowing a wide distribution among a large number of conformations.

View larger version (74K): [in this window] [in a new window]   FIGURE 4   PHS 10-ns simulation of TRH at 300 K, starting from the final 300 K PLS structure b with a restraining force of 0.87 kcal/molÅ as manifested in the trajectory of three representative torsional angles. (a) 1 Fluctuations, (b) 1 histogram and integration, (c) 2 fluctuations, (d) 2 histogram and integration, (e) 1 fluctuations, (f) 1 histogram and integration.

The results obtained in the PHS simulation can be compared to experimental results for TRH (Deslauries et al., 1973; Donzel et al., 1974; Montagut et al., 1974; Ward et al., 1986; Vicar et al., 1979; Unkefer et al., 1983). The experimental results were obtained by ¹H- or ¹³C-NMR spectroscopy and analyzed on the basis of the Karplus relationship between coupling constants and torsional angles (Karplus and Karplus, 1972; Bystrov et al., 1969). These yield several possible values from each measurement that are listed in Table 1. Comparing these results to the PHS results, we note that most of the PHS results are in the vicinity of an experimental result, and that those with the highest relative populations in the PHS simulation all match one of the experimental values within a 5° error in most cases and a 40° error in ₁, which is the most varying angle. Following the criterion for comparing experimental and computational results set previously (Yao et al., 1994a-c; Mohanty et al., 1997), i.e., that if the deviations from the experimentally suggested structure of the backbone torsional angles are smaller than 50°, the structure is assumed to be correct, the results from the PHS simulation agree very well with the experimental results, and can even aid the evaluation of the currently available inconclusive experimental results. In addition, even in the limit of the stricter criterion in which the allowed error is reduced to 40°, the prediction accuracy remains 100%.

Comparing the results obtained from the 5-copy PLS simulation for the torsional angles, Fig. 3, with those obtained from PHS simulation and with the experimental results, it is clear that the PLS simulation leads to one of the PHS possibilities for each of the torsional angles, which also appears among the experimental results. The range of results obtained from the 5-copy PLS simulations that originated from structures a and b is summarized in Table 1. In addition, the torsional angles of the TRH have also been extracted from the 1-4-copy PLS simulations and included in Table 1. The results demonstrate that in all cases the angles fall within the PHS and experimental ranges, although results obtained in different cases may differ. Comparing these PLS results to the analysis regarding the most probable conformation in the PHS simulation, it follows that the prediction accuracy of the most stable conformation in all simulations is above 71% and reaches 86% for the 1-copy PLS simulation. In view of the low energy barriers in the TRH as manifested in the PHS simulation, the deviation from the global minima is within reasonable limits and as noted above, the prediction accuracy for determining a possible conformation is 100% for all cases. This indicates that in the case of TRH, which is a small tripeptide, LES is not an essential ingredient, and that satisfactory results can be obtained without it. Nonetheless, according to the present analysis, applying LES yielded correct results, which is an important observation with respect to further application of the PLS method to more complicated molecular systems in which the application of LES is extremely useful (see, for example, Roitberg and Elber, 1994).

As discussed above, the results should be independent of the shell thickness. Thus, the combined 5-copy PLS method was applied to shells of varying thickness and the final results were compared to those obtained with a 5 Å shell. The 6 Å water shell consisted of 121 water molecules and the system was prepared in a similar way to the 5 Å shell. A snapshot at 16.7 ps was selected from the equilibration at 1800 K. The initial overall rms deviations at 1800 K between the 6 Å shell structure and structures a and b were 1.93 Å and 1.62 Å, respectively. These deviations decreased as the temperature was lowered and reached the respective values of 0.59 Å and 0.40 Å at 300 K. Further application of the PLS method with 8 Å (239 water molecules) and 10 Å (410 water molecules) water shells led to final structures with overall rms deviations from structure a of 0.64 Å and 0.54 Å, respectively, and of 0.34 Å and 0.45 Å from structure b (see Table 3). The torsional angles of the structures obtained with the 6 Å, 8 Å, and 10 Å shells all fall in the ranges obtained from the 5 Å 5-copy PLS simulation. These results lead us to conclude that the final structure does not depend on the shell thickness.

Properties of F_ref

The calibration of F_ref as a function of temperature is necessary to perform an SA simulation in the combined PLS method. Other than its general dependence on , little is known about c₁ and c₂. We have calibrated these constants with the requirement that R_shell will remain nearly constant throughout the SA from 1800 K to 300 K. An example of a calibrated change of the restraining force during the cooling process is plotted as a function of temperature in Fig. 5. This was accomplished for a 5 Å shell and 5-LES copies of TRH with c₁ = 0.0118 (kcal/molÅ)²K¹ and c₂ = 1.54 kcal/molÅ. To illustrate the stability of R_shell, the trajectory of the shell radius obtained in the simulation of structure b is included in Fig. 5. The trajectory demonstrates that this choice of varying the restraining force with temperature results in a stable shell with an average radius of 5.2 ± 0.5 Å.

View larger version (30K): [in this window] [in a new window]   FIGURE 5   Calibrated restraining force for the simulation of the 5-LES-copy structure b as a function of temperature. The resulting shell radius trajectory is also displayed.

To test the dependence of the restraining force on the number of copies in the LES, the same restraining force as above was applied to SA of TRH with a different number of copies. The variation of the shell radius in these calculations is shown in Fig. 6 a. It emphasizes the necessity of a correct calibration of the restraining force with temperature. In all cases where the number of copies was less than five, the shell radius has increased to large values during the annealing because the restraining force was too weak. In contrast, in a simulation with 10 copies the shell contracted because the force was too strong. Note that the instantaneous shell force, presented in Fig. 6 b, is almost the same for all the cases, illustrating that the shell expands or contracts to minimize the difference between the instantaneous shell force and the restraining force. In constant temperature simulation this quality is highly desirable because it enables the shell to respond to conformational changes in the solute while maintaining a constant average shell radius. However, in SA, where the changes in temperature produce large changes in the kinetic energy of the solvent and the shell radius, the restraining force has to be adjusted to maintain a nearly constant shell radius. Moreover, because in LES the attractive forces between the solute and the solvent are scaled inversely to the number of copies, the forces optimized for a given number of copies are out of balance when their number is changed. To investigate the dependence of the restraining force on the number of copies, it was calibrated separately for each of the simulations. The instantaneous shell force that reflects the restraining force is depicted for each of the cases in Fig. 7 a. These curves were fitted according to Eq. 7, and the results are summarized in Table 2. It is clear that c₁ has a weak dependence on the number of copies. In contrast, the main dependence on the number of copies is reflected in c₂. The corresponding shell radius for each of these cases is presented as a function of the temperature in Fig. 7 b and their average values are included in Table 2. The fact that the values of the shell radius do not deviate significantly from 5 Å indicates that the calibration was adequate.

View larger version (46K): [in this window] [in a new window]   FIGURE 6   The results of the application of the restraining force plotted in Fig. 4 on 1,2,3,4,10-LES-copy structures. (a) The shell radius trajectory; (b) the instantaneous shell force.

View larger version (53K): [in this window] [in a new window]   FIGURE 7   Individually calibrated multi-LES-copy simulations. (a) Instantaneous shell forces and curves fitted to Eq. 7; (b) shell radia trajectories; (c) dependence of c1 and c2 on the number of copies; (d) sensitivity of c1 and c2 on the number of copies expressed as percentage of increase compared to the values for one copy.

The dependence of c₁ and c₂ on the number of copies, n_rep, appears to be linear, as shown in Fig. 7 c. The equations of the fitted lines are c₁ = 0.0064 + 0.0010 n_rep and c₂ = 0.417 + 0.196 n_rep, both with a correlation coefficient R = 0.987. As noted above, the constants c₁ and c₂ do not show a similar sensitivity to n_rep, with c₂ being more sensitive to the number of copies, as can also be seen in Fig. 7 d. The different dependence originates from the distinct physical nature of these constants. The coefficient c₁ reflects the dependence of the kinetic energy on temperature. It is therefore related directly to the total number of unconstrained degrees of freedom and inversely to the total mass of the system. Increasing the number of copies in the LES method increases the degrees of freedom of the solute molecule and the mass. However, because increasing the number of copies affects only the degrees of freedom of the solute, the two effects nearly cancel each other resulting in an overall weak dependence on n_rep. The coefficient c₂ is an integration constant of the velocity expression and can be considered as the force required to maintain the system at a specified radius in the absence of kinetic energy. Not surprisingly, this force is negative because in the absence of kinetic energy the attraction between the molecules in the system will contract its dimensions. Practically, however, because the PLS computational procedure is based on a half-harmonic potential that acts only when water molecules reach beyond the outer radius of the shell, it is incapable of adjusting the water shell appropriately below 300 K and should therefore not be applied at this temperature range. It is interesting to compare the results from our calibration with those of Beglov and Roux (1995). In their system of 25 water molecules they used a restraining force of 0.9 kcal/molÅ. Because both the TRH and water are neutral and polar molecules, it is possible to obtain the restraining force at 300 K from our calibration procedure as the limit of n_rep = 0, which results in 0.97 kcal/molÅ.

To test the dependence of the restraining force on the thickness of the shell, systems with 6 Å, 8 Å, and 10 Å shells were used. TRH was replicated in five copies as before, and the system was heated and equilibrated in a way similar to the simulations with a 5 Å shell. An initial structure for cooling was selected after 16.7 ps of equilibration at 1800 K and cooled to 300 K with a similar protocol. The trajectories of the shell radii were stable during the SA, with averages of 6.3 ± 0.5 Å for the 6 Å shell, 8.3 ± 0.7 Å for the 8 Å shell, and 10.6 ± 0.3 Å for the 10 Å shell. The trajectories of the shell radius for the various shell thicknesses, including the result for the 5 Å shell (structure b) are plotted in Fig. 8 a, and the corresponding trajectories of the instantaneous shell forces are depicted in Fig. 8 b. Because c₁ should have a very weak dependence on shell thickness, all four curves have been fitted to the equation  =   c₂, where c₁ was selected from the 5 Å shell simulations. The fittings are excellent for all cases and the results obtained for the c₂ coefficients are summarized in Table 3. Plotting c₂ as a function of the average shell radius (see Fig. 8 c) shows that for a 5-copy PLS simulation, it can be fitted to the equation c₂ = 6.049 exp(0.63) + 1.308 with a correlation coefficient of 0.9999, and that c₂ converges to 1.308 kcal/molÅ as the shell radius approaches infinity. As noted above, the stronger restraining force required for the thicker shells originates in the weaker interaction forces that act between the TRH atoms and the outer water molecules due to a larger distance between them. The limiting value originates from the fact that as the shell thickness increases, the outer water molecules, whose motion is manifested in the changes in the shell radius, cease to interact with the solute and experience their surroundings as pure solvent.

View larger version (34K): [in this window] [in a new window]   FIGURE 8   Individually calibrated 5-LES-copy simulations with different shell thickness. (a) Shell radia trajectories; (b) instantaneous shell forces and curves fitted to Eq. 7 with c1 = 0.01183 (kcal/molÅ)2K1; (c) dependence of c2 on the number of the average shell radius.

	SUMMARY

TOP ABSTRACT INTRODUCTION THE LES, SA, AND... THE TEMPERATURE-DEPENDENT PLS... COMPUTATIONAL PROTOCOL RESULTS AND DISCUSSION SUMMARY REFERENCES

The primary hydration approach has been combined with locally enhanced sampling and with simulated annealing into a unified method that is expected to be effective in treating folding of specific solvated segments of large molecular systems. This was achieved by generalizing the PHS approach to a temperature-dependent form and calibrating the restraining force as a function of temperature. The relationship between force and temperature was derived from the statistical mechanical definition of average temperature. The model was further adjusted to take into account the number of LES copies by demonstrating that the constants in the expression for the instantaneous boundary force are linearly dependent on the number of copies. The dependence on the thickness of the shell of the restraining force was also explored, confirming that the restraining force should be increased with shell thickness. The method has been applied to predict the structure of the tripeptide TRH in solution with very good agreement to known experimental results and with those derived from a 10-ns simulation applying the PHS method at 300 K. These results suggest that the PLS method can be applied to structure prediction of macromolecules. However, further extension of the PLS method is required, because for macromolecules in which folding can significantly change the accessible surface area, maintaining a constant average shell radius during the SA may not be a desirable feature. Present efforts include the development of an appropriate description of the shell that will readjust to large changes in the surface area occurring during SA.