This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Braun-Sand, S. Articles by Warshel, A. Search for Related Content PubMed PubMed Citation Articles by Braun-Sand, S. Articles by Warshel, A.

Studies of Proton Translocations in Biological Systems: Simulating Proton Transport in Carbonic Anhydrase by EVB-Based Models

Department of Chemistry, University of Southern California, Los Angeles, California 90089-1062

Correspondence: Address reprint requests to Arieh Warshel, E-mail: warshel{at}usc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THE CA III SYSTEM METHODS SIMULATING REPRESENTATIVE MODELS SIMULATIONS OF PTR IN... ASSESSING THE APPROXIMATIONS... CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Proton transport (PTR) processes play a major role in bioenergetics and thus it is important to gain a molecular understanding of these processes. At present the detailed description of PTR in proteins is somewhat unclear and it is important to examine different models by using well-defined experimental systems. One of the best benchmarks is provided by carbonic anhydrase III (CA III), because this is one of the few systems where we have a clear molecular knowledge of the rate constant of the PTR process and its variation upon mutations. Furthermore, this system transfers a proton between several water molecules, thus making it highly relevant to a careful examination of the "proton wire" concept. Obtaining a correlation between the structure of this protein and the rate of the PTR process should help to discriminate between alternative models and to give useful clues about PTR processes in other systems. Obviously, obtaining such a correlation requires a correct representation of the "chemistry" of PTR between different donors and acceptors, as well as the ability to evaluate the free energy barriers of charge transfer in proteins, and to simulate long-time kinetic processes. The microscopic empirical valence bond (Warshel, A., and R. M. Weiss. 1980. J. Am. Chem. Soc. 102:6218–6226; and Åqvist, J., and A. Warshel. 1993. Chem. Rev. 93:2523–2544) provides a powerful way for representing the chemistry and evaluating the free energy barriers, but it cannot be used with the currently available computer times in direct simulation of PTR with significant activation barriers. Alternatively, one can reduce the empirical valence bond (EVB) to the modified Marcus' relationship and use semimacroscopic electrostatic calculations plus a master equation to determine the PTR kinetics (Sham, Y., I. Muegge, and A. Warshel. 1999. Proteins. 36:484–500). However, such an approximation does not provide a rigorous multisite kinetic treatment. Here we combine the useful ingredients of both approaches and develop a simplified EVB effective potential that treats explicitly the chain of donors and acceptors while considering implicitly the rest of the protein/solvent system. This approach can be used in Langevin dynamics simulations of long-time PTR processes. The validity of our new simplified approach is demonstrated first by comparing its Langevin dynamics results for a PTR along a chain of water molecules in water to the corresponding molecular dynamics simulations of the fully microscopic EVB model. This study examines dynamics of both models in cases of low activation barriers and the dependence of the rate on the energetics for cases with moderate barriers. The study of the dependence on the activation barrier is next extended to the range of higher barriers, demonstrating a clear correlation between the barrier height and the rate constant. The simplified EVB model is then examined in studies of the PTR in carbonic anhydrase III, where it reproduces the relevant experimental results without the use of any parameter that is specifically adjusted to fit the energetics or dynamics of the reaction in the protein. We also validate the conclusions obtained previously from the EVB-based modified Marcus' relationship. It is concluded that this approach and the EVB-based model provide a reliable, effective, and general tool for studies of PTR in proteins. Finally in view of the behavior of the simulated result, in both water and the CA III, we conclude that the rate of PTR in proteins is determined by the electrostatic energy of the transferred proton as long as this energy is higher than a few kcal/mol.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THE CA III SYSTEM METHODS SIMULATING REPRESENTATIVE MODELS SIMULATIONS OF PTR IN... ASSESSING THE APPROXIMATIONS... CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Proton translocations (PTRs) play a major role in biochemistry in general and bioenergetics in particular (Ermler et al., 1994; Gennis, 1989; Mitchel, 1961; Okamura and Feher, 1992; Wikstrom, 1998). Molecular understanding of this issue is crucial for the elucidation of the action of ATPase (Girvin et al., 1998), bacteriorhodopsin (Luecke et al., 1999; Luecke, 2000; Royant et al., 2000; Sass et al., 2000), cytochrome c oxidase (Ostermeier et al., 1997; Yoshikawa et al., 1998), and other important systems.

The considerations of the molecular details of PTR processes have been quite challenging and controversial. Many workers followed the influential model of Nagle and co-workers (Nagle and Morowitz, 1978; Nagle and Mille, 1981) and assumed, at least implicitly, that PTR in biological systems can be described as a concerted transfer across a "proton wire" where the key control is provided by the orientation of the elements that constitute the wire (Berendsen, 2001; de Groot and Grubmuller, 2001; Kong and Ma, 2001; Law and Sansom, 2002; Murata et al., 2000; Nagle and Morowitz, 1978; Nagle and Mille, 1981; Sansom and Law, 2001; Tajkhorshid et al., 2002; Zeuthen, 2001). It is important to clarify in this respect that in Nagle's model the PTR is controlled by two steps, a HOP step where the proton is being transferred and a TURN step where the water file rearranges its orientation. It is thus assumed implicitly that the overall barrier for the PTR process involves the effect of these two steps. However, the electrostatic energy of the transferred proton was not considered and the emphasis was placed on the orientational effect (see also below). In view of the reliance on the orientational process in water, we will refer to this view as the "orientational control model". This idea is consistent with the description of proton transfer (PT) in water and ice, where all the sites are equivalent. Recent interest (Schmitt and Voth, 1998; Vuilleumier and Borgis, 1999) in the identification of the exact mechanism of H⁺ diffusion in water (the so-called Grotthus mechanism; Agmon, 1995; Eigen, 1964; Zundel and Fritcsh, 1986) has probably strengthened the focus on the proton wire concept, although the issue of the reorientation of the environment has also been considered consistently.

