Microscopic Understanding of Preferential Exclusion of Compatible

Aug 3, 2007 - An unexpected finding was that ectoine exhibits preferential binding, as an overall tendency, around both peptides. However, with the ai...
1 downloads 0 Views 306KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 10231-10238

10231

Microscopic Understanding of Preferential Exclusion of Compatible Solute Ectoine: Direct Interaction and Hydration Alteration Isseki Yu, Yoichi Jindo, and Masataka Nagaoka* Graduate School of Information Science, Nagoya UniVersity, Furo-cho, Chikusa-ku, Nagoya 464-8601, Japan ReceiVed: December 6, 2006; In Final Form: June 6, 2007

Ectoine, a zwitterionic compatible solute (CS), acts as an effective stabilizer of protein function. Using molecular dynamics simulation, solvent spatial distributions around both met-enkephalin (M-Enk) and chymotrypsin inhibitor 2 (CI2) were investigated at the molecular level in ectoine aqueous solution. An unexpected finding was that ectoine exhibits preferential binding, as an overall tendency, around both peptides. However, with the aid of the surficial Kirkwood-Buff parameter, it was clearly shown that the preferential exclusion of ectoine from the peptide surface was weaker in the smaller M-Enk than in the larger CI2. It is concluded that a denser and more structured hydration layer, such as that developed on the surface of CI2, is an important factor in the exclusion of ectoine.

1. Introduction Ectoine1,2 is a zwitterionic compatible solute (CS) that protects protein function against environmental stresses, such as salinity or high temperature. A number of experimental works have indicated that addition of CS restricts the structural fluctuation of proteins and enhances their kinetic stability.3,4 It has also been indicated experimentally that some CSs, for example, ectoine, amino acids, and sugars, increase the melting temperatures Tm of proteins, thereby enhancing the thermodynamic stabilities of their folded (native) structures.5-7 These effects have been explained by the preferential exclusion model, which states that the CS molecules are expelled from the protein surface7-9 and the growth of the preferential exclusion corresponds with the increase of excess chemical potential of the protein.7-9 Therefore, if the amount of the preferentially excluded CS molecules increases with the spread of the solventaccessible surface area (SASA) of the protein, the added CS functions to make the protein conformation more compact in order to reduce the amount of such preferentially excluded CS molecules. Using the molecular dynamics (MD) simulation method, we have demonstrated that ectoine molecules do not strongly change the hydration layer near the surface of a target protein, chymotrypsin inhibitor 2 (CI2), in high-temperature (370 K) ectoine aqueous solution.10 In fact, our MD simulation indicated numerically that ectoine molecules are preferentially excluded near the CI2 surface. Further, we found that zwitterionic ectoine molecules strongly interact with water molecules via their hydrogen bondings11 and slow solvent diffusion in the bulk phase. It was, therefore, concluded that water molecules near the surface of proteins significantly slow their diffusion to the bulk phase.10 The slowdown brings about the restriction of structural fluctuation of the protein. We thus found one of the microscopic origins of the kinetic stabilization of proteins in aqueous solution of CS, which is in accordance with the preferential exclusion model. However, the problems of how CS molecules interact microscopically with biomacromolecules * To whom correspondence should be addressed. E-mail: mnagaoka@ is.nagoya-u.ac.jp. URL: http://frontier.ncube.human.nagoya-u.ac.jp/. Tel & Fax: +81-52-789-5623.

and whether the addition of CS might indirectly stabilize them, irrespective of their molecular properties, remain unresolved. In the present study, we performed long-time full atomic MD simulations of more than 150 ns for a small peptide, metenkephalin (M-Enk), in addition to a larger protein, CI2, in ectoine aqueous solution of the same concentration and in pure water at room temperature. To determine the spatial distribution of each solvent component, the atom number densities of water molecules and ectoine molecules around both solutes were analyzed. Furthermore, to specify the spatial difference in the degree of preferential exclusion, the preferential exclusion parameters of ectoine were defined and analyzed for both solutes with the aid of the Kirkwood-Buff (KB) theory. To obtain the spatial profiles of the preferential exclusion parameter, we modified the calculation procedure of the KB parameter from that in the original KB theory so as to enable systematic discussion of the dependence of the KB parameter on the distance from the molecular surface of a solute molecule. We found that one dominant structure of M-Enk does not preferentially exclude the ectoine molecules from its surface as CI2 does. To understand the reason for this difference in ectoine exclusion, in addition to the effect of direct interaction between M-Enk and ectoine, the influence of hydration (that is, property alteration of the hydration layer near the solute surface) on the development of ectoine preferential exclusion around each solute was examined at the molecular level. This Article is organized as follows. First, in the following section, the Computational Methods and the model systems are explained. Numerical results and discussion are provided precisely in the Results and Discussion section. Finally, the present study is summarized in Concluding Remarks. 2. Computational Methods 2.1. Molecular Dynamics Simulations. All molecular dynamics (MD) calculations were performed using the AMBER7 program.12 The force field parameter set, parm99,13 was used for all of the molecules in the present systems. The partial atomic charges of ectoine molecules were taken from our previous study.11 The starting conformation of the zwitterionic form of M-Enk that we adopted was an X-ray structure in which the peptide backbone is extended.14

