Selective Ion Binding by Protein Probed with the Statistical Mechanical

4588

J. Phys. Chem. B 2007, 111, 4588-4595

Selective Ion Binding by Protein Probed with the Statistical Mechanical Integral Equation Theory Norio Yoshida,† Saree Phongphanphanee,‡ and Fumio Hirata*,†,‡ Department of Theoretical Molecular Science, Institute for Molecular Science, Okazaki 444-8585, Japan, and Department of Functional Molecular Science, The Graduate UniVersity for AdVanced Studies, Okazaki 444-8585, Japan ReceiVed: December 13, 2006; In Final Form: February 18, 2007

Selective ion binding by human lysozyme and its mutants is probed with the three-dimensional interaction site model theory which is the statistical mechanical integral equation theory. Preliminary and partial results of the study have been already published (Yoshida, N. et al. J. Am. Chem. Soc. 2006, 128, 12042-12043). The calculation was carried out for aqueous solutions of three different electrolytes, CaCl2, NaCl, and KCl, and for four different mutants of the human lysozyme: wild type, Q86D, A92D, and Q86D/A92D, which have been studied experimentally. The discussion of this article focuses on the cleft that consists of amino acid residues from Q86 to A92. For the wild type of protein in the aqueous solutions of all the electrolytes studied, there are no distributions observed for the ions inside the cleft. The Q86D mutant shows essentially the same behavior with that of the wild type. The A92D mutant shows strong binding ability to Na+ in the recognition site, which is in accord with the experimental results. There are two isomers of the Q86D/A92D mutant, e.g., apo-Q86D/A92D and holo-Q86D/A92D. Although both isomers exhibit the binding ability to the Na+ and Ca2+ ions, the holo isomer shows much greater affinity compared with the apo isomer. Regarding the selective ion binding of the holo-Q86D/A92D mutant, it shows greater affinity to Ca2+ than to Na+, which is also consistent with the experimental observation.

1. Introduction Molecular recognition is the most fundamental and important function of biomolecules. It is regarded as a process in which a host molecule makes a complex with a guest molecule through noncovalent chemical bonds including hydrogen bonding, hydrophobic interactions, ionic interactions, or other interactions between the host and guest molecules. In a sense, the entire biological cycle of a living system can be viewed as a series of molecular recognitions. One of the most elementary processes of molecular recognition is the selective ion binding by protein. A variety of functions of protein is related to the ion binding: ion channels, ligand binding by a receptor, enzymatic reactions, and so on.1 For example, a calcium ion bound to calmodulin works as a coenzyme for controlling activity of other enzymes.2,3 An ion works sometime as a “clamp” in the construction of protein to stabilize its conformation. Ion binding is quite essential in a variety of ion channels too. The ion binding process by protein is usually so specific that failure of discrimination may cause disease or death to a living system. No wonder there are so many studies devoted to clarify the origin of the selectivity of ion binding by protein.4-6 A theoretical prediction of the ion binding by protein, however, is nontrivial for any type of theories in chemical physics. The reason why it is so is because the problem involves stability or the free energy of ions in extremely inhomogeneous environments of protein immersed in an electrolyte solution, * Corresponding author. Phone: +81(Japan)-564-55-7314. Fax: +81(Japan)-564-53-4660. E-mail: [email protected]. † Institute for Molecular Science. ‡ The Graduate University for Advanced Studies.

which extends over an infinite range or in the thermodynamic limit. The free energy or the chemical potential of an ion not only in protein but also in bulk solution should be accounted properly, because whether or not an ion is recognized by protein is determined by the difference of the chemical potentials in bulk and that in protein. The highest barrier for theories to overcome are processes called “hydration” and “dehydration”, which make crucial contributions to the chemical potentials of an ion. Whether an ion accommodated inside protein is hydrated or dehydrated is determined by a detailed balance of interactions among ions, water molecules, and protein.7 A successful theory of ion binding by protein should be the one that can account for the hydration and dehydration of an ion in bulk solution and in protein without making any ad hoc assumptions. It is needless to say that pure quantum chemistry is helpless for such problems involving the free energy. The problem seems also hopelessly difficult for the molecular simulation. Just imagine how to simulate the ion binding by protein in aqueous electrolyte solutions. It is easy to understand that the straightforward simulation of the ion binding does not converge in a sensible length of time, because it is essentially a matter of stability of ions inside and outside the protein or the problem of free energy difference. A usual maneuver in such a simulation is the “free energy perturbation” technique.8-11 The technique can be applied to the problem, if one knows in advance information such as at which pore the ion is bound, whether the pore was filled with water or not before the ion penetrates into, and/or whether the ion is hydrated in the pore. The information should be supplied from experiments or other sources. However, if such information is already known, there is no need to carry out the simulation. The perspective seems not so bright for the theories based on the continuum solvent