A different perspective has been put forward by Warshel and co-workers (Sham et al., 1999; Warshel, 1979, 1986), where the key factor that controls PT in proteins in general and PTR in particular has been identified as the electrostatic free energy of the transferred charge (the water reorganization was also taken into account but this was done while considering the system with the proton present on the donor or acceptor state). According to this view, (which is based on microscopic empirical valence bond (EVB) simulations of PTR in proteins; Warshel, 1979, 1991) the reorganization of the proton donor and acceptor (in the absence of the proton) are of secondary importance relative to the change in solvation free energy along the proton transport path. Interestingly, theoretical studies of PTR in bacterial reaction centers (RC) and cytochrome c oxidase (Kannt et al., 1998; Lancaster et al., 1996; Okamura and Feher, 1992) have implicitly recognized the importance of electrostatic effects by focusing on the pK_a values of ionizable groups and/or internal water molecules (Sham et al., 1999). An interesting, related point was made recently by Williams (2002).

At this point it might be important to clarify that the electrostatic control idea referred to the electrostatic energy of the transferred proton and not to the electrostatic energy of the water molecules that has been indeed considered by the proponents of the orientational control model. In other words, the electrostatic control model emphasizes the change of the energy of the transferred proton upon moving between different sites of the channel. It seems to us that calculations that did not consider the charge of the transferred proton could not evaluate its electrostatic barrier and thus did not contain the main element of the electrostatic model. Another aspect of the electrostatic model has been the prediction that the barrier for PTR is controlled by the electrostatic energy of the transferred proton and that the barrier is proportional to this electrostatic energy. This view is incompatible with the orientational model that emphasized the orientation of the water file without the proton. In other words, almost all those who adopted the orientational model focused, at least implicitly in their discussions and pictures, on the orientation of the water file. This was based perhaps on the assumption that once the water file is oriented correctly the transfer of the proton is not rate limiting (see discussion in "Assessing the approximations used"). Apparently, the idea that the electrostatic energy of the transferred proton is a crucial element in the control of the PTR was considered to be problematic (Nagle and Mille, 1981). More specifically, it was assumed that the use of the electrostatic energies of the proton to estimate the free energy diagrams for PTR is incorrect, not realizing perhaps the fact that G is close to G⁰ (and linearly correlated) in cases of PT with a small separation between the donor and acceptor (Warshel, 1981; Åqvist and Warshel, 1993). Instead it was suggested that the "mobility measured in ice" should be used "to estimate the rate constants for fundamental conduction processes along hydrogen-bonded chains in proteins." In this respect one may argue that the proton wire model formally includes energy states (Nagle et al., 1980) that resemble the states of the electrostatic model. However, the connection between these energies and the electrostatic free energy of the transferred proton were never made, nor evaluated. In fact, the states considered seem to correspond to orientational states rather than charge states (e.g., Fig. 2, a and b) in Nagle et al., 1980). Thus, we believe that there is a clear and fundamental difference between the orientational and electrostatic models.

View larger version (14K): [in this window] [in a new window]   FIGURE 2  Comparing the autocorrelation of the energy gap for the full model for in water (black line) and the corresponding underdamped models. The simulations were done using = 30 ps–1 and = 25 ps–1 for S/A (gray line) and S/B (gray line with ), respectively.

To examine whether the electrostatic control model is valid for general PTR processes, it is crucial to look for a system where the time dependence of the PTR is known experimentally (in aqauporin we only know that the PTR is blocked). Here one of the best benchmarks is provided by carbonic anhydrase III (CA III), because in this system we know the rate constant for PTR and its change with mutations (Silverman et al., 1993). Reproducing the observed trend by microscopic simulations can provide a very powerful way of distinguishing between different conceptual models. If, for example, we find that the PTR in CA III is controlled by the electrostatic energy of the transferred proton rather than the water orientation, we are one step forward. In view of this consideration we will focus here on CA III. We note, however, that this study addresses the general issue of PTR in proteins and simply considers CA III as a powerful benchmark. We believe that an approach that reproduces the rate in CA III from first principles will provide the correct tool for analyzing PTR in other systems. Conversely we believe that methods that cannot reproduce the PTR in CA III are not capable of modeling PTR in other proteins and channels (as long as these have a significant barrier). We also note that the PTR in CA III involves a chain of several water molecules and other donors and acceptors, thus making it relevant to other proton-conducting systems. Here it is significant that the free energy barrier for PTR in CA III is 10 kcal/mol, which is much lower than the 20 kcal/mol barrier in aquaporin (Burykin and Warshel, 2003) (which has been considered originally to reflect the water orientation effects) and quite close to the estimated 6 kcal/mol barrier in gramicidin A.