10.1021/jp068367z CCC: $37.00 © 2007 American Chemical Society Published on Web 08/03/2007

10232 J. Phys. Chem. B, Vol. 111, No. 34, 2007

Yu et al.

σRs(ri) ) min |ri - rRk |

TABLE 1: Numbers of Solvent Molecules, Simulation Length, and Concentration of Ectoine in Each System number of number of ectoine simulation concentration of water ectoine (M) system solute molecules molecules length (ns) MEmL MEmH MEp CI2m CI2p

M-Enk M-Enk M-Enk CI2 CI2

1095 1072 1178 6903 8374

36 59 0 227 0

50 50 50 2 2

〈

〉 〈

δNRs(r,t) δVR(r,t)

)

〉

NRs(r + ∆r/2, t) - NRs(r - ∆r/2, t)

VR(r + ∆r/2, t) - VR(r - ∆r/2, t) (1)

where 〈‚‚‚〉 denotes the time average defined as 1/T (‚‚‚)dt and the instantaneous integrated coordination number NRs(r,t) is defined as10 ∫T0

Ms

NRs(r,t) )

h(r - σRs(ri(t))) ∑ i)1

Here, h(x) is a step function, as follows

h(x) )

1.5 2.4 0.0 1.5 0.0

MD simulation of M-Enk in 1.5 M ectoine aqueous solution was performed as follows. First, the M-Enk was set in a periodic boundary box (37.7 × 34.2 × 31.3 Å3) filled with TIP3P water molecules15 and ectoine molecules (this system was named MEmL). Second, other solvation models of M-Enk in 2.4 M ectoine aqueous solution and pure water (that is, MEmH and MEp) were independently constructed. Ectoine concentrations in both MEmL and MEmH were within that generally used for experimental works of ectoine-type compatible solutes.5 After the starting structures were minimized molecular mechanically for 2000 cycles to reduce any bad contact pairs, each system was initially equilibrated for the first 200 ps under the NPT condition at 300 K and 1 atm. With these equilibration procedures, pressure, temperature, total energy, and the density in each system were saturated at constant values. Then, further simulation of 50 ns was performed for each system at 300 K under the NVT condition. Trajectories were numerically integrated by the Verlet method with a time step of 2.0 fs. The electrostatic interactions were treated by the particle-mesh Ewald (PME) method.16 All of the bonds involving hydrogen atoms were constrained by the SHAKE method.17 For equilibration, both temperature and pressure were regulated using the Berendsen algorithm.18 To investigate the difference in solvent distributions, MD simulation for CI2 in pure water (CI2p) and in the same solvent condition as that of MEmL (CI2m) were performed for 2 ns at 300 K following the same procedure as that in our previous simulation of CI2.10 The numbers of solvent molecules are summarized in Table 1. We adopted 1.5 M for the present ectoine concentration with reference to the experimental work with hydroxyectoine in aqueous solution.5 The weight molarity of ectoine in CI2m was 1.8 mol/kg (the solution density was 1.05 g/cm3), while that in the experiment of 1.5 M hydroxyectoine was 1.4 mol/kg (the solution density was 1.3 g/cm3). 2.2. Solvent Atom Number Density. In the following explanation, the subscript s of quantities takes w or e in correspondence to water or ectoine for the solvent component, while R takes m (M-Enk) or c (CI2) for the solute molecule. Then, the solVent atom number density (SAND) FRs(r) of a certain solvent component s is defined as a function of the distance r from the solvent molecule R (see Figure 1) by

FRs(r) )

(3)

k∈R

(2)

{

1 (x g 0), 0 (x < 0)

(4)

Then, h(r - σRs(ri(t))) takes 1 if the constituent atom i of the solvent component s is included within the region υR(r,t) (the interior of the red dotted closed curve in Figure 1), that is, if r - σRs(ri) > 0, where σRs(ri) is defined as the minimum distance between the atom i of the solvent component s and any solute atom k of the solute molecule R. Thus, NRs(r,t) is such an integral number that counts all of the constituent atoms of the solvent component s that are included within the region υR(r,t), whose volume and surface area are denoted by VR(r,t) and SR(r,t), respectively. Ms denotes the total number of constituent atoms in the solvent component s (for example, Mw in CI2m is 20709). δυR(r,t) indicates the region wedged between two regions, υR(r + ∆r/2,t) and υR(r - ∆r/2,t) (gray band in Figure 1). The volume and the number of atoms of the solvent component s that are contained in δυR(r,t) are denoted by δVR(r,t) and δNRs(r,t), respectively. 2.3. Excess Coordination Number and Preferential Exclusion Parameter. The excess coordination number of atoms in solvent component s around a certain atom site k in the solute ex ” is defined as follows19,20 molecule R “NRs ex ) Fs(∞)GRs ) Fs(∞) NRs

∫ [cks(r) - 1]dr

(5)

where cks(r) is the correlation function between the site k and the atoms in s when they are separated by r, and Fs(∞) is the atom number density of solvent component s in the bulk solvent phase. GRs is often referred to as the Kirkwood-Buff (KB) parameter.19,20 Originally, GRs was defined by GRks(R) in the infinity limit of the radius R that denotes a large distance including the whole inhomogeneous phase around the solute molecule R

