INTRODUCTION

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Structural, dynamic, and folding properties of multielement objects such as proteins are often conveniently referred to the complex landscape of their appropriate configurational energy. In the solvent, the appropriate landscape is that of the potential of mean force (PMF), that is of the configurational potential energy of the whole solute + solvent system, thermodynamically averaged over all solvent configurations (Dill et al., 1995; Frauenfelder et al., 1991; Bryngelson et al., 1995). As a consequence of the size of solvent-induced interactions, the landscapes of potential energy and of free energy can be expected to differ, even substantially. The difference is due, of course, to extra terms of enthalpy and entropy contributed by those solvent molecules that interact sizably with solute elements. These molecules (hydration water) act from a nonuniform distribution in space, due to the statistically favored configurations determined by constraints imposed by solutes. In the course of folding, as well as in conformational switching (frequently associated to function), the hydration pattern (and related free energy) and the protein conformation will be closely interdependent. It has been shown in previous work that hydration and related solvent-induced interactions and forces possess a strong non-pair-additive, manybody character (Brugé et al., 1996b; Martorana et al., 1996; San Biagio et al., 1998). Non-pair-additivity is responsible for very unexpected features of hydration and solvent-induced interactions, such as strong context dependence and long-range propagation (Martorana et al., 1996, 1997, 1998; Bulone et al., 1997; San Biagio et al., 1998). On the other hand, non-pair- additivity is also expected on the basis of the correlated energy landscape model (Plotikin et al., 1996, 1997; Shoemaker and Wolynes, 1999; Shoemaker et al., 1999). One of the purposes of the present work was to explore unambiguously the role of untruncated solvent-induced interactions in generating (rather than modifying) biologically interesting features in the PMF landscape of a simple solute, starting from a flat landscape of potential energy. A second, not minor purpose was to study the dynamical coupling between solute and solvent configurations, which is a question of central interest that has scarcely been explored.

For a better understanding of higher-order terms in solvent-induced interactions, responsible for the remarkable nonadditivity of hydration and related PMF, one can use Stillinger's expression of the free energy of water. Let us refer to the configurational potential energy landscape of a system of N unperturbed water molecules. This landscape contains a multitude of minima and surrounding basins. At equilibrium, only a small (still very numerous) subset of them, having essentially the same depth, is occupied with overwhelming probability. The given thermodynamic conditions determine this subset. As a result of thermodynamic averaging, the total free energy can be expressed (Weber and Stillinger, 1984; Stillinger 1988) in terms of depth, _m, and logarithmic multiplicity, (_m), of these basins and of their related vibrational contributions, f(T, _m), that is,

(1)

Interaction with solutes alters the potential energy landscape and its populated basins, causing a free energy change G_sw that is the free energy of hydration. In the case of n solutes or of a composite solute made of N_p elements (such as a protein) in a fixed configuration R₁, R₂, ... , R_Np, this hydration free energy contains individual, pairwise, and manybody terms expressing solvent-induced solute-solute interactions. It can be written as

(2)