To obtain a molecular understanding of PTR in proteins and to resolve the above controversies, it is essential to be able to develop simulation methods capable of converting the available structural information about relevant systems to the corresponding observed kinetics. The method of choice for actual simulations is probably the EVB method. This method was developed in 1980 for studies of chemical reactions in proteins (Warshel and Weiss, 1980) and has been used extensively by our group and by others for simulations of PT in proteins and solution (for a recent review, see Warshel, 2003). It may be important to note that recent EVB versions (for a list of versions, see Florian, 2002) have basically the same ingredients as the original EVB model. This includes the MS-EVB model (Schmitt and Voth, 1998, 1999; Vuilleumier and Borgis, 1997, 1998) used in recent attempts to study PTR in proteins (Ilan et al., 2004). More specifically, the original EVB model (Warshel and Weiss, 1980) and many of its subsequent applications used multistate models and the main point, adopted in all EVB studies of reactions in proteins, is the introduction of the electrostatic potential of the environment in the EVB Hamiltonian. As clarified in the Methods section, this is not much different than, for example, an early MS-EVB study of PTR within six water molecules in water (Vuilleumier and Borgis, 1997). At any rate, the EVB provides a general combined quantum mechanical/molecular mechanics (QM/MM) method that can be applied for basically any process in condensed phases. Moreover, the EVB can be used and has been used extensively in quantitative studies of activation free energies for PT in proteins (e.g., Åqvist and Warshel, 1992, 1993; Warshel, 1991, 2003). The problem, however, is to evaluate the actual kinetics (and/or proton current) in cases with high activation barriers and many binding sites for the proton. In such cases it is practically impossible to use direct microscopic EVB simulations to reproduce the time-dependent PTR process, because this process can occur on a millisecond timescale. Thus it is important to develop related but simpler approaches.