GRs ) lim GRks(R)

(6a)

Rf∞

where the R-dependent KB parameter is defined as

GRks(R) )

∫0R [gks(r) - 1]4πr2dr ) R 1 4πR3 Fs(∞) ∫0 gks(r)4πr2dr Fs(∞) 3 F (∞) s

{

}

(6b)

in the use of gks(r), that is, the radial distribution function (RDF) between the atom k and the atoms in the solvent component s. GRks(R) is hereafter called the radial KB parameter. Unlike the r used in subsection 2.2, r in eq 6b indicates the distance from the center of a certain atom k in the solute molecule R. GRks(R) multiplied by Fs(∞), that is, Fs(∞)GRks(R), should be related to ex , that is, the excess coordination number of atoms in the NRs solvent component s around the solute molecule R ex lim Fs(∞)GRks(R) ) Fs(∞)GRs ) NRs Rf∞

(7)

where GRs ()GRks(∞)) indicates the KB parameter of the solvent component s around the solute molecule R. However, to include the whole inhomogeneous phase, there is a practical difficulty in that R might exceed the half of the simulation box length used in our model systems. Moreover, if the profile of the radial KB parameter (eq 6b) depended

Preferential Exclusion of Compatible Solute Ectoine

J. Phys. Chem. B, Vol. 111, No. 34, 2007 10233

Figure 1. Schematic representation of the region υR(r,t) and key quantities. The distance r from (or between) certain atom(s) is that from (or between) the center of the atom(s). The (center of) constituent atoms ks and their bonds in the part of the solute molecule R are represented as black points and black lines. The whole shape of the solute molecule R is also represented as the right amoeba-like shape painted in blue. The boundary of υR(r,t) is defined as a surface whose minimum distance to any solute atom k of the solute molecule R is r and is represented as the red dashed curve. The volume and surface area of υR(r,t) are denoted by VR(r,t) and SR(r,t), respectively. The region wedged between the two regions υR(r + ∆r/2,t) and υR(r - ∆r/2,t), that is, δυR(r,t), is represented as a gray band. The volume of δυR(r,t) and the number of those atoms of the solvent component s contained in δυR(r,t) are denoted by δVR(r,t) and δNRs(r,t), respectively.

sensitively on the location of the atom site k selected arbitrarily in the solute molecule R, its direct estimation might not be suitable for discussion of the spatial difference between those solvent distributions around two kinds of solute molecules with significant differences of size, conformation, and surface topography, such as those between CI2 and M-Enk. Hence, in the present study, we have defined another type of R-dependent S (R), that is, the surficial KB parameter, KB parameter, GRs which enables us to distinguish the profiles of excess coordinaex tion numbers NRs (R) as a function of the distance R, not from a certain atom site in the solute molecule R but from the surface of the solute molecule R, which can be defined by a collection of the solute R atoms closest to each solvent atom S GRs (R) )

∫r)0

r)R

〈(

) 〉

1 δNRs(r,t) - 1 SR(r,t) dr (8a) Fs(∞) δVR(r,t)

Here, r denotes the distance from the solvent molecule R in comex mon with Figure 1. Thus, the NRs (R) profile can be obtained by the coordination number NRs(R,t) in the region υR(R,t), and its volume VR(R,t) defined in subsection 2.2, is as follows ex S (R) ) Fs(∞)GRs (R) ) 〈NRs(R,t) - Fs(∞)VR(R,t)〉 NRs

(8b)

By this definition, one can clearly understand the inhomogeneous profile of the R-dependent KB parameter and the excess coordination number of the solvent species, which reflect the spatial/ structural characteristics around/of each individual solute molecule. According to the KB theory,19,20 the preferential exclusion parameter of ectoine can be defined as

νRe(R) )

Fe(∞){GSRe(R)

ex where Nex Re(R) and NRw(R) are the excess coordination numbers of atoms in ectoine molecules and water molecules associated with the solute molecule R, respectively. A negatiVe value in νRe(R) means that the preferential exclusion of ectoine occurs until the distance R from the surface of the solute molecule R, while a positiVe value in νRe(R) means that the preferential binding of ectoine occurs. Although νRe(R) is referred to as a preferential binding parameter21 in some cases, in the present work, it is referred to entirely as the preferential exclusion parameter. 2.4. Hydrogen Bond Density in the First Hydration Layer. We define also the hydrogen bond density (HBD) FHB R in the first hydration layer of a solute molecule R as follows

-

S GRw (R)}

)

〈 〉 nHB R (t) AR(t)

Fe(∞) Fw(∞)

ex (R) (9) NRw

(10)

where nHB R (t) is the total number of hydrogen bonds (HBs) formed among water molecules in the first hydration layer of a solute molecule R at time t. AR(t) is the SASA of a solute molecule R at time t. The cutoff distance and the angle for the identification of HB were set to 3.5 Å and 60°, respectively. In the present work, a water molecule ω in the first hydration layer is defined such that its oxygen atom j exists within 3.0 Å from the closest non-hydrogen atom in the solute molecule R as