10.1021/jp0685535 CCC: $37.00 © 2007 American Chemical Society Published on Web 04/12/2007

Ion Binding Probed by 3D-RISM model. A conventional way to realize the molecular recognition by means of the continuum model is to evaluate the free energy difference of a guest molecule between the two states, one inside the host molecule and the other in the bulk solution, based on the Poisson-Boltzmann equation or similar theories such as the generalized Born equation. The continuum model requires essentially two types of phenomenological parameters as boundary conditions: one that separates the region of solute from that of solvent, or the size parameter, and the other the dielectric constant of each region.12-14 It is a usual maneuver to determine those parameters empirically comparing the theoretical results for phenomenological observables with experiments. The maneuver seems to work well more or less as far as a solute in bulk solutions is concerned, although the physical meaning of the parameters so determined is not necessarily so transparent. However, it will be disaster if one tries to apply the same maneuver to a solute in a cavity of protein. Are there any experimental data that can be compared with the phenomenological observables at the boundary inside cavity? We suppose there are no such data. Then, how can those boundary conditions, the size of solute, and the dielectric constant be determined by the theory? After all, the difficulty originates deeply from the nature of the theory itself. The continuum theory is founded on a big premise, that is, the number of molecules in the most elementary volume of the system is so large that one can ignore characteristics of a molecule. On the other hand, characteristics of molecules are essential in the molecular recognition problem. Quite recently, a new theoretical methodology has been introduced by Imai et al. to tackle the problem, which is based on the statistical mechanics of molecular liquids, or the threedimensional (3D) RISM theory.15 The method consists of two steps. In the first step, the site-site pair correlation functions of solvent are obtained on the basis of the ordinary (1D) RISM theory. The process corresponds to the preparation of solvent in experiments. Then, the three-dimensional distribution of solvent sites (atoms) around protein is calculated on the basis of the 3D-RISM theory in the succeeding step. The step corresponds to detecting solvent distribution by means of the X-ray or neutron diffraction measurements.16 In the second step, the atomic coordinates of protein are taken from the protein data bank (PDB) by removing all the coordinates concerned with solvent including electrolytes. So, no experimental information with respect to solvent is put in. By carrying out such a calculation for the hen egg-white lysozyme, Imai et al. could have detected four water molecules in a cavity of the protein, where those molecules are supposed to be according to X-ray crystallography. The success lies in the nature of the theory itself, which produces the equilibrium distribution of water in the thermodynamic limit, and thereby, it does not depend on how water molecules are trapped inside protein. In the previous paper, we have extended the method to the electrolyte solutions in order to realize the selective ion binding by protein.17 In particular, we looked at the selective ion binding by different mutants of human lysozyme produced by sitedirected mutagenesis, for which the neutron and X-ray diffraction studies have been carried out by Kuroki et al.18-21 In the present paper, we carry out the additional calculations and detailed analyses for the selective ion binding by human lysozyme and its mutants. 2. Method The distribution of water, cation, and anion around the lysozyme and its mutants in electrolyte solutions are calculated

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4589 by the reference interaction site model (RISM) integral equation theory.15,16 Because the RISM theory has already been described in detail, a brief outline of the theory is provided here. The Ornstein-Zernike integral equation for the mixture is written as22

hRβ(12) ) cRβ(12) +

Fγ ∫ cRγ(13) hγβ(32) d(3) ∑ γ∈All species

(1)

where the superscripts in Greek letters denote the species that composes a mixture and Fγ is the number density of species γ. h(12) and c(12) denote total and direct correlation functions, respectively. The summation in the right-hand side runs over species in a mixture. In this study, we consider a lysozyme molecule, which is immersed in electrolyte solutions. In the infinite dilution limit, the OZ equation can be decoupled into two sets of equations: one for the solute-solvent system and the other for the solvent-solvent system. The OZ equation for a solute-solvent system is written as follows,

hRs(12) ) cRs(12) +

Fs′ ∫ cRs′(13) hs′s(32) d(3) ∑ s′∈Solvent species

(2)

The OZ equation for a solvent-solvent system is

hss′′(12) ) css′′(12) +