Recently we developed semimacroscopic approaches that should allow one to simulate general PTR processes in proteins (Sham et al., 1999). This approach starts with the EVB formulation, which is then reduced to a simplified expression (a modified Marcus' formula) for the activation free energy for the PTR processes. The free energy changes needed for this treatment are then evaluated by semimacroscopic electrostatic calculations and the resulting calculated activation free energies and rate constants are used to evaluate the relevant kinetics by a master equation. However, this modified Marcus' approach involves several approximations that have not yet been verified by careful comparison to the corresponding microscopic results. Considering the effectiveness of the modified Marcus' approach, it seems to us that it is crucial to establish its general validity. It is also important to try to move from the master equation treatment to a more "molecular" time-dependent treatment that will also be closer to the full EVB treatment, but still allow us to study long-time PTR processes.

In view of the above considerations we developed in this work an approach that is halfway between the full EVB model and the modified Marcus' treatment. This model evaluates the effective free energy surfaces in terms of simplified EVB surfaces and then compares the results of Langevin dynamics (LD) on these surfaces with those obtained from explicit all-atom molecular dynamics (MD) simulations on actual EVB surfaces. The performance of our model is examined first by simulating PTR along a linear chain of EVB water molecules in water. It is shown that the rate of the PTR process depends exponentially on the free energy of the proton at the highest point of the free energy profile of a stepwise transfer process. With this finding in mind we move to the main challenge, which is the simulation of the PTR in CA III. It is found that the simplified EVB model reproduced the observed rate constant for PTR in CA III and in some mutants without adjusting any parameters for the protein calculation. This seems to establish the validity of our model and its applicability for general studies of PTR in proteins.

	THE CA III SYSTEM

TOP ABSTRACT INTRODUCTION THE CA III SYSTEM METHODS SIMULATING REPRESENTATIVE MODELS SIMULATIONS OF PTR IN... ASSESSING THE APPROXIMATIONS... CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
As stated in the introduction it is important to have a well-defined benchmark for studies of PTR processes in proteins. Unfortunately there are very few well-characterized systems where we know the starting initial and final position of the transferred proton and the corresponding transfer rate. Even in gramicidin A, it is not entirely clear what is the rate-limiting step of the PTR process (Decoursey, 2003). Fortunately, as will be shown below, CA III provides an excellent benchmark with clear structural information and well-studied rate constants for the PTR process (Silverman, 2000).

The catalytic reaction of CA III can be described in terms of two steps. The first is an attack of a zinc-bound hydroxide on CO₂ (Åqvist and Warshel, 1993; Silverman and Lindskog, 1988):

(1)

The reversal of this reaction is called the "dehydration step". The second step involves the regeneration of the OH^– by a series of PT steps (Åqvist and Warshel, 1992; Silverman and Lindskog, 1988).

(2)

where (in the notation of Silverman et al., 1993), BH⁺ can be water, buffer in solution, or the protonated form of Lys-64 (other CAs have His in position 64). A groundbreaking study of Silverman and co-workers (Silverman et al., 1993) has determined the k_B for the native enzyme and its mutants. This provided an extremely clear benchmark for PTR in proteins. Obviously, a model that can reproduce the rate constant for PTR in this system is likely to provide a useful tool for general studies of PTR in biological systems. Conversely, models that cannot accomplish this task are not likely to describe correctly PTR in proteins. Here one of the obvious questions is whether the rate constant for PTR in proteins is determined by the free energy of the transferred proton or by the orientation of the unprotonated water molecules (note that the electrostatic free energy of the proton reflects the reorganization of the protein plus water system).

The subsequent sections will examine our ability to model the PTR in CA and will also address the general implications of our findings. It is important to note here that the CA III system is highly relevant to other systems that "conduct" protons because it includes several well-defined water molecules that have been considered as a proton wire (e.g., Cui and Karplus, 2003; Isaev and Scheiner, 2001; Silverman, 2000). In addition to the benefit of using CA III as a model for PTR in proteins, it is also useful to understand the nature of the proton shuttle in CA III to gain a better understanding of the action of this enzyme.

	METHODS

TOP ABSTRACT INTRODUCTION THE CA III SYSTEM METHODS SIMULATING REPRESENTATIVE MODELS SIMULATIONS OF PTR IN... ASSESSING THE APPROXIMATIONS... CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

 
Our initial goal is to provide a general potential surface for PT between all the relevant sites in a given biological system. Here we believe that currently the most effective way to describe PT in condensed phases is the EVB method (Warshel, 1991). Alternative ab initio QM/MM methods are still too slow to allow for reliable simulations of general PTR (see discussion in Warshel, 2003). The EVB method was used in many works and by many research groups, and its reliability has been demonstrated in studies of PT in proteins and solution (e.g., Åqvist and Warshel, 1992, 1993; Bala et al., 2000; Billeter et al., 2001; Feierberg and Åqvist, 2002; Hwang et al., 1988; Kim et al., 2000; Kong and Warshel 1995; Schmitt and Voth, 1998; Vuilleumier and Borgis, 1998a,b; Warshel and Weiss, 1980; Warshel, 1991, 2003). Thus we will not give here a description of the EVB method and refer the reader to the available literature and book (Warshel, 1991). We also would like to point out that the EVB program is implemented in the MOLARIS program (Chu et al., 2004). However, we will describe here the EVB for the case of n protonation sites (which are considered formally as bases) and one excess proton. In this case we describe the EVB quantum system in terms of diabatic states

(3)

where B_iH⁺ is the protonated form of the B_i protonation site (e.g., an H₃O⁺).

Now, the i-th diagonal element of the Hamiltonian of this system is described by a force-field-like function that describes the bonding within donors, the bond of the proton to the i-th base, as well as the nonbonded interactions in the system and its interactions with the surroundings (protein or water). More specifically, the diagonal elements are described by

(4)

where the b⁽ⁱ⁾ values and ⁽ⁱ⁾ values are, respectively, the bond lengths and bond angles in the quantum mechanical system composed of the n bases and the excess protons. The r_k,k' runs over all the nonbonded distances in the quantum system, and the q values are residual charges of the indicated atom in the given resonance state. The screening of the charge-charge interaction is introduced only for short distances (bonding distances) where the correct quantum mechanical description does not follow the classical prescription. The r^–4 term represents an approximation for the inductive interaction between the solute charges. describes the interaction between the quantum system (the "solute") in its i-th state and the surrounding classical system (the "solvent") that includes water molecules and/or protein atoms. is the solvent-solvent classical potential surface. Finally, is the so-called "gas phase shift" that determines the relative energy of the diabatic states (Warshel, 1991). The off-diagonal elements (the H_ijs) are described by empirical functions that are fitted to experimental information and ab initio calculation, and the ground-state energy, E_g, is obtained by diagonalizing the EVB Hamiltonian.

(5)

More detail on the nature of the EVB matrix elements for different systems is given elsewhere (e.g., Åqvist and Warshel, 1993; Warshel, 1991). In this work we represent the cases where the bases are water molecules with a slightly modified version of the parameters used in reference (trajbl et al., 2002). In other cases (i.e., His, Lys, Asp, and Glu as proton acceptors) we adjusted previously used EVB parameters.

The EVB Hamiltonian represents the basis functions of Eq. 3, where the proton can be located on any of the sites of the "active space" (the part of the system that is described quantum mechanically). The rest of the system is described classically. At this point, it might be useful to clarify that the EVB and the so-called MS-EVB (Schmitt and Voth, 1998; Vuilleumier and Borgis, 1998a,b) that were so effective in studies of proton transport in water, are more or less identical. More specifically, the so-called MS-EVB includes typically six EVB states in the solute quantum mechanical (QM) region and the location of this QM region changes if the proton moves. The QM region is surrounded by classical water molecules (the molecular mechanics (MM)) whose effect is sometimes included inconsistently by solvating the charges of the gas phase QM region (this leads to inconsistent QM/MM coupling with the solute charges as explained in, e.g., Shurki and Warshel, 2003; Villa and Warshel, 2001). More recently the coupling was introduced consistently by adding the interaction with the MM water in the diagonal solute Hamiltonian. Now our EVB studies were performed repeatedly with multistate treatment (e.g., five states in Warshel and Russell, 1986) and this has always been done with a consistent coupling to the MM region. Thus the only difference that we can find between our EVB and the recent MS-EVB version is that our EVB studies did not update the location of the QM region during the simulations, because they dealt with processes in proteins, where the barrier is high, rather than with low barrier transport processes (so that the identity of the reacting region has not changed during the simulations. Also note that the MS-EVB simulations in proteins, which involve a limited number of quantum sites, do not change the QM region (e.g., Cuma et al., 2001) during the simulations. Finally, in cases of high barriers the change of the QM region is not useful, because the main issue is the ability to obtain a proper evaluation of the free energy associated with climbing the barrier. Thus we conclude that the EVB and MS-EVB are identical methods, although we appreciate the elegant treatment of changing the position of the QM region during simulations, which is a very useful advance in EVB treatments of processes with a low-activation barrier.

At this point, we would like to clarify that we view the full EVB treatment as a reliable, fully documented, and established method for simulating PTR in condensed phases. The MOLARIS program (Chu et al., 2004; Lee et al., 1993) can perform EVB calculations on any biological system. Thus the validity of the EVB is not the issue of this work. What we need here, however, is a simplified EVB-based model that will reproduce the results of the full EVB model (which treats the protein and solvent explicitly) in a much shorter simulation time, and allows for simulations of PTR processes in the microsecond range. The simplified model can be based on an effective potential that treats explicitly only the molecules that are involved directly in the PTR (the active space). In doing so, our first priority is to force the free energy functionals associated with each simplified EVB state (the microscopic Marcus' parabolas), at the given environment, to have the same minimum and similar curvature as the corresponding functionals of the full EVB model. In other words, using the full EVB treatment we describe the adiabatic ground-state free energy, g, for a PT between two sites by considering the free energy functionals (g₁ and g₂) associated with the diabatic energies (the ₁ and ₂) and mixing them by the off-diagonal element H₁₂. The results of such a treatment are described in Fig. 1 (the figure also defines the reorganization energy, ₁₂, which will also be discussed below). Now when we construct the simplified EVB model we try to force its free energy functionals to overlap those obtained by the full model. This can be done quite effectively by representing the simplified EVB surface by the same type of solute surface as in Eq. 4, while omitting the solute-solvent and solvent-solvent terms (the U_Ss and U_ss terms) and modifying the ⁽ⁱ⁾ values to reflect the missing effect of the solvent. With this in mind we use here an effective Hamiltonian, of the form

(6)

The free energy, associated with the ground-state potential () of the simplified Hamiltonian, (here the free energy accounts for the average over the coordinates of the active space) is treated as the effective free energy surface that includes implicitly the rest of the system. In other words, we use

(7)

where r are the coordinates of the active space. Now the treatment of Eq. 6 provides a simple and powerful way of representing the energetics of the environment implicitly by adjusting the ⁽ⁱ⁾ values to the corresponding ⁽ⁱ⁾ values while imposing the requirement

(8a)

(8b)

where ( )_eff represents the quantity obtained with the effective EVB potential and ( )_complete designates the results obtained when the EVB treatment considers the entire system explicitly. For convenience we usually determine (and the corresponding ⁽ⁱ⁾ values of the effective model) by the semimacroscopic electrostatic calculations described below. At this point we would like to clarify that what is done to satisfy Eq. 8a is not some ad hoc treatment, but a well-defined procedure that exploits the fact that the EVB free energy functionals, g, change linearly with the gas phase shift (the ⁽ⁱ⁾). To satisfy Eq. 8b we need to have a similar curvature to the EVB functional of the simplified and the full model. Here we have two options. The first is to modify the solute potential (usually it is enough to change the D in the ⁰ of Eq. 4) in the simplified model, and thus to change the solute reorganization energy to account for the effect of the reorganization energy of the missing solvent. The second option is to add effective solvent coordinates to the simplified model (both options will be used here).

View larger version (28K): [in this window] [in a new window]   FIGURE 1  A description of the relationship between the diabatic free energy functions (g1 and g2), the adiabatic ground state free energy (g), the reorganization energy (), and the mixing term (H12). The specific results are taken from EVB simulations of a PT step in CA III.

Finally, to obtain a stable result we have to consider one more step. That is, the calculation of by the EVB/umbrella sampling (EVB/US) approach (Warshel, 1991) or by the related free energy perturbation approaches can be quite challenging due to convergence problems. Instead, we chose to use the semimacroscopic PDLD/S-LRA method (Lee et al., 1993; Schutz and Warshel, 2001) used effectively in our previous studies of PTR and ion transport (Burykin and Warshel, 2003; Sham et al., 1999). This approach describes the total free energy of each proton configuration as

(9)

where is the charge of the k-th group in the m-th configuration. The relevant pK_a values in the different sites of the protein are determined semimacroscopically by evaluating the change in solvation energy of moving the charged group from water to the given site (for a review of this approach, see Schutz and Warshel, 2001). Now the free energy of a proton transfer between site i and site j is given by Warshel (1981):

(10)

where is the change in the effective charge-charge interaction between the donor and acceptor upon PT and (R) is the effective dielectric constant for this interaction (at a charge separation R).

With the effective potential defined by Eqs. 6 and 7 above it is possible to examine the time dependence of PTR processes by Langevin dynamics (LD) simulations. That is, the time dependence of the system is determined by a Langevin equation (McQuarrie, 1976):

(11)

where is the effective potential of Eq. 7, i runs over the ions, runs over the x, y, and z Cartesian coordinates of each ion, is the mass of the i-th ion, _i is the friction coefficient for the i-th ion, and A_i is a random force, which is related to _i through the fluctuation-dissipation theorem (Kubo and Toyozawa, 1955).

The above LD treatment represents the effect of the solvent by modifying the solute reorganization energy and by the friction coefficient of Eq. 11. A more physically consistent treatment should include also an explicit effective solvent coordinate. The prescription for doing so is well known in the electron transfer community (e.g., Dakhnovskii and Ovchinnikov, 1986; Warshel and Parson, 2001), where the solvent coordinate, Q, is considered as the corresponding reaction field, or more microscopically by the electrostatic component of the energy gap (Hwang et al., 1988; Warshel, 1982). This type of solvent coordinate has been introduced to general studies of charge transfer reactions in the framework of a two-state EVB formulation (Hwang et al., 1988). For the benefit of the reader, we outline some elements of this treatment, but recommend reading (Hwang et al., 1988) for the details. Now in the two-state model, one can write

(12)

where the r and q are the actual "normal" coordinates of the solute and solvent molecules, whereas R and Q are the effective dimensionless coordinates for the solute and solvent, respectively. The effective frequencies _Q and _R are evaluated by in which P() is the normalized power spectrum of the corresponding contribution to (₂ – ₁). The R are related to the regular solute reaction coordinate R' = (b₁ – b₂) by where µ_R is the reduced mass for the normal mode that is the compression of b₁ and extension of b₂. The reaction coordinate Q is defined in terms of the electrostatic contribution (₂ – ₁)_el to (₂ – ₁). Thus, we have Q = –(₂ – ₁)_el/_QQ), which is also related to the dimensional solvent reaction coordinate, Q', by Q = Q'(_Qm_Q/)^1/2. V₀ is the difference between the minima of ₂ and ₁. Here we replace the contribution from each set of coordinates by one effective coordinate. The displacements (the s) are related to the reorganization energy, , given by

(13)

An equivalent and more familiar definition of the solvent coordinate (see also Eq. 12) can be obtained in terms of the macroscopic reaction field (_R) at the solute cavity, i.e., taking Q to be proportional to _R we obtain

(14)

where µ₁ and µ₂ are the dipole moments of the solute for the corresponding diabatic states, and Q = _||/C (here _|| is the projection of on (µ₁ – µ₂)). In the simple two state case we can obtain a convenient description of the dependence of the ground-state potential (Hwang et al., 1988) on the solvent coordinate using:

(15)

where µ is the adiabatic dipole moment of the solute. The derivatives of the ground-state potential with respect to R and Q gives a pair of coupled LD equations for the solute and solvent systems (see Hwang et al., 1988). Here we would like to adopt the above treatment, but we have to extend it to a multistate model that describes a chain of water molecules. This is done by writing

(16)

where is defined in Eq. 9 and where we assign two solvent coordinates to each state that corresponds to a proton inside the chain. In this way we have one coordinate (Q_i,k(i)) assigned to the i-th oxygen (which is protonated in state i) and to the subsequent oxygen, while assigning a second coordinate (Q_i,k'(i)) to the i-th oxygen and the oxygen before it. When i is the first oxygen, we do not have a Q_i,k'(i). In the case of only one pair of oxygens (only Q_i,k(i)), we have the same treatment as in Eqs. 12 and 15, except that we use instead of /2, for reasons that will be clarified below. Now when we consider the total energy of the system we can write

(17)

is the lowest eigenvalue of Eq. 5, and the B term expresses the coupling between the solvent coordinates. Here reflects the effect of the _is of Eq. 16, and the other terms in the equation reflect the cost in solvent-solvent energy associated with moving the solvent coordinates from their equilibrium position. This treatment is slightly different than the one used by (Hwang et al., 1988), because in that treatment the solvent-solvent repulsion was already incorporated in the values. As a result, this is given by Note that this relationship is different (by 2) than the standard relationship, due to the presence of the solvent term in Eq. 17.

Substituting Eq. 16 in Eq. 5 we obtain

(18)

where is the solute Hamiltonian (with as diagonal elements and the gas phase as off-diagonal elements) and C_g is the ground state eigenvector of the simplified EVB Hamiltonian (with the diagonal elements of Eq. 16). The corresponding LD equation for the solvent coordinate is now expressed as

(19)

where , whereas and m_Q are the effective friction and effective mass of the solvent. Here we replace the coupling between the solute dipole and the solvent reaction field (the term of Eq. 15) by the alternative expression, which is exact in the two-state case. Now, determines the amount of positive charge on the i-th site, and also plays the role of In the multistate case we introduce a significant approximation and consider independently the coupling of the solvent to each charge.

The parameter in Eq. 16 can be obtained from the power spectrum of the solvent fluctuations (see Hwang et al., 1988) and is obtained from the electrostatic contribution to the solvent reorganization energy, The friction coefficient, can be obtained (see Hwang et al., 1988) from the relationship where

(20)

In this work we considered two models; model S/A (where S designates simplified) only includes the solute coordinate (Eq. 11), whereas S/B considers the solute and solvent coordinate (Eqs. 11 and 19) simultaneously. Both models should reproduce the main dynamical and energetic features of the PTR process. Here we note that the dynamics of the PT processes are controlled by the fluctuation of the energy gap (Hwang et al., 1988; Warshel, 1984; Warshel and Parson, 2001). Thus we require that the autocorrelation of the energy gap for the active and complete systems will be similar. Thus, we required that

(21)

We also required that the effective dynamics of the proton transport process will be similar (similar residence time) in both the simplified and full models. For model S/A we obtained the best fit using and a.u. The performance in terms of satisfying Eq. 21 is described in Fig. 2, whereas other aspects are considered in "Simulating representative models". The use of a large for H reflects the fact that this model tries to account simultaneously for the missing solvent fluctuations by the effective dynamics of the solute. In looking for the optimal parameters for model S/B we started by fixing the solute and considering only the solvent coordinate for both the simplified and full model. This gave with obtained from the power spectrum of the electrostatic energy gap of the full model for the solvated H₃O⁺ H₂O system (see discussion in Hwang et al., 1988), and obtained from the autocorrelation of the solvent contribution to the energy gap (Hwang et al., 1988). The fitting of both models also gave us a solvent reorganization energy, of 30 kcal/mol. The solvent effective mass was estimated by using the relationship that gave The values of and were then obtained by trying to satisfy Eq. 14 for both the total energy gap, where both the solute and solvent are free to move, and also by requiring that the time dependence of the proton transport (see below) will be similar in both models. A reasonable fit was obtained with and (see Fig. 2). Furthermore, we obtained an excellent fit between the _Q of the full model and model S/B.

In general, one could have tried to obtain different friction coefficients for different atoms. However, we obtained similar results while changing the values of the heavy atoms, and the next important improvement is the inclusion of the solvent coordinates. Another important issue is the friction used for the protein simulations. Here we used the same friction as in the water simulations, because the autocorrelation of the energy gap in the protein and in solution were found to be similar.

We would like to emphasize at this point that the exact value of the optimal dynamical parameters is not so crucial here because different parameters allow us to satisfy Eq. 21. More importantly, our main aim is to see if models that reproduce the approximated trend of the dynamics of the full model can be used to explore the dependence of the PTR rate on the free energy of the central proton. If we can establish a Boltzmann type dependence of the PTR rate on the electrostatic free energy of the transferred proton (with a reasonable range of parameters for the simplified model), and if we can capture such a trend with the full model (in the time range available for this model) we can gain some confidence in the validity of the electrostatic model. On the other hand, if we find "conduction band"-like PTR processes even with the large electrostatic barrier, then the orientational model is valid. It is also useful to note that it is harder to fit simultaneously the free energy surfaces and the dynamics of model S/A and the full model, than to do so with model S/B and the full model. Here the obvious difference is the fact that model S/B provides a much more realistic description of the solvent coordinate than model S/A. Nevertheless, having two different models allows us to be more certain about the validity of our conclusion.

The simulations that used the above parameters and the LD equations (Eqs. 11 and 19) were found to be too slow to allow us to obtain the relevant proton current in cases with high activation barriers. In such cases it might be advantageous to treat the system in the overdamped limit at the expense of having a large transmission factor. Thus we decided to use the overdamped LD approach (Schuss, 1980; van Gunsteren et al., 1981) of model S/A in studies of cases with high activation barriers. Now we obtained the optimal behavior (in terms of stability of the simulations) using For model S/B, we found the underdamped model to be sufficiently effective. In this case we took into account the incorrect "dynamical" effects by using the relationship

(22)

The scaling factor, S, reflects the change in transmission factor upon change from the underdamped to overdamped model (with the specific _H and ). The dependence of this factor on the activation barrier for the transfer process was evaluated using a model of six water (a, b, c, d, e, f) molecules where the pK_a of the two central molecules is being changed artificially to establish the desired activation barrier, g. The ratio between the transmission factor of the two models was then determined, according to the prescription of Hwang et al. (1988), by running downhill trajectories (starting with a proton at the center of the system) and evaluating the average time of relaxing to H₂O_b or H₂O_e. This procedure gave S 16 for g 0, and S 7 for g > 4. This means that for the challenging high barrier case the overdamped and underdamped models gave similar rate constants, but the overdamped model requires less integration steps.

To examine actual PTR in proteins, we considered the K64H-F198D mutant of CA III as a model system (see Schutz and Warshel, 2004; Silverman et al., 1993; Silverman, 2000 for a discussion of this system). The starting configuration was taken from the x-ray structure of S-glutathiolated carbonic anhydrase III (Protein Data Bank identification code 1FLJ; Mallis et al., 2000) that was then "mutated" to the K64H-F198D mutant and relaxed in a long relaxation run (100 ps) while being subjected to distance constraints of 9 Å between the Zn⁺² ion and the N nitrogen of His-64. This generated an "in" conformer for the histidine with a shorter chain of water molecules between the Zn⁺² ion and the histidine, and thus simplified the simulation and discussion. To construct the active space (the part included explicitly in the simplified models) we started by our regular procedure of embedding the protein in water and running a relaxation run (Lee et al., 1993). Next, we located all internal water molecules that can be involved in the PTR process and the protein groups that can participate in this process. Next we examined two options. In the first one, we kept all the internal water molecules and evaluated pK_a values of their protonated form by the PDLD/S-LRA approach. We then looked for the shortest network between the molecules with the highest pK_a (H₃O⁺). These molecules were then used as the active space of the simplified model with a Cartesian position restraint with K = 10 kcal·mol^–1·Å^–2. We also added a distance restraint of the form with K = 4 kcal·mol^–1·Å^–2 between the oxygens of nearby water molecules. In the second procedure we kept in the active space all internal water molecules that were a reasonable distance from each other, and then again, used the Cartesian restraint V' upon moving to the simplified model. It was eventually concluded that the first model is sufficiently reliable for this set of calculations. The LD simulations were done with our MOLARIS simulation program (Chu et al., 2004; Lee et al., 1993), considering the given chain of donor and acceptors. The optimal time steps were 20 fs and 1 fs for the overdamped and underdamped models, respectively. All the PDLD/S-LRA calculations, including the MD simulation needed to generate the protein configurations for the LRA calculations, were performed using the MOLARIS package. Similarly, we used MOLARIS to perform the EVB calculations with the complete model. The EVB parameters for the full model were taken from our previous EVB studies of the autodissociation of water in water (trajbl et al., 2002) with minor modifications (see Table 1). The EVB parameters of the simplified models were similar to those of the full model, and the main changes (see Table 1) were: i), the reduction of D in the Morse potential of the O-H stretch (this was done to reduce the solute reorganization energy in model S/A so that the total reorganization energy would be similar to that of the full model, which includes the solvent. We also used, for convenience, the same D in model S/B, although here we could use higher D and different H₁₂. ii), The charges of the simplified model were reduced because this model does not include explicitly the screening effect of the surrounding environment.

(23)

where x is the generalized reaction coordinate, which is given by Now we can exploit the fact that the curves can be approximated by parabolas of equal curvatures (this linear response relationship was found to be valid by many microscopic simulations (e.g., Åqvist and Warshel, 1993) and write

(24)

where the reorganization energy, , (which is illustrated in Fig. 1) includes now both the solute and solvent contributions. Using Eqs. 23 and 24, one obtains our modified Marcus relationship (Åqvist and Warshel, 1993; Hwang et al., 1988; Warshel, 1991)

(25)

where is the free energy of the reaction, and H_ij is the off-diagonal term that mixes the two relevant states whose average value at the transition state are and respectively. The first term in this expression is the regular Marcus' equation (Marcus, 1964), which corresponds to the intersection of and at The second and third terms represent, respectively, the effect of H₁₂ at and

It is useful to point out that the same approach used to derive Eqs. 23 and 25 can be used to derive an expression for a concerted path in a three-state system. This can be easily done numerically by a 3 x 3 EVB equation with identical reorganization energies and with

(26)

where the coordinate system is defined in Fig. 3. This type of system will be discussed at the end of the next section.

View larger version (92K): [in this window] [in a new window]   FIGURE 3  A two-dimensional representation of the free energy surface for a concerted PT. The figure considers the EVB results for the system as a function of the O-H and Oß-Hß distances. The spacing of the contour lines is 3 kcal/mol. The figure demonstrates that the concerted path () does not provide a lower barrier than the stepwise path (ß) in the typical case where and

	SIMULATING REPRESENTATIVE MODELS

TOP ABSTRACT INTRODUCTION THE CA III SYSTEM METHODS SIMULATING REPRESENTATIVE MODELS SIMULATIONS OF PTR IN... ASSESSING THE APPROXIMATIONS... CONCLUDING REMARKS ACKNOWLEDGEMENTS REFERENCES

View larger version (67K): [in this window] [in a new window]   FIGURE 4  Comparing the time dependence of the probabilities of the different EVB states, (the ) of the (H2Oa, H2Oc, H2Od, H2Oe, H2Of) model system (where Gß = 0 kcal/mol for all PT steps) using the full model (a), the simplified model (S/A) in the underdamped (b) and overdamped limits (c). When a proton is localized on a given site, the corresponding is unity.

View larger version (66K): [in this window] [in a new window]   FIGURE 5  Comparison of the time dependence of the probability amplitudes of the transferred proton for the (H2Oa, H2Oc, H2Od, H2Oe, H2Of) model system, with Gbc = Ged = 2 kcal/mol, using the full model (a), the simplified model (S/A) in the underdamped limit (b) and the overdamped limit (c). The coordinates are the same as those described in Fig. 4.

bd

be

(27)

where we used the fact that