Pairwise as well as higher-order terms in Eq. 2 are strongly dependent upon the R_i variables, which describe the specific configurations of solutes or of solute elements. The recently demonstrated strong nonadditivity of solvent-induced interactions (Brugé et al., 1996a,b; Bulone et al., 1997; Martorana et al., 1996, 1997, 1998; San Biagio et al., 1998) corresponds to unexpectedly large sizes of manybody terms (third and higher order) in Eq. 2. Because of their size, these terms can affect in a substantial way the configurational PMF landscape of a solute having internal degrees of freedom (e.g., a protein) and generate, for example, folding pathways that would not be practicable in their absence. This emphasizes the need for not using, whenever possible, approximations neglecting high-order terms, such as the widely used Kirkwood's superposition approximation (Hill, 1956) and related approaches. More sophisticated and efficient methods based on expansions in terms of pair and triplet correlation functions and proximity approximations (Garde et al., 1996; Garcia et al., 1997) prove adequate for relatively coarse-grained studies. As will be taken up again at the end of the Discussion, however, they may miss all-important details on the microscopic scale. It must be remarked that studies concerning, for example, relatively simple but realistically modeled solutes, such as di- or tripeptides, are based on approximate methods (Pettitt and Karplus, 1988; Perkyns and Pettitt, 1995; Pellegrini and Doniach, 1995; for a comparative discussion of different approaches, see Smith and Pettitt, 1994). Such studies have proved valuable in evidencing, already within the given approximation, significant differences between energy and PMF landscapes. Nevertheless, they do not account for the strong manybody, nonadditive, and long-range character of hydration and related solvent-induced forces (SIFs) evidenced more recently, as quoted above. Third- and higher order terms responsible for the nonadditivity of solvent-induced interactions are implicitly included in recent molecular dynamics (MD) studies of the hydration free energy of rigid hydrophobic model solutes in explicit solvent and its dependence upon solute size, shape, and charge (Wallqvist and Berne, 1995a,b; Wallqvist and Covell, 1995; Lynden-Bell and Rosaiah, 1997). The same is true for further studies of oligopeptides, focusing on specific configurations or trajectories (Tobias and Brooks, 1992; Duan and Kollman, 1998). These studies do not comprise, however, the entire free energy landscape, and they do not elicit the dynamical coupling between solute configurational changes and solvent reorganization.

In the present MD work, as in our previous research, we use an explicit molecular solvent and take into account solvent-induced interactions to all orders in Eq. 2. The solute used here has a very simplified multielement structure (including charged and apolar elements) and two internal degrees of freedom. In vacuo, the corresponding landscape of potential energy is flat. This allows us to elicit unambiguously at an elementary level the role of solvent in generating features that are very significant in the case of proteins and the dynamics of solute-solvent coupling. Specifically, we investigate 1) the full role of solvent in transforming a flat landscape of potential energy into a nontrivial landscape of configurational PMF involving PMF differences not far from those stabilizing the functional conformation of proteins (Note that in solutes of the simple type used here, as well as in more realistic ones, solvent-induced interactions propagate over the entire solute, as a consequence of their manybody character (Martorana et al., 1996, 1997, 1998).); 2) the contribution of the solute conformational reaction to reaching the state of minimal free energy, and 3) the structural and dynamic interplay between solute conformational changes and solvent reorganization. The present work also illustrates solvent-induced interactions between charged and apolar solute elements and effects related to the dependence of such interactions upon the sign of charges, a feature not present in continuum solvent modeling and not fully elicited in truncated PMF calculations. A final MD run using a different modeling of the explicit solvent (TIP4P instead of TIP3P) proves that these features are not critically dependent upon the precise modeling of the explicit water molecule.

	SIMULATION DETAILS

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We use two slightly different modifications, A and B, of one composite model solute in a bath of TIP3P water (Jorgensen et al., 1983). The basic model solute consists of six identical and fixed Lennard-Jones (LJ) hydrophobic spheres, described by the same LJ parameters of the solvent water and lying in the planar arrangement shown in Fig. 1, where the center-to-center nearest-neighbor distance is 4.6 Å (Martorana et al., 1997). Element 4 bears an electric dipole, modeled by a negative charge at its center and an off-center positive charge. The latter is free to move on a spherical surface, the radius of which (0.95 Å) is smaller than the LJ radius. In the case of solute A, the positive and negative charges are equal (0.47 a.u.), so that the total charge is zero. In the case of solute B, the negative charge value is twice that of the positive charge (1.04 and 0.52 respectively). Accordingly, solutes A and B can be taken as modeling a composite hydrophobic solute carrying an OH group or an OH group, respectively. Comparison of results relative to the two cases shows the effect of the additional negative charge, ceteris paribus. An additional run was performed to test the role of the specific potential used for modeling the water molecules. In this additional run, we used the same composite solute B and the TIP4P water-water potential (Jorgensen et al., 1983). Because results do not change the qualitative conclusions reached with TIP3P, they are not reported here in detail.