Fs′ ∫ css′(13) hs′s′′(32) d(3) ∑ s′∈Solvent species

(3)

It is an essential step in the RISM theory to express the pair correlation functions between two molecules in terms of those between interaction sites, which can be accomplished by averaging the functions over orientations fixing the distance between the interaction sites:

hRs ab(r) )

1 Ω2

∫ d(1) d(2) δ(R1 + la1)δ(R2 + lb2 - r) hRs(12) (4)

where R1 and la1 indicate position of molecule 1 in the laboratory frame and position of site a of molecule 1 in the molecular frame. The RISM equation can be derived from eqs 3 and 4 with a superposition approximation for the direct correlation function as s′′ ∑ $sac / css′′ cd (r) / $db + s′s′′ Fs′ ∑ $sac / css′ ∑ cd (rcd) / hdb (rdb) s′∈Solvent species c,d∈Site on s,s′

hss′′ ab (r) )

c,d∈Site on a,b

(5)

where ω is an intramolecular correlation function and the asterisk denotes the convolution integrals. On the other hand, it is preferable to account for the explicit orientation dependence of the solute-solvent correlation functions, because the protein is an extremely anisotropic molecule. Therefore, we employ the 3D-RISM theory for the solute-solvent system.23 As distinct from the usual (1D) RISM theory, the orientation average is taken over just for solvent molecules, not for the solute molecule. This averaging reduces eq 2 to the 3D-RISM equation:

4590 J. Phys. Chem. B, Vol. 111, No. 17, 2007

hsa(r) )

∑

∑

Yoshida et al.

cs′c (rc) / [$sca + Fs′ hs′s ca (rca)]

s′∈Solvent species c∈Site on s′

TABLE 1: Summary of Lennard-Jones Potential Parameter of Solvent Species and Structure of Solvent Watera

(6) In order to complement these equations, we employ the KH closure for both the solute-solvent and the solvent-solvent systems.24,25 The KH closure for the solute-solvent system is expressed as

gsa(r) )

{

exp(dsa(r)) for dsa(r) e 0 1 + dsa(r) for dsa(r) > 0

(7a)

dsa(r) ) -βusa(r) + hsa(r) - csa(r)

(7b)

H of water O of water ClK+ Na+ Ca2+ O-H [Å] ∠HOH [deg]

σ [Å]

 [kcal/mol]

q [e]

0.400 3.166 4.417 4.935 3.330 2.412 1.000 109.47

0.0460 0.1550 0.1178 0.0003 0.0028 0.4497

0.41 -0.82 -1.00 1.00 1.00 2.00

a

Parameter set of water employed SPC model.34 Parameters of ions were found in OPLS parameter set.35,36

SCHEME 1

where ua(r) and ga(r) ) ha(r) + 1 denote the interaction potential and three-dimensional distribution function (3D-DF) between the solute molecule and the solvent site a at position r, respectively. β ) 1/kBT, kB is the Boltzmann constant, and T is the absolute temperature. For a solvent-solvent system, a similar expression of KH closure is employed. The interaction potential is described as the sum of the electrostatic interaction and Lennard-Jones potential as follows:

usa(r) )

∑

b∈solute

qubqsa |r - rub|

+

∑

{( ) ( ) } σab

4ab

b∈solute

12

|r - rub|

-

σab

6

|r - rub|

(8)

where qa denotes a partial charge on site a and σ and  are the Lennard-Jones parameters with usual meanings. In Table 1, the potential parameters and structure codes using these calculations are summarized. We employ the Lorentz-Verthelot combination rule to combine the Lennard-Jones parameters. The schematic description of the computational procedure is given in Scheme 1. In order to obtain the distribution of solvent species, e.g., cation, anion, and water, around the lysozyme, through eqs 6-8, we begin with the calculation of the sitesite pair correlation functions in the solvent system by the RISM equation, eq 5. Three types of solvent systems, or aqueous solutions of 0.01 mol/kg CaCl2, NaCl, and KCl, are considered at 298 K under usual pressure. In this study, we employ the dielectrically consistent RISM method, which requires an experimental value of the dielectric constant.26 The dielectric constant for pure water, 78.5, is employed for all the solvent systems. The solvent-solvent total correlation functions, hss′′ ab (rab), calculated by RISM are used to perform the solutesolvent 3D-RISM calculation. In the (3D) RISM calculation, the Amber-99 parameter set was employed as potential parameters for the protein.27 Finally, we calculate the threedimensional distribution of solvent (water and electrolytes) around the solute by eqs 6 and 7. It is a straightforward task to obtain the one-dimensional distribution, or the radial distribution function (RDF), from the 3D-distribution function. The RDFs can be obtained by averaging the 3D-distribution function over the direction around a specified center:

g1D a (r,r0) )

1 4π

∫ ga(r0 + r) drˆ

(9)

where rˆ is a direction of vector r and r0 indicates a center for the averaging. In this article, the R-carbon (CR) and carbonylate

a

Details of the calculation and theory of D-RISM are in ref 26.

Figure 1. Atoms are chosen as an averaging center point of holoQ86D/A92D mutant. CR denotes R-carbon of residue 92; OD1 and OD2 are the carbonylate oxygens of aspartic acid of residue 92. In the case of the A92D mutant, CR is chosen by the same way.

oxygen (OD1, OD2) of residue 92 were chosen as the averaging center. (See Figure 1.) 3. Results and Discussion 3.1. Properties of Electrolyte Solvent. As has been stated in the previous section, our treatment requires the solventsolvent total correlation functions as an input to the 3D-RISM calculation for the solvent distribution around the protein. Therefore, we begin the present section with the results of the RISM calculation for the three types of electrolyte solutions: 0.01 mol/kg CaCl2, NaCl, and KCl aqueous solutions at the ambient condition. The purpose of the calculation is not only to prepare the input for the 3D-RISM calculation, but also to examine if the LJ parameters for the ions reasonably describe

Ion Binding Probed by 3D-RISM

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4591

Figure 2. Density dependence of mean activity coefficient of electrolyte solution compared with experimental value and RISM evaluation.

the thermodynamic properties of the electrolyte solutions well. The thermodynamic consistency of the parameters is particularly important in the present case, because the recognition of an ion by a protein is determined by the free energy difference of the ion inside protein and in bulk solutions, which is determined in turn by the interaction parameters, especially by the ion size. Therefore, it is a necessary condition for a theory and a model to be able to evaluate the thermodynamic quantity of ions in bulk solutions. In order to examine the validity of the results of the electrolyte solution, the mean activity coefficient of ion was calculated. The mean activity coefficient, γ(, is defined as28 ν- 1/(ν++ν-) γ( ) (γν+ + γ- )

(10)

Figure 3. The three-dimensional distribution of water oxygen (red), Na+ (yellow), and Cl- (green) around wild-type human lysozyme. Carbonylate oxygen of protein is colored by blue and amino nitrogen by purple. Threshold of 3D-DFs are 3.0 for water oxygen and Na+ and 5.0 for Cl-. Large yellow distribution near the center of figure is located around carbonylate oxygen of Asp-91. All figures with 3DDFs in the present paper were described by the graphical package gOpenMol.33

where ν+ and ν- are the ratio of number density of ions, respectively, and the activity coefficient γi of species is expressed as

RT ln γi ) ∆µi - ∆µ0i

(11)

where ∆µi and ∆µi0 are the excess chemical potentials of species i in electrolyte solution and in pure water, respectively. In Figure 2, the mean activity coefficient is compared with experimental values. The results in general show fair agreement with the experimental results. The theory especially discriminates the divalent ion from the monovalent ions quite well. Apparently, the concentration dependence of the two monovalent ions is not resolved well. However, it will not seriously influence the results for the ion recognition by protein, because the process is determined primarily by the free energy difference of the same ion inside protein and bulk solutions.29 3.2. Wild Type and Q86D Mutant. The 3D-RISM calculations were carried out for the wild type and Q86D mutant in aqueous solutions of 0.01 mol/kg NaCl. The structure of the wild type and the A92D mutant is designated as entries 1LZ1 and 1I1Z in the Brookhaven Protein Data Bank, respectively.18,30 In order to apply the 3D-RISM theory, all heteroatoms (water and ions) were removed from the structures. In Figure 3, 3DDFs of water oxygen, Na+, and Cl- at the surface of protein are shown. Conspicuous distribution of Na+ is seen at the center of the figure, which is bound to the carbonylate oxygen of Asp91. We also observe the oxygen distribution (red spot) near Asp91, which is one of the amino acid residues making the active site of protein. Figure 4 shows an expanded view of the water and electrolyte distributions around residues Q86, D91, and A92, which make

Figure 4. 3D distribution of oxygen of solvent water (red) and Na+ (yellow) around wild type (a) and Q86D mutant (b). Oxygen of crystallized water detected by X-ray experiment is described as gray colored sphere in (b). Threshold of 3D-DFs are 5.0 for all species.

