Journal of Theoretical and Computational Chemistry, Vol. 1, No. 1 (2002) 53–67 c World Scientific Publishing Company
CALCULATING ACCURATE REDOX POTENTIALS IN ENZYMES WITH A COMBINED QM/MM FREE ENERGY PERTURBATION APPROACH MARK S. FORMANECK, GUOHUI LI, XIAODONG ZHANG and QIANG CUI∗ Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin, Madison 1101 University Ave, Madison, WI 53706, USA
[email protected] Received 1 April 2002 Accepted 16 April 2002 An approach for computing accurate redox potentials in enzymes is developed based on the free energy perturbation technique in a QM/MM framework. With an appropriate choice of the QM level and QM/MM coupling scheme, the intermolecular interaction between the redox center and the protein environment can be adequately described; the speed of QM/MM methods also allows a sufficient configurational sampling for the convergence of free energy derivatives. Following the implementation into the simulation package CHARMM, the method was tested with an application to the first reduction potential of FAD in cholesterol oxidase (Chox). In addition to an accurate QM level and adequate conformational samplings, the effect of long-range electrostatic interactions due to the bulk solvent was also found to be essential. Using a semi-empirical density functional theory (SCC-DFTB) as the QM level, and a multi-stage charge-scaling scheme based on Poisson–Boltzmann calculations for the solvation effect, satisfactory agreements with experimental measurements were obtained. The study of Chox also indicates that large errors in the calculated redox potential might arise if changes in the conformational properties of the protein during the redox process are not taken into account, such as in energy minimization type of studies based on only the X-ray structure of the enzyme in one redox state. Keywords: QM/MM; redox potential; solvation; free energy perturbation; cholesterol oxidase.
redox potentials, which are subtly modulated by their surrounding environment. Alteration in the redox potentials for these species often results in the malfunction or damage of biomolecules and thus leads to serious diseases. For example, misregulation of the redox potential of copper ion can lead to radical species that damage lipids, proteins or DNA, and has been implicated in a number of neurodegenerative disorders such as the Wilsons and Menkes diseases.3 Therefore, predicting redox properties for species in biomolecules and understanding how the protein environment regulate these redox quantities are of great importance. Although related quantities such as ionization potential or electron affinity can be predicted
1. Introduction Oxidation-reduction reactions play a key role in many chemical and biological processes.1 Many metabolic processes involve electron transfer and thus redox reactions. As a matter of fact, living systems derive most of their free energy from redox reactions. In photosynthesis, for example, CO2 is reduced and H2 O is oxidized to yield carbohydrates and O2 . Important redox reagents in biological systems include oxygen, di-sulfide units, transition metal ions and a number of co-enzymes such as flavin adenine dinucleotide (FAD) and nicotinamide adenine dinucleotide (NAD). The rate and outcome of the reactions associated with these redox active species depend critically on their ∗Corresponding
author.
53
54 M. S. Formaneck et al.
fairly accurately with ab initio or density functional methods for small and medium size molecules in the gas phase,4 similar quantities in the condensed phase (e.g., in solution or biomolecules) are more difficult to obtain with theoretical calculations. This is because the effects due to the environment have to be described accurately also, which require both a satisfactory description of the intermolecular interactions between the redox center and the environment, and a sufficient sampling of the configurations associated with the environmental degrees of freedom. In most previous studies, either one of these two essential aspects is addressed. For example, in the study of redox properties of iron-sulfur proteins (e.g., ferrodoxin),6 the configurations of the protein were sampled adequately with molecular dynamics, but the iron-sulfur cluster in different redox states was treated with relatively simple partial charge models. Although such calculations can provide useful insights into the effects of the protein environment on redox properties, it is difficult to obtain quantitative results in general, especially for absolute redox potentials. By contrast, density functional theories have been used to describe the electronic structural change during the redox process in systems such as iron-sulfur proteins8 and Cu, Zn superoxide dismutase;9 the interaction with the protein and solution was approximated with a continuum model at the Poisson–Boltzmann level.10 Although the results from these studies are also instructive and may agree very well with experimental results, it is difficult to judge the accuracy of such calculations for general cases. In the current work, we adopted a microscopic approach that addresses both the intermolecular interaction and configurational sampling issues. In particular, we used a combined QM/MM treatment for the system such that the electronic structural change during the redox process can be described accurately. At the same time, the efficiency of the method also allows for a sufficient configurational sampling of the environmental degrees of freedom. As a result, quantitatively reliable values can, in principle (depending on the QM level and treatment of the QM/MM interaction), be obtained for both absolute and relative redox potentials in the condensed phase. In Sec. 2, we briefly describe the theory behind the QM/MM free energy perturbation approach and the implementation in CHARMM. Issues critical
to the accuracy of the results, such as the treatment of long-range electrostatic interactions, will also be discussed. In Sec. 3, we report preliminary results for a study of the first reduction potential of flavin adenine dinucleotide (FAD) in cholesterol oxidase, as an illustration of the QM/MM free energy perturbation approach and the importance of proper treatments of solvation effects. A few conclusions are summarized in Sec. 4. 2. Theory and Implementation To compute redox potentials in the condensed phase, we consider the following thermodynamical cycle (for the case of reduction potentials): Φg
E · Ox(g) + e(g) ∆GOx S
∆GeS Φs
E · Ox(s) + e(s)
E · Rd− (g) ∆GRd S
−
E · Rd− (s)
(1)
where “Ox” and “Rd− ” stand for the oxidized and reduced state of the redox center, respectively. Here “E” stands for enzymes, which are of our major interest, although the thermodynamical cycle is valid in general. Although Φs can in principle be calculated with the free energy perturbation approach for enzyme systems, it is often time-consuming to include a large number of solvent molecules (i.e., with periodic boundary conditions) especially if a QM/MM potential is to be used. In the alternative stochastic boundary set-up for enzyme systems11 (which has been adopted in the current work, see Secs. 2.2 and 2.3 for details), typically only a limited number of solvent molecules are included (e.g., 539 water molecules are present in the current work) and therefore the calculations give a value that nominally corresponds to the “vacuum” reduction potential, Φg (i.e., with the enzyme system in “vacuum”). To obtain the value with more appropriate treatment of the solvation effects, we employ Hess’s law to express ΦS in terms of Φg and the solvation free energy of E · Ox, E · Rd− and the electron, −
e − ∆GOx ΦS = Φg + (∆GRd S S ) − ∆GS Rd− /Ox
= Φg + ∆∆GS
− ∆GeS
(2)
QM/MM Redox Potential Calculations 55
The gas phase value (Φg ) can be obtained with free energy perturbation techniques, which will be described in Sec. 2.1 in more details. The solvaRd− /Ox , can be computed through tion term, ∆∆GS Poisson–Boltzmann type of calculations,10 which are described in Sec. 2.2. Finally, the solvation energy of the electron is taken into account by considering the appropriate reference for the reduction potential, which is 4.43 eV for hydrogen electrode.12
QM and the rest enzyme and solvent atoms are described with MM, appears to be a satisfactory compromise. With the familiar form of the QM/MM potential energy function,16 QM/MM
UTot
QM/MM
(R) = hΨ(r; R)|H QM + Hel
|Ψ(r; R)i
QM/MM
QM/MM (R) + Ubonded (R) + Uvan
+ U MM (R) 2.1.
QM/MM
Free energy perturbation with QM/MM methods
= Uele
QM/MM (R) + Uvan (R)
QM/MM
Evaluation of the reduction potential for the enzyme system in “vacuum” (Φg ) requires computing the free energy difference between the two redox states,13 Z 1 ∂G − Φg = Gg (Rd ) − Gg (Ox) = dλ ∂λ 0 Z 1 ∂U (R; λ) dλ (3) = ∂λ 0 λ
+ Ubonded (R) + U MM (R)
(6)
and the assumption that the redox center is described with QM, the hybrid potential energy at a given λ is Ox/MM
U (R; λ) = (1 − λ)[Uele
(ROx , RMM )
Ox/MM
+ Uvan+bonded(ROx , RMM )] Rd− /MM
+ λ[Uele
(RRd− , RMM )
−
where λ is a coupling parameter that converts the redox center from Ox to Rd− and the ensemble average is over the configurations of the enzyme system plus the included solvent molecules at a particular coupling strength, λ. The quantity U (R; λ) is the hybrid potential energy function in the presence of the coupling, which is usually taken in the linear form,13 U (R; λ) = (1 − λ)UOx (R) + λURd− (R) .
(4)
As a result, the “vacuum” reduction potential Φg can be cast into the following form, Z 1 dλhURd− (R) − UOx (R)iλ . (5) Φg = 0
Thus, an accurate evaluation of the reduction potential requires a satisfactory description of URd− /Ox (R) and sufficient configurational samplings to ensure the convergence of the free energy derivative at each λ. Although pure classical models based on partial charges are very fast and convenient for free energy calculations, they are not quantitatively accurate in general. Full ab initio or DFT treatment of the entire enzyme system is currently out of the question for sufficiently long simulations. Therefore, a combined QM/MM description of the potential energy function,16 in which the redox center is described with
Rd /MM
+ Uvan+bonded(RRd− , RMM )] + U MM (RMM ) .
(7)
Similarly, the free energy derivative is given in the following form, ∂G Ox/MM (R; λ) = h−[Uele (ROx .RMM ) ∂λ Ox/MM
Uvan+bonded(ROx , RMM )] Rd− /MM
+ [Uele
(RRd− , RMM )
Rd− /MM
+ Uvan+bonded (RRd− , RMM )]iλ . (8) In Eqs. (7) and (8), the dual topology approach18 was assumed in which the reactant and product appear as separate topological entities with different atomic properties (e.g., ROx and RRd− are independent). Alternatively, one may use the single topology approach,19 in which the reactant and product states appear in a single topological description and the atomic properties (such as partial charges and Van der Waal’s parameters in the general context) are transformed from one set of values to the other by the coupling parameter. As discussed in details in previous studies,18,19 both approaches have its merits and caveats. The single topology approach can
56 M. S. Formaneck et al.
be problematic if the two states tend to adopt very different conformations (e.g., different rotamer states for protein sidechains); one also has to be careful about contributions from scaled molecular parameters, e.g., Jacobian factors correspond to the change in bond lengths in the two states.20 The dual topology approach, on the other hand, might suffer from the socalled end-point problems.23 They arise from the fact that one of the states (e.g., reactant state when λ ∼ 1) will have a vanishingly small contribution to the total energy and force, which can cause large structural and positional fluctuations and thus conformational sampling problems.23 Divergence in free energy derivative can also occur due to the fast variation of certain interactions (e.g., Van der Waal’s terms) as a function of λ.24 For the current QM/MM redox calculations, the single topology approach seems to be particularly attractive because the coupling potential and the free energy derivative take simpler forms if the reduced and oxidized states are chosen to have the same Van der Waal’s and bonded parameters for the QM/MM interactions, Ox/MM
U (R; λ) = (1 − λ)[Uele −
Rd /MM
+ λ[Uele
(RQM , RMM )]
(RQM , RMM ]
Ox/MM
+ Uvan+bonded(RQM , RMM ) + U MM (RMM )
(9)
∂G Ox/MM (R; λ) = h−Uele (RQM , RMM ) ∂λ Rd− /MM
+ Uele
(RQM , RMM )iλ (10)
where only one set of coordinates, RQM , is involved for both the oxidized and reduced QM atoms. In other words, only the electronic part of the QM/MM energy [see Eq. (6)] has to be accumulated to obtain the free energy derivative. Moreover, Eq. (9) implies that no end-point problem arises because only one set of coordinate is involved for the QM part and therefore no large structural distortion can occur; i.e., there will always be contributions from one of the states that prevents the structure from falling apart. Thus the single topology approach is employed in the current work, and we leave a more detailed comparison with the dual topology approach to the future (G. Li, M. Formaneck and Q. Cui, unpublished).
Necessary modifications have been made to the simulation program CHARMM26 to allow both single and dual topology QM/MM free energy calculations in the thermodynamical integration framework.13 Since the QM/MM van der Waal’s and bonded-terms do not contribute in the single topology approach [Eq. (10)], the modification is very simple and involves only the QM/MM interface subroutines. For the dual topology approach, the reactant and product states can adopt different nuclear configurations [Eqs. (7) and (8)] and therefore includes contributions from QM/MM Van der Waal’s and bonded-terms. This requires further modifications to the REPLICA and BLOCK modules in CHARMM. The current implementation only applies to SCC-DFTB (see Sec. 3) as the QM level, although extension to other QM levels is straightforward. 2.2.
Effect of solvation
As mentioned in the beginning of Sec. 2, the proper treatment of solvation is critical to accurate evaluations of the redox potential. In addition to the fact that solvation free energies of the oxidized and reduced forms have to be calculated [Eq. (2)], additional care has to be taken to avoid artifectual structural changes or over-polarization of the QM wavefunctions due to the limited number of solvent molecules explicitly included in the stochastic boundary set-up. For such a purpose, a Poisson–Boltzmann (PB) charge-scaling technique27 developed previously for pure MM free energy simulations with stochastic boundary conditions have been adopted. In this scheme, the partial charges for the charged residues on the surface of the protein are scaled down according to a set of PB calculations to mimic the screening effect of the bulk solvent; such scaling factors tend to be large for charged residues exposed to the solvent and small in the interior of the protein (e.g., see Table 1). The QM/MM free energy calculations are then carried out with the scaled partial charges. Finally, a set of corrections for the charge scaling will be made for selected configurations; i.e., the final expression for the reduction potential for the solvated enzyme system is given by [compare with Eq. (2)] ΦS = Φg (q) + ∆∆GRd g −
Rd /Ox
+ ∆∆GS
−
/Ox
(q → Q)
(Q) − ∆GeS
(11)
QM/MM Redox Potential Calculations 57
Table 1. Scaling factors for the charged residues in the stochastic boundary set-up of the FAD-cholesterol oxidase system.a Residue
Distance (˚ A)b
Scaling
Residue
Distance (˚ A)b
Scaling
Lys 63
21.3
73.2
Glu 40
21.4
6.1
69 127
20.9 16.8
8.1 10.8
73 125
17.7 12.9
9.4 11.3
163
23.8
85.7
132
22.1
51.4
172 210
19.6 20.7
46.2 56.8
133 142
23.1 26.3
39.8 80.2
225
13.2
3.3
166
22.1
33.0
230 251
16.5 26.6
12.8 52.8
169 179
22.7 22.5
67.9 64.0
381 398
17.1 23.9
32.1 65.1
203 216
21.5 13.6
8.7 8.8
414
20.2
75.6
296
21.1
21.5
438 456
23.5 25.9
18.9 74.1
313 361
28.6 10.5
78.6 3.4
370
24.1
33.3
495
19.3
11.4
Arg 28
23.7
19.1
Asp 51
25.3
27.3
64
19.5
17.2
62
19.7
15.6
71 87
22.8 23.2
46.5 65.1
83 90
20.9 25.5
29.3 80.0
98
13.8
11.2
97
12.5
7.3
110 128
12.6 20.0
6.4 52.4
102 139
21.8 26.2
58.3 77.3
146 150
24.3 20.4
79.1 32.0
145 161
24.2 19.5
72.3 50.4
156
15.6
18.6
167
24.1
83.1
175 178
19.1 19.8
26.2 31.6
196 229
14.3 15.0
15.8 7.6
202
22.1
51.5
349
16.7
6.8
300 328
27.0 19.3
52.0 17.6
352 355
20.7 23.7
46.5 55.1
385 396
21.5 27.1
12.0 32.1
404 418
23.4 23.3
60.1 49.1
403
22.3
29.3
431
23.0
30.1
419 429
25.0 22.4
57.8 29.5
442 443
16.6 16.3
15.3 13.5
496
21.4
27.7
460
27.1
71.4
474
16.8
6.3
a The scaling factors were obtained with a set of Poisson–Boltzmann calculations in which the electrostatic potential at the active site due to a given residue is compared when the system is in vacuum and in solution (see Ref. 27 for details). In these PB calculations, only the charges on the particular residue under consideration are kept and all other charges were set to zero; a series of focusing calculations were made with the final grid size of 0.5. The dielectric constant is set to 1.0 and 80.0 for the protein and bulk solvent, respectively. The solvent probe radii and ionic strength were both set to zero. b The distance is between the center of geometry of a given residue to the N5 atom in FAD (the latter is taken as the origin of the stochastic boundary set-up).
58 M. S. Formaneck et al.
where q and Q indicate that the corresponding calculations are carried out with the scaled and full partial charges, respectively. Both charge-scaling correction Rd− /Ox (q → Q)) and solvation energy of the (∆∆Gg Rd− /Ox
(Q)) are carried full-charge structures (∆∆GS out with PB calculations. Although the magnitude of these corrections is rather small when there is little charge difference between the two states (e.g., in proton transfer reactions in triosephosphate isomerase),28 their contributions can be substantial when there is a change in the total charge of the QM region, as in redox potential calculations. Previous application of such a PB charge-scaling procedure on ligand binding properties with pure MM force field gave very stable results,39 and it is expected that this will work well for the current QM/MM free energy calculations and is essential for obtaining reliable redox potentials (see Sec. 3). 3. Application to the First Reduction Potential of FAD in Cholesterol Oxidase To illustrate the QM/MM free energy approach and the importance of solvation effects, we report preliminary results for a study of the first reduction potential of flavin adenine dinucleotide (FAD) in cholesterol oxidase (Chox). Chox is an enzyme found in bacteria, and its major function is to oxidize cholesterol for use as a carbon source in primary metabolism.30 The oxidation is made possible by FAD as the redox co-factor, which is bound non-covalently to the active site in some enzymes, while covalently attached to a His residue in others. Although many hints on the chemical mechanisms of cholesterol oxidation by Chox have been established by previous crystallographic,31 kinetic and mutation studies,30,32 a detailed picture for the catalytic pathways remains elusive. For example, it is not clear whether the oxidation reaction occurs through a step-wise radical based pathway or a concerted hydride transfer mechanism. Understanding the redox property of FAD is a key step in unraveling the catalytic mechanism of Chox. Recent mutation and X-ray studies35 showed that Asn 485, a conserved residue in the Glucose-Methanol-Choline (GMC) oxidoreductase family,36 modulates the redox property of FAD through an interaction reminiscent of the NH-π interaction found in a number of proteins.37 To quantify this contribution and systematically ex-
plore the effect of other conserved residues or water molecules in the active site, theoretical calculations are of great value. The redox process of FAD may involve the follow reactions: E · FAD(s) + e(s) → E · FAD− (s)
(12a)
E · FAD− + H + (s) → E · FADH(s)
(12b)
E · FADH(s) + e(s) → E · FADH− (s)
(12c)
−
E · FADH (s) + H + (s) → E · FADH2 (s)
(12d)
where E stands for the enzyme and (s) stands for solution. Rigorously speaking, steps (12b) and (12d) are not redox processes; however, proton transfer might be coupled with electron transfer in the experimentally measured redox potentials. In the measurements of Yin et al.,35 visible spectra clearly indicated that the product of the first reduction reaction is E · FAD− ; the product of the second reduction reaction is not clear, which could be either E · FADH− or E · FADH2 . In the current work, the major goal is to illustrate the QM/MM free energy approach, and therefore we will restrict ourselves to the first reduction potential of FAD in Chox [Eq. (12a)] and leave more detailed analysis of other processes to the future. 3.1.
Set-up of the enzyme system: Cholesterol oxidase with bound FAD
In the calculations, the X-ray structure33 for Streptomyces at 1.5 (PDB code 1B4V) with only the cofactor FAD bound to the enzyme active site was used with the standard stochastic boundary set-up11 of radius 25 centered around the flavin ring. The final model included 6175 protein atoms and 539 water molecules [Fig. 1(a)]. A Poisson–Boltzmann (PB) charge-scaling scheme27 (see Sec. 2.2) was introduced to account for solvent shielding in addition to that from the explicit water molecules in the model; a dielectric constant value of 1 was used for the protein interior in the PB calculations.39 The scaling factors are summarized in Table 1. As expected, the scaling factors tend to be large in magnitude (>50) for charged residues on the surface of the protein and tend to be small (∼5) for residues in the interior. The flavin ring of FAD was treated with QM and the rest of FAD and protein residues were described
QM/MM Redox Potential Calculations 59
Glu361
N10
Asn485
C12
C11
His447
(a)
(b)
Fig. 1. The FAD-cholesterol oxidase system based on the X-ray structure from Streptomyces at 1.5 (PDB code 1B4V). (a) Stochastic boundary set-up of 25 around the flavin ring in FAD. The FAD molecule is shown in space-filling models, several conserved residues (Glu 361, Asn 485, His 447) are shown in line forms, the protein backbone is shown in ribbon (color coded as the following: red — acidic; blue — basic; green — polar; white — non-polar), and the water molecules are shown in dots. (b) A closer view of the active site of cholesterol oxidase shown in a similar format, except that FAD is also shown in line form; several active site water molecules close to FAD are shown as red Van der Waal’s spheres. The flavin ring system of FAD is treated with SCC-DFTB, and the rest of FAD and protein atoms are described with the CHARMM 22 force field for proteins; the dash-line indicates the QM/MM boundary. The figure was prepared with VMD, W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph. 14, 33 (1996).
with the CHARMM22 force field.40 The latter include the conserved amino acids in the active site, such as Glu 361, Asn 485 and His 447 [Fig. 1(b)]; in such a way, their contribution to the reduction potential can be easily analyzed through a perturbation analysis. The QM region involved more than 30 atoms, and therefore a fast QM method has to be employed in the free energy calculations. We choose to use the SelfConsistent-Charge Density-Functional-Tight-Binding (SCC-DFTB) method,41 which was recently introduced into CHARMM in a QM/MM framework.43 It can be considered as a semi-empirical implementation
of density functional theory, and therefore offers both speed and accuracy. Benchmark calculations on several systems of biological interests, and applications to proton and hydride transfer reactions in several enzymes have been successful41–44 (also see Sec. 3.1). Link atoms were introduced between the C11 and C12 atoms in FAD to saturate the valence of the boundary SCC-DFTB atom [see Fig. 1(b)]. The link atoms interact with the MM atoms, except the “link host” MM atom (the C11 atom in this case), through electrostatic terms; no Van der Waal’s interactions are included. This scheme has been shown to be a satisfactory way
60 M. S. Formaneck et al.
to treat the QM/MM interface, particularly when the charges of the atom in the neighborhood of the link atom are small;45 this is true in the present case. For the free energy perturbation calculations, 11 windows were chosen with λ ranging from 0.0 to 1.0 with an 0.1 interval. Because the single topology approach was used, no end-point problem arises which permitted molecular dynamics calculations at λ = 0.0 and 1.0. For each window, the system was first equilibrated for 10 ps and then a production run was carried out up to 300 ps; i.e., a total of more than 3.5 ns trajectories were used for the free energy calculations. The free energy derivative was closely monitored during the calculations to ensure convergence. In the MD calculations, a time-step of 1.0 fs was used and all the bonds involving hydrogen atoms were constrained with SHAKE. Temperature was set to be 300 K with mix Langevian dynamics/Newtonian dynamics unique to the stochastic boundary procedure.11 Following the free energy calculations, Poisson– Boltzmann (PB) calculations were carried out to obtain the charge-scaling corrections and solvation free energies [Eq. (11)]27 for each window where more than twenty configurations were taken for each λ value. In the PB calculations, Mulliken charges from SCCDFTB/CHARMM calculations were used for the QM atoms. Dielectric constants of 1.0 and 80.0 were used for the protein and bulk solvent, respectively; a grid size of 0.4 was used. Much discussion has been made in the literature about the appropriate dielectric constant for the protein to use in PB calculations,46 and different values have been proposed.10 A value of 1.0 rather than the popular value of 4.0 was chosen here because conformational fluctuations (and therefore part of the protein dielectric response) were included in the MD calculations; such a value was also chosen in previous PB correction calculations in the context of ligand binding.39 We also note that the choice of using PB calculations to evaluate the charge-scaling corrections in vacuum is consistent with the fact that no cutoff was used in computing the electrostatic component of the SCC-DFTB/CHARMM interactions. Finally, perturbation analysis was used to examine the contribution of a conserved residue, Asn 485, to the reduction potential; it was found to be important in the recent measurements of Yin et al.35 One hundred and fifty configurations were taken from each window, and free energy derivatives were re-calculated
with the partial charges on the side-chain of Asn 485 set to zero; this gives a qualitative estimate on the contribution of Asn 485. Since Asn 485 is a chargeneutral residue and is very close to the redox center, the effect of long-range electrostatic interactions due to the bulk solvent is not expected to be important; Rd− /Ox Rd− /Ox i.e., ∆∆GS (q → Q) and ∆∆GS (Q) were not considered in the perturbation calculations. 3.2.
Test of SCC-DFTB for the flavin system
Although SCC-DFTB was shown to be far more accurate than AM1 for a number of hydrogen bonding systems and proton transfer reactions,41–43 its semiempirical nature requires careful tests before the application to each new system. Here we have used the flavin ring in the gas phase for such a purpose. Geometries and energies were calculated for the flavin ring in the gas phase with different protonation and redox states at AM1, SCC-DFTB and B3LYP/631+G(d, p) levels. As shown in Fig. 2 (only “FAD” and “FAD− ” geometries were shown because only the first reduction potential of FAD was studied here), the geometries are very similar at all three levels, although SCC-DFTB has a slightly smaller RMS difference in bond lengths (0.01) than AM1 (0.02) relative to the B3LYP results. The energetics, however, are rather different at the AM1 level. As shown in Table 2, SCC-DFTB gives very satisfactory results for the redox reactions in the gas phase and has a RMS error of 3.9 kcal mol−1 compared to B3LYP/6-31+G(d, p); the largest error is associated with the first redox step. By contrast, the standard AM1 approach has a much larger RMS error of 12.1 kcal mol−1 ; for the first reduction potential, AM1 has an error as large as nearly 20 kcal mol−1 . Thus, the SCC-DFTB method is a more appropriate choice as the QM level in the current QM/MM calculations. 3.3.
Free energy simulations
As shown in Fig. 3(a), the free energy derivatives are all essentially converged after 150 ps of MD simulations (excluding 10 ps of equilibration); to ensure convergence, MD simulations were extended to 300 ps when the plateau of the free energy derivative have lasted at least for 50 ps. The converged values as
QM/MM Redox Potential Calculations 61
1.219/1.237 (1.240/1.251) [1.229/1.250] 1.299/1.359 (1.305/1.364) [1.312/1.366]
O
1.381/1.398 (1.389/1.401) [1.385/1.395]
C
1.369/1.355 (1.394/1.354) [1.376/1.354]
N
C
N
1.415/1.403 (1.419/1.416) [1.413/1.405]
C
C
C
C
C
C
1.463/1.426 (1.504/1.456) [1.458/1.427] C
C O
N
1.221/1.241 (1.244/1.259) [1.229/1.252]
N
C
1.388/1.400 (1.401/1.405) [1.381/1.392]
1.309/1.333 (1.323/1.350) [1.319/1.340]
C
1.470/1.454 (1.445/1.436) [1.444/1.437]
C
Fig. 2. Optimized structures of a model flavin complex in the neutral and anion states in the gas phase at different levels; numbers without parentheses are B3LYP/6-31+G(d, p) values, and these with parentheses and brackets are AM1 and SCCDFTB values, respectively. Distances are given in Angstroms.
Table 2. A comparison of B3LYP, SCC-DFTB and AM1 for related redox processes in flavin in the gas phase.a Process FAD + e
→ FAD
B3LYP/6-31+G(d, p)
SCC-DFTB
AM1
−44.3 −333.5 −377.8 −49.4 −331.5 −380.9
−37.3 −333.8 −371.1 −50.3 −330.8 −381.1
−63.5 −314.5 −378 −54.2 −321.5 −375.7
3.9
12.1
−
FAD− + H+ → FADH FAD + e + H+ → FADH
FADH + e → FADH− FADH− + H+ → FADH2 FADH + e + H+
→ FADH2
RMS errorb
—
Values are given in kcal mol−1 , and no zero point energy has been included in the energetics. See Fig. 2 for structures. b RMS error measured relative to the B3LYP results. a
a function of λ are plotted in Fig. 3(b). It is seen that ∂G/∂λ increases in magnitude nearly linearly as a function of the coupling strength, which is reminiscent of the linear-free energy relation.48 The result can be qualitatively understood with Scheme 1, in which
the free energy of each chemical state is approximated by a parabola as a function of a collective environmental (solvation) coordinate; switching λ from 0.0 to 1.0 changes the environmental coordinate from that which favors the oxidized state (sO ) to that which favors
62 M. S. Formaneck et al.
20
(a)
0
0
-20
-20
∂G/∂λ (kcal/mol)
∂G/∂λ (kcal/mol)
20
-40 -60
-60
-100
-100 0
50
100
150 t(ps)
200
250
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ
300
-50
(c)
(d)
-55
-70
Corrected ∂G/∂λ (kcal/mol)
∆∆GSRd-/Ox(Q,λ) (kcal/mol)
-40
-80
-80
-60
(b)
-80 -90
-100 -110 -120 -130
-60 -65 -70 -75 -80 -85 -90
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 λ
Fig. 3. Results from SCC-DFTB/CHARMM free energy perturbation calculations. (a) The convergence behavior of free energy derivatives as a function of production simulation length. (b) The free energy derivative as a function of the coupling parameter, λ, from “vacuum” simulations. Note the nearly linear relationship (see text for discussions). (c) The solvation Rd− /Ox
free energy difference between the reduced and oxidized states, ∆∆GS
(Q; λ), as a function of λ. (d) The solvent Rd− /Ox
corrected free energy derivative [Eq. (11)] including the effect of the solvation (∆∆GS Rd− /Ox
(∆∆GS
[(q → Q); λ]) as a function of λ.
(Q; λ)) and charge-scaling
QM/MM Redox Potential Calculations 63 G
the trend is consistent with the free energy derivatives in Fig. 3(b). In other words, an appropriate configuration sampling is an essential step in obtaining reliable redox potentials, and free energy perturbation calculation is a robust approach in this regard.
Ox
kO(s-sO)2 2 VO λ=0
λ=1
3.4. Rd
-
kR(s-sR) 2
VR sO
sR
2
s
Scheme 1.
the reduced state (sR ). The free energy difference between the two can then be approximated as − ∂G = h∆U Rd /Ox iλ ∼ (VR − VO ) ∂λ λ kR s2R − kO s2O + (kO sO − kR sR )s + 2 1 + (kR − kO )s2 . 2
(13)
If the curvatures of the two free energy surfaces are similar, the quadratic dependence in Eq. (13) disappears and one is left with a linear dependence on s and a constant, which is qualitatively the observed trend in Fig. 3(b). Such a significant variation of ∂G/∂λ as a function of λ indicates that the protein environment of the redox center changes substantially when an additional electron is added to FAD. This suggests that it is not appropriate to calculate redox potentials based on the X-ray structure of enzymes at only one chemical state.58 As a matter of fact, we initially tried to estimate the reduction potential by energy minimization calculations for the oxidized and reduced states starting from selected protein configurations, and it was found that the results depend very sensitively on the protein conformations; it is about −10 kcal mol−1 when the conformations were collected from a MD calculation for the neutral state, but nearly −70 kcal mol−1 if the conformations were from MD trajectories for the anion state. Note that
Charge-scaling corrections, solvation free energies and reduction potential in the solvated system
Taking the free energy derivative in Fig. 3(b) and numerically integrating with the Trapezoid rule gives a free energy difference of −50.9 kcal mol−1 ; the results are −49.9 kcal mol−1 and −50.3 kcal mol−1 if only the first 150 ps and 200 ps trajectories are used, respectively, which clearly indicates a satisfactory convergence. The value of −50.9 kcal mol−1 corresponds to −2222 mV relative to the hydrogen electrode, which is very different from the experimental measurement of −498 mV.35 Such a discrepancy clearly indicates that the charge-scaling correction and solvation energies have substantial contributions, which is consistent with the consideration that the two states differ by a unit charge. This is confirmed by the Poisson– Boltzmann results [Eq. (11)]. The charge-scaling corRd− /Ox
[(q → Q); λ], is of the rection in vacuum, ∆∆Gg order of +60 kcal mol−1 and has relatively small variations as a function of λ (thus not shown in Fig. 3). The Rd− /Ox (Q; λ), by contrast, solvation correction, ∆∆Gs varies significantly as a function of λ. The value is the largest, about −120 kcal mol−1 , for λ = 0 (oxidized neutral state), and reduces to nearly half at λ = 1 (reduced anion state). In other words, the strong preference of the protein environment (which has a total charge of −6) over the neutral state at λ = 0 is partially compensated by the bulk solvent. Adding these two corrections to the free energy results in vacuum, one obtains the corrected free energy derivative. As shown in Fig. 3(d), the free energy derivative including the solvent corrections is more linear as a function of λ compared to the uncorrected values [Fig. 3(b)]. As summarized in Table 3, the integrated conRd− /Ox (q → Q) and tribution (over λ) of ∆∆Gg Rd− /Ox
(Q) are +67.5 and −96.1 kcal mol−1 , re∆∆GS spectively. It is instructive to compare the latter value to an estimate based on the generalized Born model for the solvation of spherical ions;49 using the radius of
64 M. S. Formaneck et al.
Table 3. Contributions to the calculated first reduction potential of FAD in cholesterol oxidase.a Quantity
good agreement with the experimental measurement of −498 mV for an enzyme system.
Value −1
Φg (q ) (kcal mol ) Rd− /Ox ∆∆Gg (q → Q) (kcal mol−1 )
−50.9 +67.5
Rd− /Ox ∆∆GS (Q) (kcal mol−1 )
−96.1 −983(−679)c
b
ΦS (mV)
−498
ΦExp (mV)b d
∆ΦS (mV) Asn 485
∆ΦExp (mV) N485L
e
−191 −78
a
The gas phase value (Φg (q)) was calculated with SCCDFTB/CHARMM free energy perturbation calculations; the charge-scaling correction and the solvation free energies were obtained with Poisson–Boltzmann calculations (see text) for each window and integrated over λ. The solution result (ΦS ) was obtained based on Eq. (11). b The reference was taken to be the standard hydrogen electrode. Note that the experimental value was reported relative to Safranin T in Ref. 35, which has a redox potential of −249 mV relative to the hydrogen electrode. c The value in parentheses is with an correction of 7 kcal mol−1 , which is the difference between SCC-DFTB and B3LYP/631+G(d, p) results for the electron affinity of the model flavin in the gas phase (see Fig. 2 and Table 2). d The change in the reduction potential when the partial charges on the side-chain of Asn 485 were set to zero. Only Φg (q) was considered in the perturbation calculations, because the solvation corrections are not expected to have a large contribution to the effect of a charge-neutral residue close to the active site. e The change in the reduction potential when Asn 485 was mutated to a leucine residue (from Ref. 23).
25 for the stochastic boundary set-up, and the charges for the enzyme system in the two redox states (−7 and −6, respectively), the difference in the solvation energy is −85.8 kcal mol−1 , which is fairly close to the PB results of −96.1 kcal mol−1 . Taking the two solvent correction terms Rd− /Ox Rd− /Ox (q → Q) and ∆∆GS (Q)) into (∆∆Gg account, the calculated reduction potential for FAD in Chox becomes −983 mV, which is in a closer agreement with the experimental value compared to the vacuum value of −2222 mV. Moreover, we note that SCC-DFTB has an error of 7 kcal mol−1 relative to B3LYP/6-31+G(d, p) for the first electron affinity of model FAD in the gas phase (see Table 2). Assuming the same error in the SCC-DFTB/CHARMM free energy results, our current best estimate for the first reduction potential of FAD in Chox becomes −679 mV, which can be considered to be in a rather
3.5.
Contribution of Asn 485
When the partial charges on the side-chain of Asn 485 were set to zero, the free energy derivatives were found to decrease in magnitude by about 4 kcal mol−1 , which corresponds to a negative shift in the reduction potential (see Table 3). The change is nearly invariant as a function of λ, which is consistent with the consideration that Asn 485 is very close to FAD and thus the contribution is not sensitive to the change in the protein environment or solvation effect (in terms of Scheme 1, this corresponds to that the two free energy surfaces merely have a relative shift that is independent of the “solvent” coordinate); this may not be true for a charged residue (X. Zhang and Q. Cui, unpublished). The contribution of 4 kcal mol−1 (i.e. 191 mV) due to Asn 485 is consistent with the previous estimate for such an N–H · · · π interaction37 (see Fig. 1(b) for the relative position of Asn 485 and the flavin ring in FAD). The value is also in qualitative agreement with the recent measurement of Yin et al.,35 who found that the first reduction potential of FAD shifted by −78 mV in the N485L mutant. Further comparison of the mutant X-ray structure with that of the wild type found that the side-chain of the close-by Met 122 adopted a different rotamer state in the mutant, which allowed the presence of an additional water molecule in the active site.35 Therefore, the absolute value of 78 mV is likely to be a lower bound for the contribution from Asn 485. To further verify the calculated contribution, reduction potentials for the mutant are being calculated with different X-ray structures; contributions from other residues are also being analyzed systematically with similar perturbation techniques (X. Zhang and Q. Cui, unpublished). The calculations in this section clearly illustrate that an accurate estimate of redox potential in enzymes requires the proper treatment of many aspects of the system. This includes a reliable description of the electronic structure of the redox center, its interaction with the enzyme environment, the response of the enzyme conformations to the redox reaction as well as the effect of long-range electrostatic interactions due to the bulk solvent. It is very satisfying to see
QM/MM Redox Potential Calculations 65
that after all these are taken into account, reliable redox potentials can be derived from simulations that are based primarily on “first principles” from quantum chemistry and statistical mechanics. 4. Conclusions Redox reactions play an important role in many chemical and biological processes, and therefore an accurate prediction of redox potentials in complex systems, such as enzymes, is of great importance. Although related quantities such as ionization potential and electron affinity can be estimated accurately for small to medium-size gas phase molecules,4 it is highly challenging to obtain reliable redox potentials for condensed phase systems. In the current work, a useful approach for computing redox potentials in enzymes based on free energy perturbation calculations in the QM/MM framework was developed. With an appropriate choice of the QM level and QM/MM coupling scheme, the interaction between the redox center and its environment in different redox states can be adequately described; the speed of QM/MM methods also allows sufficient configurational samplings for the convergence of free energy perturbation calculations. The method was illustrated with an application to the first reduction potential of FAD in cholesterol oxidase (Chox). With the SCC-DFTB approach as the QM level, which is both fast and fairly accurate, good agreements with recent experimental measurement were obtained provided that the effect of solvation is included properly; the latter was treated through a multi-stage charge-scaling scheme based on Poisson– Boltzmann calculations, which was developed previously for ligand-binding applications.27,39 The study of Chox also indicates that the conformational properties of the protein can change substantially during the redox process, and therefore it is essential to include a sufficient amount of configurational sampling in the calculations. For example, large errors might occur in energy minimization type of studies based only on the X-ray structure of the enzyme system in one redox state. To further improve the accuracy and stability of the redox calculations, a number of developments will be pursued in future works. First, more efficient and systematic approaches for describing solvation effects in QM/MM calculations will be
developed. These can be based on a straightforward coupling of QM/MM methods with Poisson– Boltzmann calculations,52 or with more sophisticated solvent boundary treatments.51 Moreover, a perturbation correction (rather than a constant error shift based on gas phase results used in Sec. 3.4) can be made to improve the QM/MM free energies. The fact that SCC-DFTB results are in general very close to the more sophisticated calculations (such as B3LYP) suggests that a perturbative treatment is likely to be appropriate. For the specific application to cholesterol oxidase, the current results for the first reduction potential of FAD is an encouraging start. In future calculations, other redox processes [Eq. (12)] will be studied with similar techniques developed and tested here; comparison will also be made with solution calculations. Free energy component analysis53 will be used to systematically analyze the contribution of protein residues to the reduction potentials in both the wild type and mutant enzymes. With detailed computational analysis and further experimental support, mechanisms for the regulation of FAD redox potentials in cholesterol oxidase and related flavin enzymes will emerge. Acknowledgments The current research is supported by a start-up fund from the Department of Chemistry and College of Letters and Science of University of Wisconsin, Madison. Q.C. would like to thank Dr. W. Yang for the interesting discussions, and Prof. N. Sampson for sending us the preprints of several references on cholesterol oxidase. We also appreciate the invitation from Prof. J.Z. Zhang to submit the current work to the inaugural issue of the Journal of Theoretical and Computational Chemistry. References 1. D. Voet and J.G. Voet, Biochemistry, 2nd ed. (John Wiley & Sons Inc., New York, 1995) 2. D.A. Harris, Bioenergetics at a Glance (Blackwell Science, London, 1995). 3. D.J. Waggoner, T.B. Bartnikas and J.D. Gitlin, “The role of copper in neurodegenerative diseases,” Neurobiol. Diseases 6, 221 (1999). 4. L.A. Curtiss, P.C. Redfern, K. Raghavachari and J.A. Pople, “Assessment of Gaussian-2 and density functional theories for the computation of ionization
66 M. S. Formaneck et al.
5.
6.
7.
8.
9.
10. 11.
12.
13. 14. 15. 16. 17.
18.
19.
20.
21.
22.
potentials and electron affinities,” J. Chem. Phys. 109, 42 (1998). W. Koch and M.C. Holthausen, A Chemist’s Guide to Density Functional Theory, 2nd. ed. (Wiley-VCH, New York, 2001). T. Ickiye, “Simulations of electron transfer proteins,” in Computational Biochemistry and Biophysics, eds. O.M. Becker, A.D.Jr. Mackerell, B. Roux and M. Watanabe (Marcel Dekker Inc., New York, 2001). P.J. Stephens, D.R. Jollie and A. Warshel, “Protein control of redox potentials of iron-sulfur proteins,” Chem. Rev. 96, 2237 (1996). J.M. Mouesca, J.L. Chen, L. Noodleman, D.A. Bashford and D.A. Case, “Density-functional Poisson– Boltzman calculations of redox potentials of ironsulfur proteins,” J. Am. Chem. Soc. 116, 11898 (1994). R. Konecny, J. Li, C.L. Fisher, V. Dillet, D. Bashford and L. Noodleman, “Cu, Zn-superoxide dismutase geometry optimization, energetics and redox potential calculations by density functional and electrostatic methods,” Inorg. Chem. 38, 940 (1999). B. Honig and A. Nicholls, Science 268, 1144 (1995). C.L.III. Brooks and M. Karplus, “Solvent effects on protein motion and protein effects on solvent motion: Dynamics of the active site region of Lysozyme,” J. Mol. Biol. 208, 159 (1989). H. Reiss and A. Heller, “The absolute potential of the standard hydrogen electrode: A new estimate,” J. Phys. Chem. 89, 4207 (1985). R.J. Zwanzig, J. Chem. Phys. 22, 1420 (1954). C.L.III. Brooks, M. Karplus and M.B. Pettitt, Adv. Chem. Phys. 71, 1 (1987). P.A. Kollman, Chem. Rev. 93, 2395 (1993). M. J. Field, P.A. Bash and M.J. Karplus, J. Comput. Chem. 11, 700 (1990). J. Gao, “Methods and applications of combined quantum mechanical and molecular mechanical potentials,” in Reviews in Computational Chemistry, eds. K.B. Lipkowitza and D.B. Boyd (VCH, New York, 1996), 7, 119. D.A. Pearlman, “A comparison of alternative approaches to free energy calculations,” J. Phys. Chem. 98, 1487 (1994). T. Simonson, “Free energy calculations,” in Computational Biochemistry and Biophysics, eds. O.M. Becker, A.D.Jr. Mackerell, B. Roux, M. Watanabe (Marcel Dekker Inc., New York, 2001). S. Boresch and M. Karplus, “The Jacobian factor in free energy simulations,” J. Chem. Phys. 105, 5145 (1996). S. Boresch and M. Karplus, “The role of bonded terms in free energy simulations: I. Theoretical analysis,” J. Phys. Chem. A103, 103 (1999). S. Boresch and M. Karplus, “The role of bonded terms in free energy simulations: II. Calculation of their
23. 24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
influence on free energy differences of solvation,” J. Phys. Chem. A103, 119 (1999). B. Roux, M. Nina, R. Pomes and J. Smith, Biophys. J. 71, 670 (1996). T. Simonson, “Free energy of particle insertion: An exact analysis of the origin singularity for simple liquids,” Mol. Phys. 80, 441 (1993). T.C. Beutler, A.E. Mark, R.C. van Schaik, P.R. Gerber and W.R. van Gunsteren, “Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations,” Chem. Phys. Lett. 222, 529 (1994). B.R. Brooks, R.E. Burccoleri, B.D. Olafson, D.J. States, S. Swaminathan and M.J. Karplus, Comp. Chem. 4, 187 (1983). T. Simonson, G. Archontis and M. Karplus, “Continuum treatment of long-range interactions in free energy calculations. Application to protein-ligand binding,” J. Phys. Chem. B101, 8349 (1997). Q. Cui and M. Karplus, “Triosephosphate Isomerase (TIM): A theoretical comparison of alternative pathways,” J. Am. Chem. Soc. 122, 2284 (2001). Q. Cui and M. Karplus, “QM/MM studies of the Triosephosphate Isomerase (TIM) catalyzed reactions: Verification of methodology and analysis of the reaction mechanisms,” J. Phys. Chem. B106, 1768 (2002). N.S. Sampson, “Dissection of a flavo-enzyme active site: The reaction catalyzed by cholesterol oxidase,” Antioxidants Redox Signaling 3, 839 (2001). J. Li, A. Vrielink, P. Brick and D.M. Blow, “Crystal structure of cholesterol oxidase complexed with a steroid substrate: Implications for flavin adenine dinulceotide dependent alcohol oxidases,” Biochem. 32, 11507 (1993). I.J. Kass and N.C. Sampson, “Evaluation of the role of His447 in the reaction catalyzed by cholesterol oxidase,” Biochem. 37, 17990 (1998). Q.K. Yue, I.J. Kass, N.S. Sampson and A. Vrielink, “Crystal structure determination of cholesterol oxidase from Streptomyces and structural characterization of key active site mutants,” Biochem. 38, 4277 (1999). N.S. Sampson and I.J. Kass, “Isomerization, but not oxidation, is suppressed by a single point mutation, E361Q, in the reaction catalyzed by cholesterol oxidase,” J. Am. Chem. Soc. 119, 855 (1997). Y. Yin, N.S. Sampson, A. Vrielink and P.I. Lario, “The presence of a hydrogen bond between asparagines 485 and the π system of FAD modulates the redox potential in the reaction catalyzed by cholesterol oxidase,” Biochem. 40, 13779 (2001). D.R. Cavener, “GMS oxidoreductases: A newly defined family of homologous proteins with diverse catalytic activities,” J. Mol. Biol. 223, 811 (1992). M. Levitt, and M.F. Perutz, J. Mol. Biol. 201, 751 (1988).
QM/MM Redox Potential Calculations 67
38. M.F. Perutz, G. Fermi, D.J. Abraham, C. Poyart and E. Bursaux, J. Am. Chem. Soc. 108, 1064 (1986). 39. G. Archontis, T. Simonson, D. Moras and M. Karplus, “Specific amino acid recognition by Aspartyl-tRNA synthetase studied free energy simulations,” J. Mol. Biol. 275, 823 (1998). 40. A.D.Jr. MacKerell, “All-atom empirical potential for molecular modeling and dynamics studies of proteins,” J. Phys. Chem. B102, 3587 (1998). 41. M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai and G. Seifert, “Self-consistent-charge density functional tight binding method for simulations of complex materials properties,” Phys. Rev. B58, 7260 (1998). 42. M. Elstner, T. Frauenheim, E. Kaxiras, G. Seifert and S. Suhai, Phys. Stat. Sol. B217, 357 (2000). 43. Q. Cui, M. Elstner, E. Kaxiras, T. Frauenheim and M. Karplus, “A QM/MM implementation of the self-consistent-charge density functional tight binding (SCC-DFTB) method,” J. Phys. Chem. B105, 569 (2001). 44. Q. Cui, M. Elstner and M. Karplus, “A theoretical analysis of the proton and hydride transfer in Liver Alcohol Dehydrogenase (LADH),” J. Phys. Chem. B106, 2721 (2002). 45. N. Reuter, A. Dejaegere, B. Maigret and M. Karplus, “Frontier bonds in QM/MM methods: A comparison of different approaches,” J. Phys. Chem. A104, 1720 (2000). 46. T. Simonson and D. Perahia, Proc. Natl. Acad. Sci. 92, 1082 (1995). 47. T. Simonson and D.J. Perahia, J. Am. Chem. Soc. 117, 7987 (1995).
48. A. Fersht, Structure and Mechanism in Protein Science (W.H. Freeman and Company, New York, 1999). 49. M. Born, Z. Phys. 1, 45 (1920). 50. Q. Cui, J. Chem. Phys. to be submitted. 51. D. Beglov and B. Roux, J. Chem. Phys. 100, 9050 (1994). 52. W. Im, S. Berneche and B. Roux, “Generalized solvent boundary potential for computer simulations,” J. Chem. Phys. 114, 2924 2001. 53. S. Boresch and M. Karplus, “The meaning of component analysis: Decomposition of the free energy in terms of specific interactions,” J. Mol. Biol. 254, 801 1995. 54. G. Archontis and M. Karplus, “Cumulant expression of the free energy: Application to free energy derivatives and component analysis,” J. Chem. Phys. 105, 11246 (1996). 55. A. Warshel and W.W. Parson, “Computer simulations of electron transfer reactions in solution and in photosynthetic reaction centers,” Ann. Rev. Phys. Chem. 42, 279 (1991). 56. J. Aqvist and A. Warshel, “Simulation of enzymereactions using valence-bond force-fields and other hybrid quantum-classical approaches,” Chem. Rev. 93, 2523 (1993). 57. C.N. Schutz and A. Warshel, “What are the dielectric ‘constants’ of proteins and how to validate electrostatic models,” Proteins: Structure, Function and Genetics 44, 400 (2001). 58. I. Muegge, P.X. Qi, A.J. Wand, Z.T. Chu and A. Warshel, “The reorganization energy of cytochrome c revisited,” J. Phys. Chem. B101, 825 (1997).