View larger version (29K): [in this window] [in a new window]   FIGURE 1   Configuration of the six LJ spheres of the planar basic model solute. The center-to-center nearest-neighbor distance is 4.6 Å. Sphere 4 bears a negative charge, fixed in the center, and a positive charge free to move on the surface of an inner sphere, as shown. The negative charge value is equal to that of the positive charge in the A case, and it is twice the positive charge in the B case. The z axis is orthogonal to the figure plane. The spherical coordinates  and  of the positive charge point, i.e., of the dipole orientation, are defined as usual and are referred to this orthogonal system.

Simulations were performed using the Amber package, version 5.0 (Case et al., 1997). The thermodynamic ensemble was N, V, T, at 298 K. The simulation box was 33.9 × 28.4 ×24.9 ų and contained 797 water molecules, so that no less than three water layers surrounded the solute. We used a 12-Å cutoff of the interaction potentials, and periodic boundary conditions. Starting from a configuration (obtained by replicating an equilibrated box of 216 molecules), we performed a further 20-ps equilibration at constant pressure and temperature. After equilibration, each trajectory was at least 2.6 ns long. Time steps were 1 fs, and one configuration every 20 fs was stored for data analysis.

In our presentation the angular orientation of the dipoles is given in terms of the  and  angles defined as usual, with respect to the orthogonal axes shown in Fig. 1. The related angular distribution function, g(, ), was computed from the number n(, ) of times in which the dipole orientation fell within a 10° × 10° box, around a particular (, ) point on the configurational surface. This value was normalized as

(3)

(where N is the total number of analyzed configurations). The related free energy term was computed as

(4)

This is the - and -dependent part of the G_SW quantity expressed by Eq. 2, and it contains contributions to all orders coming from all elements of our solute. Hydration was computed as the distribution of space occupancy probability, p, of water oxygen and hydrogen atoms, normalized to that of bulk water.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Our basic model solute is schematically shown in Fig. 1. In vacuo, the dipole orientation distribution is homogeneous and the potential energy surface is flat, because electric dipoles and charges do not interact with apolar LJ particles. In the solvent, two well-developed regions of maximum probability appear, as shown in the three-dimensional view of Fig. 2. The reflection symmetry of the solute is reproduced in the g(, ) raw data within an overall "noise" that decreases as expected for increasing trajectory lengths. In our case this noise is 30% or less, and it is in large part removed from data in Figs. 2-4 by imposing reflection symmetry with respect to the solute plane. A first comparison of data relative to cases A and B is possible from Figs. 2 and 3. We see that the effect of the additional negative charge is a somewhat higher localization (with no measurable shift) of the statistically favored orientations of the dipole.

View larger version (28K): [in this window] [in a new window]   FIGURE 2   Contour plot representations of g(, ) evidencing the two positions of maximum g(, ). (Left) Solute A (uncharged). (Right) Solute B (negatively charged).

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Normalized distribution function of the -coordinate of the dipole. The values of n()/N indicate the number of configurations corresponding to a -value in a given 10° interval, divided by the total number of analyzed MD configurations. Thin line, solute A; thick line, solute B.

The free energy contribution related to data of Fig. 2 is shown in Fig. 4 in its dependence upon  and . It is due to the interaction of the dipole with explicit water, as anticipated in the Introduction. The landscape exhibits two minima, symmetrically located with respect to the plane of the solute. For symmetry reasons, these minima of course would not exist if the dipole-bearing element 4 were alone in the solvent. Therefore, they are due to the interaction of the dipole with the solvent perturbed by all other (apolar) solute elements. Equivalently, they can be viewed as being caused by solvent-induced interactions between the dipole and the remaining LJ spheres, as described by Eq. 2.

View larger version (74K): [in this window] [in a new window]   FIGURE 4   Contour plots (left) and three-dimensional views (right) of the free-energy landscape in the configurational , plane. Top, Solute A (uncharged) in TIP3P water; center, solute B (negatively charged) in TIP3P water; bottom, solute B (negatively charged) in TIP4P water.