ω ∈ {w | σRw(rj(t)) < 3.0Å}

(11a)

σRw(rj) ) min |rωj - rRl |

(11b)

where l∈R,l*H

)

Nex Re(R) -

FHB R

is the minimum distance between the oxygen atom j of the water molecule ω and any atom l of the solute molecule R. It should be noted that the density FHB R defined by eq 10 is the number of

10234 J. Phys. Chem. B, Vol. 111, No. 34, 2007

Yu et al.

Figure 2. The end-to-end distance d in (a) MEp, (b) MEmL, and (c) MEmH.

Figure 3. Radial distribution functions between one of the two oxygen atoms in the ectoine COO group (OE1) and backbone amide nitrogen atoms of M-Enkb during the bent period in MEmL (NT, the amide nitrogen in the N terminus; Nn (n ) 2-5), the amide nitrogen in the nth residue). A typical snapshot of an ectoine molecule binding to M-Enkb is depicted in the inset (M-Enkb: backbone (red tube), amide group (orange ball and stick), and side chain (red line); ectoine: carbon (light blue), nitrogen (blue), oxygen in COO group (red), and hydrogen (white)).

HBs per unit area, while the density FRs(r) defined by eq 1 is the number of atoms per unit Volume. 3. Results and Discussion 3.1. Structural Fluctuation Restriction of Met-Enkephalin by Ectoine Addition. To investigate the stability of the peptide backbone dynamics from the structural viewpoint, the end-toend distance d of M-Enk was introduced as an order parameter and monitored independently in pure water and 1.5 and 2.4 M ectoine aqueous solutions. Figure 2 shows the time dependences of d(CR1-CR5) in each solution (the residue number was assigned serially from the N terminus). It was seen that most of the ordered M-Enk conformations with d(CR1-CR5) ≈ 5 Å are bent structures, which are denoted collectively as M-Enkb. Hereafter, those time intervals during which M-Enk preferentially adopts such a bent structure, that is, M-Enkb, are called bent periods. Further analysis showed that M-Enk adopted various conformations in all of the solutions. However, during the bent periods, for example, 20-30 ns in MEp, 8-34 ns in MEmL, and 28-50 ns in MEmH, the peptide backbone structure of M-Enkb was maintained essentially in the same conformation, that is, a round bent conformation forming a salt bridge between

the oppositely charged N and C termini (see the inset in Figure 3). In fact, in previous MD simulations, enkephalins also exhibited a preference for the bent structure in pure water and in DMSO.22,23 Despite these common M-Enk behaviors observed in all three solution systems, it was found that there is a significant change in the dynamics of M-Enk in ectoine aqueous solutions. Namely, during the bent periods in these solutions, M-Enkb remains more stably in the backbone conformation with less fluctuation (Figure 2b and c). In contrast, in pure water, the fluctuation is considerably larger in the bent periods (Figure 2a). It is true that our previous study showed that the ectoine addition significantly and directly slows the solvent diffusion.10 However, that effect is also likely an important factor to indirectly restrict the M-Enk dynamics itself. In fact, in the denser 2.4 M ectoine solution, that is, MEmH, it was recognized that M-Enk required a longer time (about 30 ns) to transition from the initial stretched structure into the bent one (Figure 2c). It should be noted that these results are essentially consistent with current experimental work on protein dynamics in aqueous solutions with CSs.3,4,24 3.2. Stabilization of Met-Enkephalin by Direct Interaction with Ectoine. In addition to the indirect effect of ectoine addition, that is, the slowdown of solVent diffusion, our analysis of the direct interaction between M-Enk and ectoine molecules implies that one should take into consideration the influence of direct interactions between a solute and ectoine molecules not only on the coordination of ectoine but also on the dynamics of such small biomacromolecules as M-Enk. To investigate the M-Enk stabilization in ectoine aqueous solutions from the microscopic viewpoint, the direct interaction between M-Enkb and ectoine molecules was analyzed in the bent periods of MEmL (Figure 2b). From this analysis, we found a characteristic coordination of an ectoine molecule to M-Enkb that is likely to reduce the structural fluctuation of the solute. In this coordination (see inset in Figure 3), at least one ectoine molecule was bound for a long time (more than several nanoseconds) to the center of the bent peptide backbone of M-Enkb since almost all of the backbone amide hydrogen atoms interacted simultaneously with the ectoine COO group via the electrostatic interaction. This ectoine binding is probably one of the microscopic origins restricting the structural fluctuation of the bent structure. Figure 3 shows the profiles of radial distribution functions (RDFs) of ectoine oxygen around amide nitrogen in M-Enkb (we utilized 11-20 ns in MEmL for the analysis). The large and the sharp first peaks in the RDFs indicate a strong electrostatic interaction between positively charged backbone amide hydrogen atoms and the oxygen atoms in the COO group of ectoine. This result attracted our interest since it is well-known phenomenologically that zwitterionic CS aqueous solutions do

Preferential Exclusion of Compatible Solute Ectoine