the active site of the protein. Hereafter, we call the active site a “cleft”. In the figure, the distribution function, g(r), greater than five is painted with a color: red for oxygen and yellow for Na+. g(r) > 5 implies that the probability of finding water molecules at the position is 5 times greater than that in the bulk. Taking a glance at the figure, it is obvious that there is no Na+ distribution in the cleft of the wild-type protein, because we see only the red spots, which indicate the distribution of water oxygen. We have also performed the calculation for aqueous solutions of the other electrolytes: CaCl2 and KCl. However, we have not observed the distribution for any of those ions in the cleft. The Q86D mutant exhibits essentially the same behavior with that of the wild type but with the water distribution changed slightly. (There is a trace of yellow spot that indicates a slight possibility of finding a Na+ ion in the middle of the binding site, but it would be too small to make a significant contribution to the distribution.) Instead, the distribution cor-

4592 J. Phys. Chem. B, Vol. 111, No. 17, 2007

Yoshida et al.

Figure 5. RDFs of oxygen of solvent water and Na+ around wild type and Q86D mutant. The R-carbon of residue 92 was chosen as the origin for averaging the 3D-DFs.

responding to the water oxygen is observed as is shown in the red color in Figure 4b. The distribution covers faithfully the region where the crystallographic water molecules have been detected, which are shown with the spheres colored gray. There is a small difference between the theory and the experiment, which is the crystallographic water bound to the backbone of Asp-91. The theory does not reproduce the particular water molecule by unidentified reasons. Except for the difference, the observation is consistent with the experimental finding, especially that the protein with the wild-type sequence binds neither Na+ nor Ca2+. In Figure 5, the 1D distributions or RDFs are shown. The RDF is evaluated by eq 9. In this article, we employ the R-carbon of residue 92 as an averaging center to evaluate the RDFs. The RDF of the water oxygen around the wild type and Q86D have the first minimum at 6 Å, which coincides with the size of the cleft. Although no peaks are found in the RDF of ions inside the cleft, a conspicuous peak is observed at 9 Å. The peak implies the binding of the cation by the carbonyls of Asp-91 and Asp-86 that point outward the cleft. Such a feature also appears in Figure 4. It may be of interest for some purpose to find out how many water molecules are accommodated inside the cleft. The number can be found from the running coordination number of the water oxygen, which is calculated from the RDF of the water oxygen using the equation

CN(r) ) 4πFγ

∫0r gγ(s)s2 ds

(12)

The equation defines the number of water molecules inside the sphere of radius r. Shown in Figure 6 is the running coordination number of the water oxygen calculated from the RDF (Figure 5) of the Q86D mutant using eq 12. The CN(r) obtained from the X-ray structure (1I1Z) is also plotted for comparison. (The X-ray result shows digitized increments, because the plot is reproduced from the gray spots in Figure 4b.) The number of waters in the cleft was slightly underestimated compared with the experimental results. Such an error might have originated from the closure relation, eq 7, we have employed.24,25 In fact, it has been reported by Imai et al. that the HNC closure gives better results for CN(r) in the case of protein in pure water.15 However, in the present calculation, we stay with the KH closure, because the HNC closure shows an extremely poor convergence for the electrolyte solution. There might be another reason for this discrepancy, which is the “frozen” protein structure used in our calculation. The results may change if we relax the protein structure upon binding ions.

Figure 6. Running coordination number of solvent water around R-carbon of residue 92 of Q86D mutant compared with X-ray results. The theoretical values are estimated from the RDFs in Figure 5.

Figure 7. 3D distribution of oxygen of solvent water (red), K+ (purple), and Na+ (yellow) around A92D mutant in KCl (a) and NaCl (b) solution. Threshold of 3D-DFs are 5.0 for all species.