The landscape in Fig. 4 is such that, starting from any dipole orientation, the system is preferentially driven in either of the two minima. A well-defined connecting path between minima is also visible in the figure. Comparison of landscapes relative to cases A and B with that relative to solute B in TIP4P water (also shown in Fig. 4) leads to the conclusion that these essential features are qualitatively independent of the specific modeling of water molecules. In quantitative terms, the presence of the excess negative charge enhances the depth of the minima (3.6 kJ mol¹ in the case of solute B, compared to 3 kJ mol¹ in the case of solute A). The related "free energy of activation" G relative to configurational switching from one minimum to the other is ~2.5 kJ mol¹ in case A and 2.8 kJ mol¹ in case B. The corresponding values for solute B in TIP4P water are 4.5 kJ mol¹ for the depth of the minima and 4.5 kJ mol¹ for the switching. Notably, and notwithstanding the elementary structure of our solute, these values are not far from those of the free energy stabilizing proteins' functional conformations.

As visible from Fig. 2, the favored dipole orientations correspond to the positive charge pointing more toward other solute elements than toward the solvent. This allows a larger exposure of the negative charge to solvent. This is related to the stronger interaction of negative charges with water (Bulone et al., 1997), in agreement with computational and experimental data (Migliore et al., 1988; Straatsma and Berendsen, 1988; Friedman and Krishnan, 1973; Lynden-Bell and Rosaiah, 1997; Marcus, 1994) on the hydration free energy of positive and negative ions, and it is caused by the asymmetrical charge distribution on water molecules.

Hydration patterns related to the G(,) contributions can be computed as space distributions of the (normalized) occupancy probability, p, of water's oxygen and hydrogen atoms. If all configurations along the entire MD run are used to this purpose, the patterns reflect the symmetry of the solute. Important information concerning the interplay between solute configuration and solvent organization, as well as its dynamics, is obtained by computing such patterns separately from two different sets of segments of the MD trajectory. The two sets correspond, respectively, to the dipole pointing toward the upper or lower half-space cut by the solute plane. Patterns so obtained are shown in Figs. 5 and 6. Hydration is seen to occur preferentially around the element bearing electric charges, with a higher localization in the case of excess negative charge (solute B). The occupancy probability (normalized as specified in Simulation Details) reaches values as high as 3.5 (and, locally, even much more) around element 4, while the corresponding values in the neighborhood of the purely hydrophobic elements are in the 1.5-2.5 range. Notably, the hydration pattern even around distant hydrophobic elements is considerably altered by the addition of the excess negative charge on element 4. This agrees with the long-distance propagation and collective context-dependent nature of hydration and related free energy and forces found in simple solutes of the type used here, as well as in more realistic ones (Martorana et al., 1997, 1998; San Biagio et al., 1998).

View larger version (56K): [in this window] [in a new window]   FIGURE 5   Representation by SciAn graphics (Pepke and Lyons, 1993) of hydration isosurfaces obtained from the set of configurations corresponding to the dipole vector pointing upward for three different p values (2, 2.5, and 3). Gray, oxygen; blue, hydrogen. (Left) Solute A (uncharged). (Right) Solute B (negatively charged). Note that the distribution of the hydrogen atoms of hydration water is not faithfully rendered in this figure, as a consequence of their large vibrational disorder (see also Fig. 6).

View larger version (25K): [in this window] [in a new window]   FIGURE 6   As for Fig. 5 (note the different viewpoint). Gray clouds represent oxygen isosurfaces at p = 3 for solute A and p = 4 for solute B. The hydrogen distribution is visualized by its projection on a plane below the solute. Note that hydrogens are much closer to the charged solute element 4 than to other purely apolar solute elements. Moreover, they are more localized (higher p values) in the case of solute B.

Patterns in Figs. 5 and 6 and the just described procedure used to obtain them imply that hydration is profoundly rearranged concurrently with orientational switching events of the dipole. The dynamics of these events is evidenced by the time evolution of the  orientational coordinate of the dipole, shown in Fig. 7. Notably, the time between switching events is of the order of 100 ps, about two orders of magnitude longer than the structural relaxation time of bulk water, and of the same order of magnitude as experimentally measured structural relaxation times of hydration water (Franks, 1973). This shows that even if the reconfigurational times of the solute taken alone and those of the unperturbed solvent are individually short, their coupled dynamics can be dramatically slower.