J. Phys. Chem. B, Vol. 111, No. 34, 2007 10235 during the CI2m simulation, the entire time duration of the CI2m simulation was utilized for the analysis. In contrast, because the structure of M-Enk changed clearly, just a part of the bent period (11 ∼ 20 ns) in the MEmL simulation was utilized (Figure 4). δNRs(r,t) was then calculated from 1000 evenly extracted snapshots for CI2 and from 10000 snapshots for M-Enkb. Calculation of δVR(r,t) for such numbers of snapshots was computationally demanding. However, according to our analysis, the variance (the relative standard deviation) of δVR(r,t) for each solute peptide was found to be less than 3% at any distance r. In other words, it was understood that both δVR(r,t)s were almost constant during the utilized time durations. In order to approximate FRs(r), therefore, δNRs(r,t) of each snapshot was divided by the average value of δVR(r,t)s, that is, δVR(r)s calculated from 10 CI2 snapshots and from 100 M-Enkb snapshots, respectively

Figure 4. Solvent atom number densities (Fcw(r), blue solid; Fmw(r), blue dotted; Fce(r), red solid; and Fme(r), red dotted). The entire time duration of the CI2m simulation and a part of the bent period (11 ∼ 20 ns) in the MEmL simulation was utilized for the analysis.

not favor the exposure of the peptide backbone unit.25 However, from a further finding that the direct interaction is significantly reduced in the other structures of M-Enk in both MEmL and MEmH, we understood that this ectoine binding is induced by the characteristic conformation of M-Enk, in which the close gathering of multiple backbone amide hydrogen atoms results in the concentration of their positive charges. This is consistent with the analytical model proposed by Bolen and co-workers that shows a correlation between the polar surface area of osmolytes (compatible solutes) and their favorable interaction with the peptide backbone.26 Moreover, a similar interaction was previously indicated theoretically in the MD simulation of bis(enkephalin) in 1.0 M sodium acetate aqueous solution.27 In this case as well, the simultaneous binding of an acetate anion to backbone amide groups influenced the dynamics of the peptide conformation.27 These results indicate that an ectoine molecule could behave as an ionic rather than a zwitterionic species when the solute peptide forms a charge-concentrated surface as M-Enkb does. Especially for a smaller peptide such as M-Enk, this electrostatic effect likely becomes more effective if, as may be expected, the influence from the other residues is weaker than that from residues in a larger protein such as CI2. Although our MD simulation indicated that the ectoine addition also reduces the structural fluctuation of M-Enk, as we explained above, the mechanism is different from that in the case of CI2 where the slowdown of the solVent diffusion is the main contribution. This was one of our motivations in the present study to investigate quantitatively how and to what extent the exclusion of ectoine should be developed around not only CI2 but also M-Enk (see subsections 3.3 and 3.4). We checked the reproducibility of the following key aspects. First, several simulations at 370 K were previously executed (data not shown), and the restriction of the transition between extended and bent structures of the M-Enk with the addition of ectoine was reproduced. Second, the same characteristic ectoine binding shown in Figure 3 was also frequently observed. 3.3. Difference between the Solvent Distributions around Met-Enkephalin and Chymotrypsin Inhibitor 2. Solvent atom number densities (SANDs) FRs(r) (eq 1) are shown in Figure 4. To obtain their numerical values, ∆r was set to 0.05 Å, and δVR(r,t) was then calculated with the resolution of 0.008 Å3. Because CI2 did not show any significant structural change

FRs(r) =

〈

〉

δNRs(r,t) δVR(r)

(12)

In Figure 4, the profiles of FRs(r) are found to show reasonable plateau values from r ≈ 6.0 Å, indicating that the ectoine and water distributions reach their bulk values. It is significant that Fmw(r) is smaller than Fcw(r) around the first and second peaks, and this fact indicates that M-Enkb does not have such a dense hydration layer as CI2. In contrast, the larger intensity of Fme(r) until r ≈ 5 Å compared with that of Fce(r) indicates a denser accumulation of ectoine molecules near the surface of M-Enkb. 3.4. Development of the Preferential Exclusion near the Surface of Chymotrypsin Inhibitor 2. Here, we discuss the preferential exclusion of ectoine using both the surficial KB S parameters GRs (R) (eq 8a) and the preferential exclusion parameter VRe(R) (eq 9) given by the KB theory.19,20 The time durations utilized for the analysis are the same as those in subsection 3.3. The values of several parameters in the KB theory defined originally at infinite distance from the solute molecule R were approximately substituted at R ) 10.0 Å and are presented in Table 2 with the average SASA of both solutes. S Hereafter, assuming that 10.0 Å is large enough, GRs (R) (R ) 10.0 Å) will be identified with the original KB parameter GRs. First, let us remember again that KB parameters GRw or GRe itself are related to the extent to which the individual solvent component (water or ectoine) is excluded or bound around a solute molecule R. It is understood that the negative values of both Gmw and Gcw in Table 2 indicate that the exclusion of water molecules mainly comes from the van der Waals volume. The present numerical value of Gcw()-5827.35 cm3/mol) could be reasonably compared to that of a similar size of protein in “pure” water obtained in a previous theoretical work;28 Imai et al. developed a theoretical treatment of the KB theory combined with 3D-RISM theory and calculated the partial molar volume for several proteins in pure water.28 According to their results, the KB parameter of water for ubiquitin, which has almost the same SASA and slightly larger amino acid chain length as CI2, was evaluated as -6081 cm3/mol. Although it is meaningless to directly compare the values of the KB parameter of two solutes with different molecular sizes, we would like to emphasize that the opposite sign of KB parameters means a clear opposite tendency of the ectoine distributions around two solutes, that is, the binding for M-Enkb (positive value in Gme) and the exclusion for CI2 (negative value in Gce). Next, we further discuss the preferential exclusion of ectoine using the preferential exclusion parameter νRe (eq 9). Since the positive value of νRe indicates the preferential binding of ectoine