3.3. Ion Binding Mutant. The distributions of water oxygen and cations, Na+ and K+, around the cleft of the A92D mutant are compared in Figure 7. The structure of the A92D mutant is taken from the Brookhaven Protein Data Bank with the entry number, 1I20.18 The A92D mutant shows marked distribution of Na+ in the cleft. The Na+ ion seems to be bound to the negatively charged carbonylate oxygen of Asp-92. The oxygen atoms of water are also distributed in the cleft along with the Na+ ion, indicating that the ion is not entirely naked there. The distribution of the ion is also seen around the carbonyl of Asp91, which is pointing toward the bulk solution. On the contrary to the case of sodium ions, the K+ ion is not found near the cleft. Because Na+ and K+ have the same charge, the selectivity should originate from the difference in size of the ion species. The Lennard-Jones diameters of ions are summarized in Table 1. It is rather obvious that a K+ ion does not penetrate into the cleft because of the size. The Ca2+ distribution around the A92D mutant in CaCl2 aqueous solution was also examined. The Ca2+ ion was found in the cleft contrary to the experimental result. The distributions of Na+ and Ca2+ look quite similar. However, if one examines the RDF of those cations in Figure 8, the first peak for Na+ is slightly higher and broader compared with Ca2+, and Ca2+ has a much higher second peak. Those observations suggest that Ca2+ has less affinity to the cleft in spite of its greater charge. This may be because Ca2+ should pay for greater dehydration penalty upon penetration compared with Na+.31 The same procedure applied to a doubly mutated lysozyme, Q86D/A92D, which has been also investigated experimentally by Kuroki and Yutani. There are two isomers for the Q86D/ A92D mutant, “apo-Q86D/A92D” and “holo-Q86D/A92D”. In biochemistry, the word “holo” generally means the complete and operative form of an enzyme, which consists of protein

Ion Binding Probed by 3D-RISM

Figure 8. RDFs of K+, Na+, and Ca2+ for A92D mutant. The R-carbon of residue 92 was chosen as the origin for averaging the 3D-DFs.

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4593

Figure 10. RDFs of Na+ and Ca2+ in apo- and holo-Q86D/A92D mutant. The R-carbon of residue 92 was chosen as the origin for averaging the 3D-DFs.

Figure 9. 3D distribution of oxygen of solvent water, Na+, and Ca2+ of apo- and holo-Q86D/A92D mutant: (a) apo in NaCl, (b) apo in CaCl2, (c) holo in NaCl, (d) holo in CaCl2. Threshold of 3D-DFs are 5.0 for all species.

and cofactors, and “apo” is an enzyme without its cofactor. Here, we call the isomers with and without calcium binding ability “holo” and “apo,” respectively. The structures of the apo- and holo-Q86D/A92D mutants are designated as entries 2LHM and 3LHM in the Brookhaven Protein Data Bank, respectively.19 In the case of apo-Q86D/A92D, the carbonyl groups of both Asp-86 and Asp-91 are directing to bulk solution, and in the case of holo-Q86D/A92D, carbonyls of Asp-86 and Asp-91 are pointing inward the cleft. In Figure 9, the ion distributions around the cleft are depicted for both apo- and holo-Q86D/ A92D. In both cases, the K+ ion shows a tiny peak at the center of the cleft (not shown). On the other hand, marked distributions of Na+ and Ca2+ ions were found for both isomers. Although the distribution of Na+ and Ca2+ seems to look similar, the height of peaks in the cleft is quite different: 10.6 and 19.1 for Na+ and CA2+ in apo-Q86D/A92D, respectively, and 32.2 and 63.3 for Na+ and Ca2+ in holo-Q86D/A92D, respectively. Because three carbonylate oxygen atoms of the holo-Q86D/ A92D mutant are pointing inward the cleft, the holo mutant can strongly bind both cations, Na+ and Ca2+. According to the experimental results, the Q86D/A92D mutant undergoes isomerization from apo to holo conformers in the presence of

Figure 11. RDFs of oxygen and hydrogen of water and calcium ion around holo-Q86D/A92D mutant. OD1 and OD2 (see Figure 1) were chosen as the averaging center for parts a and b, respectively.

Ca2+ ions. As a result of this isomerization, holo-Q86D/A92D selectively binds Ca2+ in the cleft. In Figure 10, RDFs of Ca2+ and Na+ for apo- and holo-Q86D/A92D are shown. The first peak of Ca2+ is twice as high as that of Na+ for the holo mutant. Therefore, the holo-Q86D/A92D mutant shows selective ion binding of Ca2+, which is in good agreement with experiment. The apo mutant also shows affinity to both Na+ and Ca2+ ions. However, the binding ability of the apo mutant is considerably lower than that of the holo mutant as can be inferred from the height of the first peaks. It may be quite difficult to clarify experimentally whether the apo isomer has the Ca2+ binding ability or not, because the Q86D/A92D mutant isomerizes into the holo mutant in the presence of Ca2+ in the real experimental condition. Our results suggest a mechanism for the ion binding by the Q86D/A92D mutant. The mutant has considerable

4594 J. Phys. Chem. B, Vol. 111, No. 17, 2007

