1Institute for Protein Research, Osaka University, Yamadaoka, Suita, Osaka 565-0871 and 2Institute of Molecular and Cellular Biosciences, University of Tokyo, Bunkyo-ku, Tokyo 113-0032, Japan 4Present address: Riken Harima Institute, HTPF Kout, Mikazuki-chou, Sayon-gun, Hyogo 679-5148, Japan
3 To whom correspondence should be addressed. e-mail: kitao{at}iam.u-tokyo.ac.jp
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Keywords: component analysis/free energy perturbation/human lysozyme/modeling denatured-state structure/protein stability
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Recently, a more general relationship between the native structure and the stability of proteins has been found by exhaustive mutational studies on human lysozyme (Funahashi et al., 1999, 2000). In this empirical relationship, the Gibbs free energy change is expressed as a sum of several free energy components as functions of the ASA of polar and non-polar residues, the distance of hydrogen bonds and the number of introduced water molecules due to the substitution. Although this relationship has been derived by the least-squares fitting of a large number of experimental data in mutants of human lysozymes, it was also successful in explaining the stability changes caused by the amino acid substitutions of other proteins (Funahashi et al., 2001
). Hence the empirical relationship is considered to be very useful for understanding general tendencies in protein stability. However, there are a few mutants whose stability change cannot be explained by the relationship. This suggests that more detailed analysis including other important stabilization factors, such as the denatured structure, may be necessary in such cases.
Here, we used the free energy calculation based on the all-atom molecular dynamics (MD) simulation (Kollman, 1993) in order to study the effect of the amino acid substitutions on conformational stability in detail. In the free energy calculation method, the free energy difference,
G, is calculated as a difference
G =
GN
GD, where
GN and
GD are defined as the free energy change from wild-type to mutant along the alchemical path in the native state and in the denatured state, respectively. Therefore, it is necessary to carry out the simulations not only in the native state but also in the denatured state in order to obtain
G. Since there is no structure experimentally determined in atomic detail in the denatured state, we have proposed a few models for the denatured state. Because calculation of
GN is highly reliable, we can select the possible model structure in the denatured state by comparing
G determined by different denatured-state models with the experimental
G (
Gexp). After obtaining the possible structure of the denatured state, we investigate the mechanism of the stabilization or destabilization of the mutant proteins by detailed atomic analysis, such as free energy component analysis as discussed later. This is the computational strategy that we proposed in previous papers in order to understand both the mechanism of the protein stability and the structure of the denatured state (Sugita and Kitao, 1998a
). In the free energy calculation, the obtained
G can be decomposed into several components (Shi et al., 1993
; Smith and van Gunsteren, 1994
). This analysis, which is often referred to as free energy component analysis, gives us useful information, which cannot be obtained directly by experimental studies. The free energy components are intrinsically dependent on the pathways as to how the structure and parameters change from the wild-type to the mutant proteins (Shi et al., 1993
; Smith and van Gunsteren, 1994
). However, it has been shown that this analysis is meaningful when the path is carefully chosen and clearly defined (Boresch and Karplus, 1995
). Since the free energy component is more sensitive than the total free energy changes, we employ the accurate treatment of the electrostatic interactions and the boundary condition of a protein and surrounding water molecules. We have already shown that the software introducing these advanced treatments can give us more reliable results than conventional methods (Sugita and Kitao, 1998b
).
The thermal stability of human lysozyme has been exhaustively studied experimentally using calorimetry and X-ray crystallography (Takano et al., 1995, 1997, 1998, 1999a,b, 2001; Funahashi et al., 1996, 1999, 2000, 2001, 2002
; Yamagata et al., 1998
). In theoretical studies, the relative stability of four mutant human lysozymes, I23V, I56V, I89V and I106V, has already been reported (Sugita et al., 1998
). The study was concentrated on how the relative stability is changed by the same type of mutations (I to V) at the different mutation sites. In the present study, we focused on the stability changes due to the different types of mutations at the same site. We carried out free energy calculations on two mutations (I56A and I56F) and compared the stability changes caused by these mutants including I56V reported previously. This mutation site is one of the key residues for elucidating the mechanism of amyloidosis (Pepys et al., 1993
). In addition, thermodynamic and structural data for several hydrophobic mutants on the 56th residue (see Figure 1) are available and have already been examined in the framework of the empirical relationship between structure and stability (Funahashi et al., 1999
). In this paper, we compare the results of free energy calculations studied here with those obtained by an empirical relationship to elucidate the detailed mechanism of protein stability caused by different types of mutations at the same position of human lysozyme. We will show that the effect of the long-range interaction on the stability changes may be underestimated in empirical relationships when the structural change caused by mutation is relatively small. It will also be suggested that estimation of the change in accessible surface area,
ASA, may be overestimated if structure around the mutation site in the denatured state is native-like, which would cause the overestimation of the free energy change.
|
![]() |
Materials and methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
The mutational free energy changes in the native, GN and denatured state,
GD, are calculated by the thermodynamic integration (TI) equation:
where H() is an intermediate Hamiltonian and
is a coupling constant between the wild-type (
= 0) and mutant (
= 1). The bracket < >
denotes a canonical ensemble average for the system characterized by the Hamiltonian H(
). Free energy change calculated by the TI equation can be divided into energy components such as
G =
Gcov +
GLJ +
Gel =
Gcov +
Gnonbd(2)
where Gcov,
GLJ,
Gel and
Gnonbd are the covalent, LennardJones, electrostatic and non-bonded free energy components, respectively. Furthermore, we divide the non-bonded free energy changes into the residual free energy component,
Gres, as
G =
Gres(3)
In the current study, we employed the computer program for free energy calculation developed previously (Sugita and Kitao, 1998a,b; Sugita et al., 1998
). High-precision algorithms such as the cell multipole method (CMM) (Ding et al., 1992
) for electrostatic interaction and the spherical solvent boundary potential (SSBP) (Beglov and Roux, 1994
) for boundary conditions of the system including a protein and surrounding water molecules are incorporated in this program. Compared with the conventional approximation, not only the hysteresis of the free energy changes but also the errors of free energy components become greatly reduced (Sugita and Kitao, 1998b
).
Simulation procedures
The AMBER all-atom energy function (Weiner et al., 1986) is used for the protein and peptide models and the TIP3P model (Jorgensen et al., 1983
) is employed for water molecules. Although the AMBER force field developed by Weiner et al. is not the latest version of the AMBER force field, we use it to compare the current results with those in our previous studies on the same protein. The CMM is employed to evaluate the electrostatic interactions within the sphere and SSBP is applied to include the effect of reaction field and to reduce the artificial effect near the boundary of the sphere. The radii of the spheres including the protein and surrounding water molecules are initially 34 Å in the native state and 20 Å in the denatured state. The system is simulated in an isothermalisobaric ensemble (300 K, 1 atm), achieved by the NoseHoover algorithm (Nose, 1984
; Hoover, 1985
) and SSBP. Water molecules are treated as a rigid body. The time step of MD,
t, is 0.5 fs. Since thermodynamic parameters were measured in low salt conditions (0.05 M Gly.HCl), we did not add large numbers of counter ions around the protein. Only sodium ion, the atomic position of which was determined experimentally, was considered in the simulation.
To equilibrate the system in the isothermalisobaric ensemble, 100 ps MD simulations were performed before the free energy calculation. Twenty intermediate states were prepared by changing the coupling parameters, , as i = 0.025, 0.075, ..., 0.975. The free energy changes,
GN and
GD, are calculated by thermodynamic integration (TI) and double-wide sampling (Jorgensen and Ravimohan, 1985
). In each intermediate state, 5.0 ps MD simulation was performed to equilibrate the system, then 5.0 ps MD simulation was carried out to collect the data for free energy calculations. In total, 200 ps calculations were performed in both the forward and reverse directions. Since I56V mutation had been calculated in our previous study (Sugita and Kitao, 1998a
), we calculated I56A and I56F mutations in this study and compared the effects of three mutations at the same mutation site. In the case of I56A mutation, we started from wild-type human lysozyme and changed the parameters and topology of the 56th Ile residue of human lysozyme into Ala in the forward direction. To obtain the free energy differences caused by I56F mutation, we used I56F mutant as an initial structure and changed the parameters and topology of the 56th residue into alanine (F56A mutation). Then, by considering the thermodynamic cycle, we obtained the free energy differences for I56F from the results obtained for I56A and F56A mutations. In free energy perturbation, we employed the single topology method. During the process of Ile to Ala mutation, two C
atoms change to Hß atoms. H
, H
and C
atoms change to dummy atoms, which do not interact with other atoms. C
changes to Hß and other atoms in the aromatic ring to dummy atoms in I56F mutation. Free energy changes were calculated in reverse directions (Ala to Ile and Ala to Phe mutations) to examine the hysteresis error.
All the MD/FEP/TI calculations were performed by using the program package (Sugita and Kitao, 1998b) developed based on the framework of the Minimization/MD program PRESTO version 2.0 (Morikami et al., 1992
).
![]() |
Results and discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
The free energy calculation method based on MD simulations has now become a common tool in understanding protein stability and functions (Kollman, 1993). In the study of protein stability (Sugita and Kitao, 1998a,b
; Sugita et al., 1998
; Kono et al., 2000
; Pan and Daggett, 2001
), the free energy difference between the wild-type and mutant proteins (
Gcal) is defined as the difference between the mutational free energy changes in the native state,
GN, and in the denatured state,
GD.
Gcal is equal to the experimental free energy differences (
Gexp). It is shown by using the following thermodynamic relations:
G =
GW
GM =
Gexp =
GN
GD =
Gcal(4)
where GW and
GM are the free energy changes between native state and denatured state in the wild-type and the mutant proteins determined experimentally, respectively.
To calculate GD, a five-residue peptide including the mutation site in the middle is employed. Two representative structures for the denatured state are prepared for the calculations in order to examine the structural dependence of the free energy differences on the denatured state structure. In the first model, which we refer to as a native-like model (NTV model), the initial structure was taken from the native structure of the corresponding region. To sample conformational space in the vicinity of the initial conformation, we added a mild constraint for the backbone atoms during the simulations. The second model is an extended model (EXT model) in which the main-chain dihedral angles are set to be 180°. A more detailed discussion on these models can be found in our earlier paper (Sugita and Kitao, 1998a
).
In Table I, the calculated free energy changes in the native and denatured states calculated by two models are listed. When the EXT model is employed, Gcal becomes larger than
Gexp in the cases of I56V and I56F mutations. On the other hand,
Gcal of the I56A mutation is much smaller than the experimental value when the EXT model is employed. The calculated free energy differences are in good agreement with the experimental free energy differences in three mutations when the NTV model is employed. This result suggests that the possible structure near the 56th residue is native-like even in the thermally denatured state (Sugita and Kitao, 1998
). Although we also tried randomly generated initial structures in the case of I56V in our previous study, the native-like structure has been the only model that can reproduce experimental value. In the previous studies, possible structures in three other mutation sites (I23, I89 and I106) were also estimated to be native-like in the same denaturation conditions. Judging from the present and previous results of free energy calculations, the thermally denatured state of human lysozyme is predicted from the native-like structure around the 56th residue.
|
Relation between structural and stability changes
The structures of wild-type and mutants were compared in order to examine the relation between the structural and stability changes. The root mean square deviations from wild-type to I56V, I56A and I56F in the crystal are 0.1, 0.1 and 0.4 Å, respectively (Funahashi et al., 1999). Corresponding values for the final structures in MD simulations are 1.2, 1.1 and 1.2 Å, respectively. They are all in the range of the atomic fluctuations at room temperature. Although
Gexp for I56F is much greater than the others, a significantly large difference was not observed in the native structure of I56F. This is the case where the difference in protein stability caused by the amino acid substitutions cannot be explained only by the structural change in the native state. In the next section, we examine how free energy component analysis can explain the differences in stability caused by the amino acid mutations at the 56th residue.
Comparison of the relative stability of three mutations
In Table II, free energy components, covalent, LennardJones and electrostatic components, are shown. In the cases of I56V and I56A mutations, the major contributions are those originating from LennardJones interactions. They are negative values in both GN and
GD, which means free energies of mutants are lower than wild-type. Since
GD is lower than
GN,
G are positive, i.e. the mutants are less stable than the wild-type protein. In the case of the I56F mutation, the major contribution is that originating from the electrostatic interaction. In this mutation, the LennardJones component is negligible in
G, because those in
GN and
GD are cancelled out. Although Ile and Phe are both hydrophobic residues, they have very different characters, i.e. Ile is an aliphatic residue whereas Phe is an aromatic residue whose partial charges of side-chain atoms are much larger. In this calculation, the partial charges of the H
and H
atoms in the Ile residue are 0.0270.029e, whereas that of hydrogen atoms in the aromatic ring of Phe is 0.15e. Partial charges of carbon atoms are also very different between these two residues. The large difference in partial charges between Ile and Phe is seen not only in the employed force field but also in other force fields. The electrostatic interactions between the mutation site and the other residues are significantly different between I56F mutant and wild-type protein. Large changes in the charge distribution in the mutation site may be the main cause of the electrostatic contribution
|
|
|
|
In this section, we examine the free energy components obtained by the empirical general relationship between structure and stability determined using experimental values (Funahashi et al., 1999, 2000) by comparing them with the free energy components obtained by MD simulation. As pointed out in the Introduction, this method is very useful for understanding the general tendencies of the thermal stability of various proteins although it has been derived from experimental data for human lysozyme.
This empirical relationship was described in detail in the original paper; however, we summarize it here to make the discussion of this component analysis clearer. The changes in the hydrophobic effect, conformational entropy, hydrogen bonding and water molecule embedded in the interior of a protein are often considered to be important stabilization factors. In this treatment, the difference in Gibbs free energy between the wild-type and mutant proteins (Gemp) is expressed by the following equation:
Gemp =
GHP +
Gconf +
GHB +
GH2O(5)
where GHP,
Gconf,
GHB and
GH2O represent the changes in
Gemp due to the hydrophobic effect, the conformational entropy of the substituted residue, the formation and removal of hydrogen bonding and the introduction of water molecules, respectively. In the three mutations studied, the number of hydrogen bonds does not change and no introduced water molecules are found around the mutation sites. Therefore, only the first two terms in Equation 5 are meaningful.
GHP is further decomposed as follows:
GHP =
ASAnon-polar + ß
ASApolar(6)
where ASAnon-polar and
ASApolar represent the difference in the total ASA of all non-polar atoms in the protein between the wild-type and mutant proteins upon denaturation and that of polar atoms, respectively. The coefficients
and ß are fitting parameters, which were determined in a previous study (Funahashi et al., 1999
). It is assumed that the structure around the mutation site is the extended structure.
The results are shown in Table IV. A significant difference between Gexp and
Gemp is seen in the cases of I56V and I56F. Here we further discuss the possible reasons for the difference. One possible reason is that the conformation near the mutation site is native-like even in the denatured state as suggested in the free energy calculation, whereas it is assumed to be extended in the empirical relationship. This may cause the overestimation of
Gconf,
ASAnon-polar and ß
ASApolar. In this case, these values should be smaller because the conformations in the native and denatured states are similar. In I56V, overestimation of
ASAnon-polar may be the main cause of the error because
Gconf is relatively small.
|
Conclusions
The effects of the hydrophobic mutations at the 56th residue in human lysozyme have been studied by free energy calculations. To study the dependence of the calculated free energy difference on the denatured-state structure, we employed two different models for the denatured state, the native-like model and the extended model. From the free energy analysis, it is shown that free energy difference G calculated by the native-like model is in good agreement with the experimentally determined
G. The results strongly suggest that the structure in the vicinity of the mutation site is native-like even in the thermally denatured state, which is also consistent with the previous studies (Sugita and Kitao, 1998a
; Sugita et al., 1998
). This indicates that the model of the denatured state is very important for carrying out free energy analysis. By using free energy component analysis, the changes in thermal stability caused by three mutations were studied. The results indicate that two of the mutants (I56A and I56V) are destabilized by the change in the LennardJones interaction, whereas I56F is largely destabilized by the change in the electrostatic interactions.
We also examined the free energy components determined in the empirical relationship obtained with a large number of structural and thermodynamic data in human lysozyme. A significant difference between Gexp and
Gemp is seen in the cases of I56V and I56F. In the case of I56V, a possible reason for the difference may be the overestimation of
ASAnon-polar because the conformation in the denatured state is assumed to be extended in the empirical relationship. It is also suggested that the effect of the electrostatic interaction on the stability change in I56F mutation is underestimated by the empirical relationship.
The empirical relationship is a powerful method for reproducing and understanding the free energy difference caused by single amino acid mutations in proteins. In this method, the free energy difference can be obtained by relatively simple calculation. Since parameters are determined by the statistical treatment of experimental data, there should be some cases where deviation from the average is significantly large. Since such cases should be limited in number, the free energy simulation based on the all-atom MD simulation can be very useful for understanding them, although it is computationally much more demanding. We conclude that the combined use of free energy calculation based on MD simulation and the empirical relationship is a favorable and promising approach for understanding the mechanism of protein stability in detail.
![]() |
Acknowledgements |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Boresch,S. and Karplus,M. (1995) J. Mol. Biol., 254, 801807.[CrossRef][ISI][Medline]
Buckle,A.M., Henrick,K. and Fersht,A.R. (1993) J. Mol. Biol., 234, 847860.[CrossRef][ISI][Medline]
Ding,H.-Q., Karasawa,N. and Goddard,I.W.A. (1992) J. Chem. Phys., 97, 43094315.[CrossRef][ISI]
Eriksson,A.E., Baase,W.A., Zhang,X.J., Heinz,D.W., Blaber,M., Baldwin,E.P. and Matthews,B.W. (1992) Science, 255, 178183.[ISI][Medline]
Funahashi,J., Takano,K., Ogasahara,K., Yamagata,Y. and Yutani,K. (1996) J. Biochem. (Tokyo), 120, 12161223.[Abstract]
Funahashi,J., Takano,K., Yamagata,Y. and Yutani,K. (1999) Protein Eng., 12, 841850.
Funahashi,J., Takano,K., Yamagata,Y. and Yutani,K. (2000) Biochemistry, 39, 1444814456.[CrossRef][ISI][Medline]
Funahashi,J., Takano,K. and Yutani,K. (2001) Protein Eng., 14, 127134.
Funahashi,J., Takano,K., Yamagata,Y. and Yutani,K. (2002) J. Biol. Chem., 277, 2179221800.
Hoover,W.G. (1985) Phys. Rev. A, 31, 16951697.[CrossRef][ISI][Medline]
Jorgensen,W.L., Chadrasekhar,J. and Madura,J.D. (1983) J. Chem. Phys., 79, 926935.[CrossRef][ISI]
Jorgensen,W.L. and Ravimohan,C. (1985) J. Chem. Phys., 83, 30503054.[CrossRef][ISI]
Klein-Seetharaman,J., Oikawa,M., Grimshaw,S.B., Wirmer,J., Duchardt,E., Ueda,T., Imoto,T., Smith,L.J., Dobson,C.M. and Schwalbe,H. (2002) Science, 295, 17191722.
Kollman,P.A. (1993) Chem. Rev., 93, 23952417.[ISI]
Kono,H., Saito,M. and Sarai,A. (2000) Proteins, 38, 197209.[CrossRef][ISI][Medline]
Koradi,R., Billeter,M. and Wuthrich,K. (1996) J. Mol. Graph., 14, 5155.[CrossRef][ISI][Medline]
Matthews,B.W. (1993) Annu. Rev. Biochem., 62, 139160.[CrossRef][ISI][Medline]
Morikami,K., Nakai,T., Kidera,A., Saito,M. and Nakamura,H. (1992) Comput. Chem., 16, 243248.[CrossRef][ISI]
Nose,S. (1984) Mol. Phys., 52, 255268.[ISI]
Pace,C.N. (1992) J. Mol. Biol., 226, 2935.[ISI][Medline]
Pan,Y. and Daggett,V. (2001) Biochemistry, 40, 27232731.[CrossRef][ISI][Medline]
Pepys,M.B. et al. (1993) Nature, 362, 553557.[CrossRef][ISI][Medline]
Schwalbe,H., Fiebig,K.M., Buck,M., Jones,J.A., Grimshaw,S.B., Spencer,A., Glaser,S.J., Smith,L.J. and Dobson,C.M. (1997) Biochemistry, 36, 89778991.[CrossRef][ISI][Medline]
Serrano,L., Kellis,J.T.,Jr, Cann,P., Matouschek,A. and Fersht,A.R. (1992) J. Mol. Biol., 224, 783804.[ISI][Medline]
Shi,Y.Y., Mark,A.E., Wang,C.X., Huang,F., Berendsen,H.J. and van Gunsteren,W.F. (1993) Protein Eng., 6, 289295.[Abstract]
Smith,P.E. and van Gunsteren,W.F. (1994) J. Phys. Chem., 98, 23662379.
Sugita,Y. and Kitao,A. (1998a) Biophys. J., 75, 21782187.
Sugita,Y. and Kitao,A. (1998b) Proteins, 30, 388400.[CrossRef][Medline]
Sugita,Y., Kitao,A. and Go,N. (1998) Fold. Des., 3, 173181.[ISI][Medline]
Takano,K., Ogasahara,K., Kaneda,H., Yamagata,Y., Fujii,S., Kanaya,E., Kikuchi,M., Oobatake,M. and Yutani,K. (1995) J. Mol. Biol., 254, 6276.[CrossRef][ISI][Medline]
Takano,K., Funahashi,J., Yamagata,Y., Fujii,S. and Yutani,K. (1997) J. Mol. Biol., 274, 132142.[CrossRef][ISI][Medline]
Takano,K., Yamagata,Y. and Yutani,K. (1998) J. Mol. Biol., 280, 749761.[CrossRef][ISI][Medline]
Takano,K., Yamagata,Y., Funahashi,J., Hioki,Y., Kuramitsu,S. and Yutani,K. (1999a) Biochemistry, 38, 1269812708.[CrossRef][ISI][Medline]
Takano,K., Yamagata,Y., Kubota,M., Funahashi,J., Fujii,S. and Yutani,K. (1999b) Biochemistry, 38, 66236629.[CrossRef][ISI][Medline]
Takano,K., Funahashi,J. and Yutani,K. (2001) Eur. J. Biochem., 268, 155159.
Weiner,S.J., Kollman,P.A., Nguyen,D.T. and Case,D.A. (1986) J. Comput. Chem., 7, 230252.[ISI]
Yamagata,Y., Kubota,M., Sumikawa,Y., Funahashi,J., Takano,K., Fujii,S. and Yutani,K. (1998) Biochemistry, 37, 93559362.[CrossRef][ISI][Medline]
Received April 24, 2003; revised July 7, 2003; accepted July 18, 2003.