10236 J. Phys. Chem. B, Vol. 111, No. 34, 2007

Yu et al.

TABLE 2: Key Quantities Indicating the Difference between the Solvent Distributions around M-Enkb and CI2a R

AR (Å2)

GRw (cm3/mol)

GRe (cm3/mol)

ex NRw

ex NRe

νRe

νRe/AR (Å-2)

m (M-Enkb) c (CI2)

880.67 4769.12

-728.45 -5827.35

612.45 -4614.92

-100.94 -805.70

17.91 -147.19

39.22 38.67

0.046 0.0081

ex Symbols. AR: average SASA of the solute molecule R. GRs: KB parameter of the solvent component s around the solute molecule R. NRs : excess coordination number of atoms in the solvent component s around the R. VRe: preferential exclusion parameter of ectoine around the R. ex S VRe/AR: preferential exclusion parameter of ectoine per unit SASA of the R. (GRs, NRs , and VRe were evaluated by the surficial KB parameter GRs (R) (eq 8a) using R ) 10.0 Å.) SASAs were calculated with a probe sphere having a radius of 1.4 Å. a

Figure 5. Differences of the preferential exclusion parameter of ectoine; (a) δνme(r) around M-Enkb and (b) δνce(r) around CI2.

as νme around M-Enkb, in the case of CI2, despite the local exclusion of ectoine around CI2 shown by Gce, the present positive νce value indicates that ectoine molecules do not show the preferential exclusion but the preferential binding as an overall tendency in the CI2m simulation. However, it could be said that a considerable difference exists between two solutes in the characteristics of ectoine exclusion. For instance, while νRes become very similar positive values for both solutes, their value per unit SASA for CI2, that is, νce/Ac, is about 6 times smaller than that for M-Enkb (see the last column in Table 2). This fact is clear evidence that the surface of CI2 excludes ectoine molecules somewhat more strongly than does the surface of M-Enkb. Furthermore, another significant difference was found in the ectoine exclusion style between those two νRe(R) profiles. In Figure 5, δνRe(r), the difference of νRe(R) defined as

δνRe(r) ) νRe(r + ∆r/2) - νRe(r - ∆r/2)

(13)

is shown, where a negative and a positive value in δνRe(r) mean that the preferential exclusion (black) and the preferential binding (gray) of ectoine occur at the distance r from the surface of solute molecule R, respectively. The large black area in Figure 5b shows that the preferential exclusion of ectoine is developed near the surface of CI2. On the other hand, δνme(r) shows the preferential binding of ectoine at almost any distance from the surface of M-Enkb to the bulk solvent phase (Figure 5a). From these analyses, it seems probable that a stronger preferential exclusion of ectoine is developed around CI2 than around M-Enkb. It was, therefore, unexpected that the computational result, that is, νRe values, showed that not only M-Enk but also CI2 exhibited preferential binding of ectoine, considering the generally believed functions of the CSs on proteins, that is, preferential exclusion. We therefore discuss the factors weakening the preferential exclusion of ectoine, focusing on the molecular force field used in the present study. Smith et al. previously compared the KB parameters of several other solutes (urea, acetone, amides, and NaCl) in

water29-32 with experimental data, and it is well-known that when some generally used force fields are applied to these solutes, they tend to show too much self-association in TIP3P and SPC33 water. At the same time, in most of these cases, the KB parameters of water molecules around the cosolvent molecules are considerably smaller than the experimental values. The excessive self-association of cosolvent molecules is attributed to their lack of affinity for water molecules. Smith et al. reported that many solute models using ab initio gas-phasederived charge distributions might not accurately represent polarized solution charge distributions.29,30 These facts may be one reason that excessive association between the ectoine molecules and the solute seemed to occur in our simulation, in which the partial atomic charges of ectoine were optimized by ab initio methods in the gas phase.11 In other words, these findings suggest that, in reality, ectoine molecules may accumulate water molecules more strongly, making access of ectoine molecules to the protein surface more difficult. One may expect, therefore, that ectoine molecules with more appropriate charge distributions would increase the water density a little more than what is estimated by present MD simulations and reduce the νce value reasonably so as to reproduce the preferential exclusion of ectoine for CI2 (cf. Equation 9). 3.5. Difference between the Property of Hydration Layer around Met-Enkephalin and that around Chymotrypsin Inhibitor 2. According to the analysis of surficial KB parameters, it was clearly shown that preferential exclusion develops near the surface of CI2 while it does not develop around M-Enkb. As shown in Figure 3, there are strong electrostatic interactions between polar backbone amide hydrogen atoms of the solute and the ectoine COO groups in the case of the bent structure of M-Enk in MEmL. In addition to the presence of such characteristic interactions, we suggest the hydration difference is also important in inducing the change in the ectoine preferential exclusion. In fact, as summarized in Table 3, with the amino acid chain length nres R , total number of constituent atoms nR, average SASA AR, and the average number of HB