Figure 12. Comparison of Ca2+ position between X-ray result and theoretical estimation. Threshold of 3D-DF is 25.0 for Ca2+ in the righthand side figure.

binding ability to Na+ and Ca2+ already in its apo form. After binding the cations, it would undergo isomerization from the apo to holo conformers in order to get more stability. Further investigation to clarify the mechanism is in progress in our group. Depicted in Figure 11 are RDFs of oxygen, hydrogen, and Ca2+ averaged around the OD1 and OD2 (see Figure 1). Both oxygen atoms in the carbonyl groups of Asp-92 strongly bind water hydrogen. A salient peak of Ca2+ is found only around OD1. On the other hand, the peak of water oxygen around OD2 is greater than that around OD1. Consequently, water is bound to OD2, and Ca2+ is bound to OD1, located at center of the cleft. In Figure 12, the position of the Ca2+ binding site is compared with the result of X-ray crystallography.19 In the figure, the distribution function, g(r), greater than 25 is painted. The highest peak of Ca2+ around the holo-Q86D/A92D mutant is located at almost the same position as the experimental results. The results indicate that the Ca2+ binding site empirically observed is the one that has the highest probability of finding the ion in the theoretical description. 4. Conclusion In this article, we have presented theoretical results for the ion binding by human lysozyme on the basis of the 3D-RISM theory. The ion distribution around the wild type, Q86D, A92D, and Q86D/A92D mutants in several electrolyte solutions, KCl, NaCl, and CaCl2, has been evaluated. The doubly substituted mutant, Q86D/A92D, has two isomers distinguished with whether it has a Ca2+ ion or not: apo, without Ca2+; holo, with Ca2+. Because the difference between wild type and mutants is only the active site, the discussion has been focused on the active site, which consists of amino acid residues from Q83 to A92. The wild type and the Q86D mutant show no cation binding ability in accord with the experimental results. The 3D-RISM theory indicates that the A92D and Q86D/A92D mutants have cation binding ability. Na+ and Ca2+ ions are bound by the active site of the A92D and Q86D/A92D mutants, though K+ ions are not found in the active site. Although the experimental result indicates that A92D binds Na+ only, the theoretical prediction shows that the mutant has the ability to recognize both Na+ and Ca2+. Both isomers of the Q86D/A92D mutant, apo and holo, exhibit strong affinity to Na+ and Ca2+. In the case of holo-Q86D/A92D, Ca2+ shows a peak of distribution that is higher than Na+. It implies that the holo-Q86D/A92D binds Ca2+ preferably in accord with the experimental results.19 In the present study, we could have demonstrated the extraordinary ability of the 3D-RISM method to “detect” ions that are bound selectively by different mutants of a protein. However, the achievement is just a half of the way to our