View larger version (53K): [in this window] [in a new window]   FIGURE 7   Time evolution of the -coordinate of the dipole. Top, Solute A (uncharged); bottom, solute B (negatively charged).

In closing, it is of interest to compare the present approach to others where the use of Kirkwood's superposition approximation and expansions in terms of pair and triplet correlation functions (see, e.g., Hirata et al., 1982; Pettitt et al., 1986; Kitao et al., 1991; Klement et al., 1991; Pellegrini and Doniach, 1995; Garde et al., 1996) have allowed, for example, reproduction with reasonable accuracy hydration patterns and related free energies (Hummer et al., 1996; Garcia et al., 1997). For our comparison it is useful to remember that on a very local scale, such as that of individual protein residues, solvent-induced interactions are the result of differences of large terms, each of which is sizably affected by nonadditivity (see, e.g., Lazaridis et al., 1995). Advantages of either approach depend on the type of problem addressed and on the related, relevant scale of details. Recalling cases of either type will help our comparison. Let us first consider the case of the switching between T and R conformations of hemoglobin, which plays a crucial role in its functional oxygen transport properties. This is the first protein functional process that has been shown to be energetically dominated by changes in hydration and related free energy (Bulone et al., 1992, 1993). In this case, the conformation and related hydration changes cover the length scale of the entire protein. Hydration details are irrelevant, and inaccuracies of, say, 20% in the number of water molecules statistically involved in the hydration patterns of T and R conformations and related free energy would not affect our level of understanding. A complementary situation occurs instead, in cases where much finer details are needed, such as in functional interactions of individual protein residues. At this level of detail, work on simple models as well as on realistic systems has shown that high-order terms in Eq. 2 can even reverse the sign of forces expected to occur on individual residues (Martorana et al., 1996, 1997, 1998; San Biagio et al., 1998). In such cases, given the large free energies involved in hydration, even a mere 5% inaccuracy in hydration calculations could upset predictions concerning forces acting on individual residues. This would clearly imply all-important differences, e.g., in protein folding and interactions. It follows that detailed and costly (in terms of computer time) calculations accounting for interactions to all orders, as in the present work, are needed in all such cases.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION SIMULATION DETAILS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

For the MD studies reported in this work we have used a very simple multielement solute. The use of highly simplified solutes in realistic, explicit solvent with interactions accounted for to all orders, offers unambiguous views of the role of solvent on a detailed microscopic scale, such as that necessary to deal with specificity and recognition. These views complement and add to the coarser, large-scale perspective provided by approximate methods. Moreover, they can provide tentative "building blocks" for knowledge-based potentials. In the case of a solute like the present one, but without electric dipoles or charges, solvent-induced interactions are known to be strongly nonadditive and to propagate end to end (Martorana et al., 1997) (which is not always predictable in terms of overall additive PMF). Our solute includes charged and apolar elements with two degrees of freedom. Its configurational landscape of potential energy is flat in vacuo. However, the sole untruncated interactions with explicit molecular solvent prove to be capable of generating a rich configurational landscape of PMF and protein-like features. The depth of minima in this landscape is not far in size from free energies that stabilize protein functional conformations, notwithstanding the simplicity of our solute. The landscape visualizes solvent-induced interactions between apolar and charged groups as well as their dependence on charge sign. These findings, of course, are relevant to protein conformation and folding.

The twin observed minima reflect the symmetry of our model solute, so that much information is recovered from either of the two half-spaces. However, the actual switching of the solute + solvent system between the two minima that we have observed and discussed is a first instructive metaphor for the dynamical coupling between protein conformational switching and hydration reconfiguration. As we have seen, switching occurs on a time scale longer (by two orders of magnitude) than that of reconfigurational times of the solute taken alone or that of the unperturbed solvent. This longer time scale, however (~100 ps), is still shorter than that of actual protein switching or folding times. A relation of this finding to protein function is again illustrated by the case of human hemoglobin recalled at the end of the preceding section, because the functional switching between R and T conformations is similarly coupled to a reconfiguration of hydration water molecules. As just discussed, the related change in free energy overwhelms that of intramolecular interactions and dominates the oxygen transport function (Bulone et al., 1992, 1993; Palma et al., 1994).