Preferential Exclusion of Compatible Solute Ectoine

J. Phys. Chem. B, Vol. 111, No. 34, 2007 10237

Figure 6. Hydrogen-bonding (HB) (red dotted lines) networks of water molecules (blue triangles) in the first hydration layer of (a) M-Enkb and (b) CI2 in pure water.

TABLE 3: Key Quantities Related to the Sizes of Solute Molecules M-Enkb and CI2 and the Hydrogen-Bonding (HB) Networks around Them in the Pure Watera R

nres R

nR

AR (Å2)

nHB R

3 FHB R × 10 (Å-2)

m (M-Enkb) c (CI2)

5 65

75 1076

880.67 4769.12

2.2 36.4

2.5 ( 1.5 7.6 ( 1.2

a

Symbols. nres R : amino acid chain length. nR: total number of

constituent atoms. AR: average SASA. nHB R : average number of hydrogen bonds. FHB : hydrogen bond density for the solute molecule R R. HB nHB R for each solute molecule R in the pure water, the Fm was HB found to be about three times smaller than Fc . Snapshots of HB networks around each solute in pure water are shown in Figure 6a and b. The loose water HB network around M-Enkb should make the penetration of ectoine molecules in the layer easier. We suggest that the structural difference, that is, the molecular size difference, between CI2 and M-Enkb is one of the important factors that changes the properties of their hydration layer. Because the first hydration layer of the smaller M-Enk needs much larger curvature to encapsulate it, M-Enkb cannot form as dense a layer as that of CI2 (see Fmw(r) and Fcw(r) in Figure 4). This inversely increases the ratio of ectoine atoms in the layer. Further, the increase in the curvature of the first hydration layer should also reduce some tension in the water HB networks wrapping around the solute peptides. To investigate this, we analyzed the HB networks in the first hydration layer of the bent structures of M-Enkb and CI2. To eliminate the influence of ectoine addition, we investigated them in pure water systems, that is, MEp and CI2p (Table 3). As was defined by eq 10, the HBDs’ FHB R s for each solute were obtained by averaging over

five different snapshots from CI2p and the bent period (21∼30 ns) of MEp. According to the scaled particle theory, the surface tension on the hard sphere cavities, created in the liquid molecules, would increase with the cavity radius until around 10 Å.34 This is consistent with the fact that CI2 has a denser and more structured first hydration layer than that of M-Enk (with the average radius ∼8 Å). Finally, let us discuss the influence of the “surface” topography on the preferential exclusion of ectoine. Because almost all of the atoms in M-Enk are highly exposed to the solvent phase even when M-Enk is in its bent structure, M-Enk does not have as many internal surfaces (cavities) as CI2 might have. If these cavities could hold water molecules and would exclude only ectoine molecules, the presence of such cavities should enhance the preferential exclusion resulting from the size effect, that is, the cavity effect. However, according to the observation that CI2 did not have any significant cavities during the simulations, it was clarified that water molecules are also excluded from the interior of CI2, implying the above cavity effect is insignificant for CI2. Next, we discuss, further, the total area of the external surface of the protein that might be related to enhance the ectoine preferential exclusion. According to our calculations, among the whole SASA of CI2, which was calculated with a probe sphere having a radius of 1.4 Å, it was found that only about 300 Å2 of the area cannot be accessed by a probe sphere with a radius of 4.0 Å (that is, the size of an ectoine molecule) and is nothing but 6% of the whole SASA of CI2. We conclude, therefore, that the effect of surface topography must be less important than the other two on the preferential exclusion when CI2 is in its native structure. 4. Concluding Remarks It is generally considered that CSs are preferentially excluded from the surface of proteins.7-9 Similarly, our MD simulations

10238 J. Phys. Chem. B, Vol. 111, No. 34, 2007 of CI2 also showed the preferential exclusion of ectoine near its surface. On the other hand, the preferential exclusion of ectoine was found to be significantly weakened near the surface of the much smaller solute M-Enk in its bent structure. To obtain the spatial profiles for the preferential exclusion of ectoine, we modified the calculation procedure of the KB parameter from that in the original KB theory to define a new R-dependent KB parameter, that is, the surficial KB parameter, and systematically examined the dependence of the KB parameter on the distance from the molecular surface of a solute molecule. Accordingly, the origin of this weakening of the preferential exclusion around M-Enkb can be attributed to the following three factors: (i) the characteristic strong interaction between polar backbone amide hydrogen atoms and ectoine COO groups, (ii) the lack of a dense and structured hydration layer, which is partially due to the smallness of M-Enk, and (iii) the lack of the “surfaces” that enhance the ectoine exclusion by the steric exclusion. Of these three factors, (i) and (ii) must be coupled and enhance each other to result in the weakening of preferential exclusion. As for factor (iii), although we first considered the steric exclusion to be an essential factor, it might not be so important to increase the amount of preferential exclusion of ectoine around CI2. In summary, it is concluded that the presence of a dense and structured hydration layer, which should be developed on a protein surface such as that of CI2, is an important microscopic factor enhancing the preferential exclusion of zwitterionic CSs. At the same time, it is also suggested that the influence of the direct interaction between a smaller peptide and a zwitterionic CS molecule, which might be induced by the charge-concentrated conformation of the peptide, should be more seriously considered. In the present microscopic investigation with the aid of the surficial KB parameter, it was shown that the preferential binding of ectoine occurs around M-Enkb, while the preferential exclusion occurs clearly in the immediate vicinity of CI2. However, the preferential exclusion parameters νRe were themselves found to be positive, showing preferential binding not only around M-Enkb but also around CI2. Although these results seem to be contradictory, they may also be viewed as evidence that there is still room for improvement of the molecular force field in the estimation of νRe with the use of molecular dynamics simulations. We suspect that the electrostatic interaction between the ectoine molecules and the water molecules plays an essential role in the distributions of ectoine molecules around proteins. It may be conjectured that improvement of the partial atomic charge distribution of ectoine and alteration of the water force field, for example, from TIP3P to some other water model, would bring about an increase in the water accumulation around the ectoine molecules and would contribute to increasing Gew and decreasing GRe. This would improve the preferential exclusion parameters νRe so that they are in much better agreement with the experimental ones. Acknowledgment. This work was supported partly by Grantin-Aids for Scientific Research and for Young Scientists (B) (1975004) from the Ministry of Education, Culture, Sport, Science and Technology in Japan and also by a Grant-in-Aid

Yu et al. for the 21st Century COE program “Frontiers in Computational Science” at Nagoya University and the Core Research for Evolutional Science and Technology (CREST) “High Performance Computing for Multi-scale and Multi-physics Phenomena” from the Japan Science and Technology Agency. References and Notes (1) Galinski, E. A. Experientia 1993, 49, 487-496. (2) Galinski, E. A.; Pfeiffer, H. P.; Truper, H. G. Eur. J. Biochem. 1985, 149, 135-139. (3) Lamosa, P.; Turner, D. L.; Ventura, R.; Maycock, C.; Santos, H. Eur. J. Biochem. 2003, 270, 4604-4614. (4) Foord, R. L.; Leatherbarrow, R. J. Biochemistry 1998, 37, 29692978. (5) Knapp, S.; Ladenstein, R.; Galinski, E. A. Extremophiles 1999, 3, 191-198. (6) Yancey, P. H.; Clark, M. E.; Hand, S. C.; Bowlus, R. D.; Somero, G. N. Science 1982, 217, 1214-1222. (7) Lee, J. C.; Timasheff, S. N. J. Biol. Chem. 1981, 256, 7193-7201. (8) Arakawa, T.; Timasheff, S. N. Arch. Biochem. Biophys 1983, 224, 169-177. (9) Arakawa, T.; Timasheff, S. N. Biophys. J. 1985, 47, 411-414. (10) Yu, I.; Nagaoka, M. Chem. Phys. Lett. 2004, 388, 316-321. (11) Suenobu, K.; Nagaoka, M.; Yamabe, T.; Nagata, S. J. Phys. Chem. A 1998, 102, 7505-7511. (12) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohlke, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman P. A. AMBER7; University of California: San Francisco, CA, 2002. (13) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049-1074. (14) Griffin, J.; Langs, D.; Smith, G.; Blundell, T.; Tickle, I.; Bedarkar, S. Proc. Natl. Acad. Sci. U.S.A. 1986, 83, 3272-3276. (15) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (16) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 1008910092. (17) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327-341. (18) Berendsen, H. J. C.; Postma, J. P. M.; van Gusteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (19) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774-777. (20) Ben-Naim, A. Statistical Thermodynamics for Chemists and Biochemists; Plenum Press: New York, 1992. (21) Aburi, M.; Smith, P. E. J. Phys. Chem. B 2004, 108, 7382-7388. (22) Spoel, D. V.; Berendsen, H. J. C. Biophys. J. 1997, 72, 2032-2041. (23) Shen, M.; Freed, K. F. Biophys. J. 2002, 82, 1791-1808. (24) Jacob, M.; Schindler, T.; Balbach, J.; Schmid, F. X. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 5622-5627. (25) Auton, M.; Bolen, D. W. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 15065-15068. (26) Street, T. O.; Bolen, D. W.; Rose, G. D. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 13997-14002. (27) Smith, P. E.; Marlow, G. E.; Pettitt, B. M. J. Am. Chem. Soc. 1993, 115, 7493-7498. (28) Imai, T.; Kovalenko, A.; Hirata, F. J. Phys. Chem. B 2005, 109, 6658-6665. (29) Weerasinghe, S.; Smith, P. E. J. Phys. Chem. B 2003, 107, 38913898. (30) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 118, 1066310670. (31) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 119, 1134211349. (32) Myungshim, K.; Smith, P. E. J. Comput. Chem. 2006, 27, 14771485. (33) Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (34) Ka, L.; Chandler, D.; Weeks, J. D. J. Phys. Chem. B 1999, 103, 4570-4577.