Yoshida et al. ultimate goal, because our prediction is a priori dependent on the protein structure, which has been obtained from experiments. Our ultimate goal is to predict the ion recognition by protein without depending on the conformation determined by experiments. It can be done in principle by performing the conformational sampling of protein in electrolyte solutions, or the folding simulation, on the free energy surface produced by the 3D-RISM theory.16,32 It is, however, rather a long-term goal. On the way to the ultimate goal, we can set an intermediate target for predicting the selective ion binding and the conformational change of protein simultaneously. Let us take the human lysozyme as an example. We know the wild type of the protein does not bind cations in its active site. Let us substitute an amino acid with other ones to create a mutant, say A92D. The conformation of the protein is not only unstable due to the distortion introduced by the mutation, but also in nonequilibrium state with the aqueous solution of sodium chloride. Then, we perform a conformational search on the free energy surface in the limited space around the mutant structure in order to find the stable conformation in the electrolyte solution. If the protein structure obtained is close enough to the experimental one and if we find the distribution of cation in the active site, it means that our intermediate goal is attained. We can apply the method to an unknown mutant to predict its ability of ion binding. Acknowledgment. We are grateful to Dr. Irisa for bringing the problem of calcium binding by lysozyme. Authors are also grateful to Dr. Maruyama for his computational support and to Dr. Imai for his helpful discussion. This work is supported by the Grant-in-Aid for Scientific Research on the priority area of “Water and Biomolecules” from the Ministry of Education in Japan, Culture, Sports, Science, and Technology (MONBUKAGAKUSHO). We are also grateful to the support by the grant from the Next Generation Supercomputing Project, Nanoscience Program of the ministry. References and Notes (1) Hille, B. Ionic Channels of Excitable Membranes; Sinauer Associates: Sunderland, MA, 2001. (2) Bachs, O.; Agell, N. Calcisum and calmodulin function in the cell nucleus; R. G. Lande: Austin, TX, 1995. (3) Calcium and Cell Function; Cheung, W.-Y, Ed.; Academic Press: New York, 1982. (4) Aqvist, J.; Luzhkov, V. B. Nature. 2000, 404, 881-884. (5) Luzhkov, V. B.; Aqvist, J. Biochim. Biophys. Acta. 2005, 1747, 109-120. (6) Luzhkov, V. B.; Almlof, M.; Nervall, M.; Aqvist, J. Biochemistry 2006, 45, 10807-10814. (7) Chong, S.-H.; Hirata, F. J. Phys. Chem. B 1997, 101, 3209-3220. (8) Torrie, G. M.; Valleau, J. P. Chem. Phys. Lett. 1974, 28, 578581. (9) Computer Simulation of Biomolecular Systems; van Gunsteren, W. F., Weiner, P. K, Eds.; ESCOM: Leiden, Holland, 1988. (10) Karplus, M.; Petsko, G. A. Nature 1990, 347, 631-639. (11) Kollman, P. Chem. ReV. 1993, 93, 2395-2417. (12) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027-2094. (13) Cramer, C. J.; Truhlar, D. G. Chem. ReV. 1994, 99, 2161-2200. (14) Tomasi, J.; Menncci, B.; Cammi, R. Chem. ReV. 2005, 105, 29993093. (15) Imai, T.; Hiraoka, R.; Kovalenko, A.; Hirata, F. J. Am. Chem. Soc. Commun. 2005, 127, 15334-15335. (16) Hirata, F. Molecular Theory of SolVation; Kluwer: Dordrecht, Holland, 2003. (17) Yoshida, N.; Phongphanphanee, S.; Maruyama, Y.; Imai, T.; Hirata, F. J. Am. Chem. Soc. 2006, 128, 12042-12043. (18) Kuroki, R.; Yutani, K. J. Biol. Chem. 1998, 273, 34310-34315. (19) Inaka, K.; Kuroki, R.; Kikuchi, M.; Matsushima, M. J. Biol. Chem. 1991, 266, 20666-20671. (20) Kuroki, R.; Kawakita, S.; Nakamura, H.; Yutani, K. Proc. Natl. Acad. Sci. 1992, 89, 6803-6807. (21) Kuroki, R.; Taniyama, Y.; Seko, C.; Nakamura, H.; Kikuchi, M.; Ikehara, M. Proc. Natl. Acad. Sci. 1989, 86, 6903-6907. (22) Gray, C. G.; Gubbins, K. E. Theory of molecular fluids; Clarendon Press: Oxford, 1984.

Ion Binding Probed by 3D-RISM (23) Kovalenko, A.; Hirata, F. J. Chem. Phys. 1999, 110, 10095-10112. (24) Kovalenko, A.; Hirata, F. J. Phys. Chem. B 1999, 103, 7942-7957. (25) Kovalenko, A.; Hirata, F. J. Chem. Phys. 2000, 112, 10391-10402. (26) Perkyns, J.; Pettitt, B. M. J. Chem. Phys. 1992, 97, 7656-7666. (27) Wang, J.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049-1074. (28) Fried, V.; Blukis, U.; Hameka, H. F. Physical Chemistry; Machillan, New York, 1977; pp 870-874. (29) Moore, W. J. Physical Chemistry; Prentice Hall: Englewood Cliffs, NJ, 1972; pp 444. (30) Artymiuk, P. J.; Blake, C. C. F. J. Mol. Biol. 1981, 152, 737-762.

J. Phys. Chem. B, Vol. 111, No. 17, 2007 4595 (31) Ikuta, Y.; Matsugami, M.; Maruyama, Y.; Hirata, F. Chem. Phys. Lett. 2007, 433, 403-408. (32) Mitsutake, A.; Hinoshita, M.; Okamoto, Y.; Hirata, Y. J. Phys. Chem. B 2004, 108, 19002-19012. (33) (a) Laaksonen, L. J. Mol. Graphics 1992, 10, 33. (b) Bergman, D. L.; Laaksonen, L.; Laaksonen, A. J. Mol. Graphics Model. 1997, 15, 301. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunstern, E. F.; Hermans, J. Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981. (35) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 903-910. (36) Aqvist, J. J. Phys. Chem. 1990, 94, 8021-8024.