Metalloenzymes catalyze many difficult and critical reactions in biological systems.¹ The ability of these biological catalysts to perform these reactions under mild conditions with selectivity and stability of the catalytic center generally exceeds chemists' current capabilities in developing related synthetic catalysts. Of known proteins, about 40% contain metal ions based on recent analysis of the protein data bank (PDB).² Fe, Zn, Mn, and Cu are the transition metals that occur most frequently among the catalytic sites in metalloenzymes, and Mg is prominent among the nontransition metals.^1,2

Fe has been implicated in early electron transfer going back to prebiotic processes of geological and perhaps biological significance, including the light driven generation of H₂ gas which could serve as reducing equivalents for organic synthesis, other direct oxidative or reductive reactions catalyzed by ferrous iron in water and by iron-oxo (or hydroxo) and iron-sulfur minerals, and formation of banded iron-oxo deposits.³ For example, it has been shown experimentally that carboxylic acids can be generated by reduction of CO₂ by magnetite (Fe₃O₄) but not by Fe(II) in water alone.⁴ Similarly, a stoichiometric reaction starting from FeS and H₂S generates FeS₂ (pyrite) and HCOO^- + H₂O from bicarbonate (HCO₃^-).^5,6 Such stoichiometric reactions can be considered as precursors to catalytic reactions once a mechanism is found for regenerating the initial form of the catalyst.

Fe, Mn, and Cu are redox active and can bind to a variety of different ligands. In biological systems, metal ion oxidation states Fe(II,III,IV), Mn(II,III,IV), and Cu(I,II) are most prominent; the spin states of Fe and Mn are also variable, with high-spin mainly preferred.¹ In many dinuclear or polynuclear enzymes and related electron-transfer proteins, these metal sites can also be spin-coupled.⁷ By contrast, Zn(II) is redox inactive, although its Lewis acid character allows it to participate in many reactions involving charge flow.⁸ In some Zn-Fe and Cu-Zn enzymes, metal sites are linked electronically by bridging groups.^9,10 Mg(II) is a hard divalent cation with a strong affinity for water and phosphate coordination, and a catalyst for many hydrolysis reactions.¹¹ All of these metals, in the 2+ or higher oxidation states, can promote the acidification of H₂O and for Fe(III,IV), Mn(III,IV) of OH^- as well.¹² Some of the main ligand players on this stage are carboxylates (Asp,Glu), thiolate (Cys), imidazole (His), phenol (Tyr), and guanidinium (Arg), which are highly versatile. The first four can undergo protonation/deprotonation, while Arg can act as a cationic assistant to metal centers.¹² The role of Arg has been demonstrated in carboxypeptidase A, alkaline phosphatase, and liver alcohol dehydrogenase. We will see later in the section on protein tyrosine phosphatases that Arg can promote catalytic activity independent of metals as well when assisted by Asp, Cys, and other amino acid residues acting in concert. Carboxylates and coordinated water or hydroxyl can also show considerable geometric flexibility.¹³ Cysteine thiolate is a soft ligand due to the sulfur, which exhibits strong covalency and charge transfer particularly with Fe.^7,12 Cysteine can be in the anionic or neutral (protonated) form. Like tyrosine, Cys can form radicals in some cases.^14,15 An unusual example of versatility is observed in the enzymes galactose oxidase and cytochrome c oxidase, where tyrosine side chains at the active site are modified by covalent linkage to the cysteine and histidine side chains, respectively.^15a-d This covalent linkage, in conjunction with coordination to the metal ion (Cu(II)), may modulate the energetics of tyrosine radical formation during the catalytic cycle.^15a-d

Our goal in this review is to develop some general themes for catalytic reaction pathways in enzymes and synthetic systems via several examples, relating electron transfer, proton transfer, and charge flow to energetics and structural transformations. Most of the examples will be drawn from reactions of metalloenzymes, but we will also examine a single nonmetal containing enzyme (protein tyrosine phosphatase), a hydrolytic ribozyme (hammerhead ribozyme where the catalytic center is an aquated Mg(II) ion), and Lewis acid- and Bronsted acid-catalyzed synthetic reactions utilizing click chemistry. In both enzymatic and effective synthetic catalysts, the match between substrate and catalyst leads to both control and reactivity.

In click chemistry, the substrates are selectively reactive because they possess particular common functional groups, so that certain metal ions in solution (and other Lewis acids) are effective catalysts. These multiply bonded substrates are largely unreactive in most biological systems unless the protein cavity has a shape that is highly complementary to the assembled substrates. Because many derivatives of the two combining reactants can be prepared synthetically, both appropriate reactivity and high selectivity become feasible. Both in vivo and in vitro applications have been realized.

There are a number of general themes in this review: (1) Where and how does proton transfer facilitate bond cleavage, and how is this linked to charge transfer? (2) How are electron transfer and proton transfer coupled or gated? (3) How do the catalytic metal sites and amino acid residues provide orientation and guidance for the substrate along reaction pathways? (4) What are the comparative roles of active site flexibility versus rigidity, and how do these and related factors affect mechanism? When is flexibility important, and when is it counterproductive? What are the roles of substrate and protein strain? (5) What are the important protein environmental influences from the first coordination shell, second (or third) coordination shell, and the more extended protein/solvent environment? What are the "sensitive features" of the coordination environment that can account for effective or ineffective catalysis?

We will analyze reaction mechanisms in some detail for a number of illustrative examples; however, we will not be comprehensive in covering metalloenzyme catalysis because there are excellent reviews on related subjects, particularly on radical enzymes.^15d-f Our group has also written a number of recent reviews covering the electronic structure, properties, and energetics of iron-sulfur proteins,^16,17 structural, energetics, and reaction path issues for a number of metalloenzymes,¹⁸ and applications of density functional theory.¹⁹ In the specific areas of this review, we will focus on newer developments in these fields while trying to maintain coherence with important earlier work and fundamental concepts. Also, for many of the systems we will discuss, the catalytic reaction pathways are only partially understood, and certain steps are missing or possess uncertainties. There is considerable scope for further theoretical and experimental work on all of the systems analyzed here.

The main subjects are: (1) protein tyrosine phosphatases, nonmetal containing enzymes which hydrolyze tyrosine phosphate monoesters; (2) hammerhead ribozyme, which hydrolytically cleaves a phosphodiester backbone in RNA; (3) click chemistry, where Zn, Lewis acids, or Bronsted acids catalyze ring formation from multiply bonded precursors (azides, nitriles, and acetylenes); (4) MnSOD, FeSOD, and CuZnSOD which dismute superoxide to molecular oxygen and hydrogen peroxide; (5) the nature of the catalytically active site and the reaction pathways of the nitrogenase MoFe protein, where dinitrogen is reduced catalytically to ammonia; (6) intermediates of iron-oxo dimer enzymes (MMO, RNR); and (7) the reaction pathways of lipoxygenase, which are a class of Fe-containing enzymes that catalyze the hydroperoxidation of fatty acids.

The principal computational tools used here are density functional theory (DFT) methods either with the hybrid (B3LYP) potential²⁰ or with "pure" DFT (typically these are of generalized gradient approximation (GGA) type) exchange-correlation potentials, particularly Becke-Perdew86 (BP86)^21,22 or Perdew-Wang91 (PW91).²³ We have also used the more rapid local VWN-Stoll²⁴ exchange-correlation potential productively in the PTPase work, calibrating reaction energies there versus B3LYP, PW91, and MP2.²⁵

Hybrid potentials (of which B3LYP is one example) are derived from a total energy equation where the exchange-correlation energy is a linear combination of the Hartree-Fock exchange energy and "pure" DFT-GGA exchange and correlation energies.^19-26 The coefficients of these terms are typically determined from fits to thermochemical data. Application of the variational principle yields a set of mean field one-electron equations giving finally one-electron energies and orbitals as well as the corresponding self-consistent-field (SCF) potential(s). These SCF potentials are orbital dependent for hybrid potentials, while for "pure" DFT-GGA methods, all one-electron orbitals with the same spin share the same SCF potential. However, the SCF potential for spin electrons can be different from the SCF potential for spin electrons. For a discussion of both the fundamental and the practical aspects of using DFT methods, we refer to a very extensive literature.^19,26

We want to concentrate on a few significant points. DFT methods can be used for calculations on large systems with good accuracy for structures, properties, and energetics.^19,27 This applies even to transition metal complexes,^18,19,27-29 and the typical quantitative accuracy of these methods is often quite good for identifying and systematizing important energetic features of metal-containing enzymes, and for distinguishing feasible versus unlikely reaction pathways. A rough rule is to expect reaction energies and barriers intrinsically accurate to about 2-5 kcal/mol for B3LYP and somewhat larger for GGA-DFT (BP86, BLYP, PW91, BPW91).^19,28 Reaction energies and barriers for transition metal complexes typically involve bond making or breaking of only a few bonds, or reorganization of noncovalent interactions, which can also be reasonably strong. Typical barrier or reaction energies are often 0-50 kcal/mol; these are much less than total atomization energies, or dissociation energies for the complex as a whole, and so are the expected errors. Very often, the construction of the quantum mechanical model cluster (i.e., deciding what to include or exclude from the cluster), and finding a good representation of the surrounding protein and solvent environment, or solvent and counterions for synthetic systems, is equally important. Further, for transition metal complexes, the basis set must be of high quality, at least double or triple plus polarization on the metal, and double on the ligands.³⁰

Many transition metal complexes have high-spin (HS) metal sites, and spin-polarized methods are required with different SCF potentials for different spin. These distinguish between the orbitals and corresponding electron densities of the (spin-up) versus (spin-down) electrons, and the one-electron energies _i for the orbitals are also spin dependent.¹⁶ These electronic spin polarization effects can be very large when there are high-spin metal sites (as in many Fe and Mn enzymes) and can profoundly change the energy level scheme as confirmed by photoelectron and optical spectroscopies.³¹ The versus densities may be different either locally or globally. For high-spin states of the molecule as a whole, a single simple spin unrestricted DFT calculation may suffice: this covers a single spin-polarized metal site, or a number of parallel spin sites including spin delocalization to ligands. However, dinuclear and polynuclear transition metal complexes are often antiferromagnetically (AF) coupled, and metal-ligand radical complexes can also be spin-coupled.^16,18 In this case, the AF-coupled state is represented by a "broken symmetry" (BS) state, where the and electron densities occupy different regions of space, and the spin sites are oppositely aligned. The energy and properties of the BS state can then be compared to the corresponding HS state where the spin vectors of the sites are parallel aligned, also called the F (ferromagnetically coupled) state. For a Heisenberg exchange coupling Hamiltonian of the form, H_spin = -2JS₁^.S₂, the energy difference between the HS and BS states is

All of the results above apply in the "weak coupling regime" which is typical for dinuclear and polynuclear F- or AF-coupled complexes. Further, the results above are readily generalized to polynuclear spin-coupled complexes; in that case, there may be several different BS states to consider, representing different spin alignments, and replacing a far larger number of pure spin states of the spin-coupled system.

We note also that the individual metal site spins can also be variable, ranging from HS to intermediate-spin (IS) or low-spin (LS), depending on the ligand field and other electronic influences, which changes the quantum numbers S₁, S₂. The total system energy difference between high-spin and broken symmetry states (or between the high-spin and "true" AF-coupled E(S_min)) can be a chemically significant quantity, of order 2-16 kcal/mol for oxo, hydroxo, or water bridged Fe(III)-O-Fe(III), or Fe(IV)-O-Fe(IV), and this can have important implications for reaction pathways.

The phenomenon of valence delocalization introduces additional complexity in mixed-valence dimers, or for mixed-valence pairs of sites in polynuclear complexes, when the site geometries are similar and the metal site energies are nearly equal.²⁷ For each of the states of the Heisenberg spin ladder, there is then an additional resonance delocalization splitting term, which is linearly dependent on the total dimer spin (S₁₂) of the mixed-valence pair, as B(S₁₂ + ¹/₂). The energy of bonding or antibonding electron delocalization is dependent on the spin alignment of the mixed-valence pair. This phenomena, called "resonance delocalization coupling" or "double exchange" or "spin-dependent-electron-delocalization", occurs often in polynuclear complexes, particularly of 3Fe4S and 4Fe4S type where there are mixed-valence pairs of Fe sites. It is more rare in dimers, but a number of cases have been spectroscopically confirmed. For the large MoFe₇S₉X supercluster in the nitrogenase catalytic site, the current evidence supports some delocalization over multiple sites, but probably less than in 4Fe4S complexes. Additional studies are needed here in view of the newly discovered central ligand X. The identity of ligand X and the oxidation states and reactions of the MoFe superclusters will be analyzed later.

Active site models include at least the first coordination sphere for the metal site and often considerably more within the "quantum cluster" region. The problem of full geometry optimization is briefly analyzed below. Performing partial geometry optimization is also quite feasible by related methods, either with some atoms frozen in Cartesian space, or by converting to internal coordinates, and then freezing appropriate internal degrees of freedom. There can be a number of reasons for doing this. Certainly, the size of the geometry optimization problem is considerably reduced; this saves computer time, and one may know that certain degrees of freedom will not change much. Another reason is that one may want to embed the quantum cluster in a more extended protein (or synthetic) environment, and so freezing certain atoms in a frame, or utilizing internal coordinate constraints, may be needed to be compatible with the surrounding protein atom framework (known experimentally, for example). These mixed constrained/free clusters are used frequently, and geometry optimization is entirely feasible.³² While the forces in the unconstrained directions are zero at extrema, the forces in the constrained directions are ill-defined, and therefore obtaining frequencies, and clearly identifying minima versus transition states, is not feasible by frequency calculations. Another way of getting at transition states and reaction paths is to use linear transit methods between minima to obtain approximate saddle points, and then to refine these where needed by detailed transition state searching (see below).³³ Alternatively, one can start by guessing the position of transition states lying between selected intermediates. This guess is often based on the extended Hammond postulate which says that the transition state of a reaction is closer in the position of its reaction coordinate to the (initial or final or intermediate) state having the higher energy.³⁴ The Hessian matrix is defined as the matrix of second derivatives of the total energy with respect to the geometric coordinates and is highly geometry dependent. Because transition state searches require a good quality initial Hessian matrix, and are finally successful only when one finds the "right region" in the multidimensional conformational space, linear or quadratic transit methods are very useful. Also, all minima (or extrema, including saddle points, transition states) are determined by a "local" variational equation to obtain a zero value for the gradient vector (the gradient of the total energy with respect to geometry), and extrapolation is required, so searching by guess or transit paths is required in practice. The "local model" is usually a multidimensional quadratic function and cannot represent the true complexity of the potential energy (PE) surface. Further, the PE surface is typically more complicated near transition states than near intermediates.

The use of prior structural information is certainly valuable and often essential, where at least the "resting state" structure is available from medium or high-resolution X-ray crystallography.³⁵ In some cases, X-ray structures of intermediates or states complexed to inhibitors are available.³⁶ Also, spectroscopic evidence may provide valuable (partial) structural information where X-ray structures are absent and can complement and (even correct) X-ray structures in many cases.³⁷ All of these methods relieve the working scientist from the intractability of global searches for the lowest of multiple minima or multiple transition states in large molecules. Even a medium-sized molecule can have over 100 conformational minima connected by several hundred transition states.³⁸ Having protein structures of good quality (about 2 Å resolution) available solves for the minima of the "soft modes" in the protein conformational space and provides the extended protein/solvent framework in which to embed the active site transition metal complex or "quantum cluster". Still it should be remembered that X-ray structures can have both some heterogeneity of structure as well as structural uncertainties.³⁹ The structure may correspond to a mixture of redox states, a mixture of protonation states, or a mixture of tautomeric or conformational states. There may be uncertainties about redox and/or protonation states as well.^18,19,37

The question of whether the X-ray structure that is observed from the crystallized protein is the relevant "catalytically active" form of the protein is highly complex and very system dependent. We cite a few relevant examples to illustrate the possibilities. MnSOD and FeSOD are found either in dimeric or in tetrameric form.⁴⁰ The dimeric interface between the A and B chains in these proteins contains a Glu-His hydrogen bond where the His is a ligand to the Mn or Fe. This appears to be very significant catalytically based on DFT/electrostatic redox potential calculations, and in accord with the common appearance of Glu-His-Metal and Asp-His-Metal motifs in metalloproteins.^41,42 Further, this dimer interface provides access to the O₂^- substrate to get near the active site pocket. In other cases, the X-ray crystal structure may display highly interesting, although nonphysiological, interactions, even fairly close to the active site transition metal complex. A very recent example is a very-high resolution (1.3 Å) structure of the water-soluble Rieske iron-sulfur protein fragment from a thermophilic eubacterium expressed in E. coli.⁴³ The A and B molecules, each containing an (Cys)₂FeS₂Fe(His)₂ cluster, display H-bonding of shared protons between the His ligands from each molecule. By contrast, the full cytochrome bc₁ complex is still a dimer, but it is enormous.⁴⁴ Each half contains two heme b centers and a heme c₁ center in addition to one Rieske protein. Here, the Rieske iron-sulfur centers do not contact one another during the physiological redox and proton pumping cycle. Instead, they H-bond to hemes and/or to ubiquinones, and the Rieske protein fragments are evidently highly mobile during the catalytic cycle. The full complex is embedded in a membrane, but with a water-soluble part as well, and the resolution of the relevant X-ray structure is not as high (3.5 Å). A third example is the X-ray structure of the hydroxylase component of methane monooxygenase. Here, the available protein structures are of medium resolution (about 2 Å).^45,46 This resolution is good enough to obtain a general layout of the active site diiron cluster with respect to the remaining protein, but not sufficient to define the protonation state(s) of the bridging solvent-based (water or hydroxyl or oxo) ligands. Here, optical spectroscopy and DFT calculations help to define the active site coordination environment.^47,48 A general problem, cited in all of this literature, is the movement of protein side chains, which may or may not have a unique conformation within an X-ray structure. Here again, spectroscopy and quantum chemical calculations may help to sort out the problem. Also, solvent molecules may enter or leave the active site (or change protonation states), and these can play a mechanistic role. While medium-resolution X-ray structures (about 2 Å) are often adequate for constructing reasonable active site models and including the remaining protein environment, very-high-resolution structures (1.3 Å or lower) are highly desirable, particularly if the state in the X-ray structure can be closely monitored by spectroscopy.

Another point to remember in some mixed quantum/classical electrostatics models is that when the model of the surrounding environment allows for the evaluation of interaction energies but does not include analytic gradients, frequency calculations and precise transition state searches cannot be done. However, approximate transition states can be found by linear transit and related two-dimensional grid searching, and these can be quite valuable (see the discussion of Poisson-Boltzmann models in section 2.4).

If the quantum cluster is not too large, it is often feasible to perform full geometry optimizations on the ground state, intermediates, and transition states.⁴⁹ One finds where the gradients of the potential energy surface are zero (that is, the forces are zero). Most DFT programs can calculate energy gradients analytically, but analytic second derivatives are not available. Instead, an approximate inverse of the Hessian matrix is constructed and updated, particularly by the BFGS (Broyden-Fletcher-Goldfarb Shanno) method.⁵⁰ In the Newton method, the quadratic function approximating the potential energy surface is used to extrapolate the point where the gradient vector will be zero. Because the quadratic model surface is only an approximation to the actual surface, and because the Newton method requires analytic second derivatives, more flexible quasi-Newton methods are typically used.⁵⁰ In quasi-Newton methods, the energy minimization is performed along a line between the current and previous point, and after the approximate inverse Hessian is updated, the inverse Hessian matrix and gradient vector are used to take the following step, possibly with a smaller step size than in the Newton method. If convergence in the gradient and geometric displacement is not obtained, this process is performed iteratively.

Frequency calculations can be used to confirm that the ground state and intermediates are true minima within the framework of the Born-Oppenheimer surface of energy E({R_j}) versus geometry {R_j} (where {R_j} denotes the set of nuclear position coordinates j = 1,....,M for M nuclei), and because the Hessian (converted to the equivalent mass-weighted coordinate force constant matrix F) is positive definite for minima (after diagonalization), all eigenvalues from the equation⁵¹

There are several ways to represent the protein and solvent environment, which differ in assumptions, accuracy, and ease or difficulty of computation. The models and related problems include the following:

(1) Gas-Phase Calculations with or without Constraints. A "large enough" quantum model system is then needed. The choice of a physically realistic quantum model system becomes more difficult if the active site quantum cluster has a net charge or is quite polar (for example, having charged parts).

(2) Stability and Convergence. Many active site complexes in electron-transfer proteins and enzymes are multiply charged anions. Dinuclear and polynuclear iron-sulfur clusters are the clearest and most dramatic examples, because these are often multiply charged (Fe₂S₂(SR)₄^2-,3-, Fe₄S₄(SR)₄^2-,3- and simpler analogues Fe₄S₄(X)₄^2-, X = Cl^-, Br^-).⁷ In this respect, they (and related synthetic complexes) share common features with "textbook" dianions and trianions CO₃^2-, SO₄^2-, PO₄^3-.⁵⁴ While all of these systems are stable in condensed phases, they are unlikely to be stable in the gas phase because of the strong electrostatic repulsion of the excess electrons. The strength of this repulsion from classical electrostatics is roughly proportional to Q²/2R, where Q is the total cluster charge and R is the average radius, and this simple electrostatic contribution to the ionization potential (IP) in the gas phase becomes more negative as (2Q + 1)/2R (where Q is a negative integer, representing the cluster charge prior to loss of an electron, and the numerator is in units of e² where e is the electron charge). Clearly, the cluster becomes more unstable to loss of an electron as the cluster becomes more negatively charged.⁵⁵ Such clusters in the gas phase are also often unstable to loss of a ligand, or in the case of Fe₄S₄ (X)₄^2-, X = Cl^-, Br^-, to symmetric fission into two equal dinuclear fragments, 2[Fe₂S₂X₂^-].⁵⁶ From the experimental photoelectron spectroscopy, the first IP of Fe₄S₄(X)₄^2-, X = Cl^-, Br^-, is near 1 eV (positive), so these states are stable. However, with less electronegative terminal ligands, X = SCH₃^-, our recent DFT calculations give an IP of nearly 0 eV (slightly negative), so this cluster is unstable (or possibly marginally stable). For the Fe₄S₄(X)₄^3-, X = SCH₃^-, cluster, we calculate IP = -4 eV, which is very unstable.⁵⁷

It is surprising then that some multiply charged anions can be observed even when their measured ionization potentials are negative.⁵⁸ This has been seen by photodetachment photoelectron spectroscopy. It reflects the metastable character of these states. While one-electron loss is energetically favorable, the excess electron must surmount a repulsive Coulomb barrier arising from the outer shell of anions. Wang et al.⁵⁸ have studied the behavior of the copper phthalocyanate tetrasulfonate tetraanion in some detail [Cu(II)Pc(SO₃)₄] ^4-; the electron binding energy is negative (-0.9 eV), while the repulsive Coulomb barrier of the outer (SO₃) ^- groups is positive (3.5 eV), so that the photoelectron kinetic energy exceeds the energy of the incoming photon by 0.9 eV. At photon energies below 3.5 eV, electron tunneling through the Coulomb barrier is seen. We raise these issues not only because of the increasing role of electrospray ionization of large systems which can be combined with photoelectron spectroscopy detection. These metastable states also have computational analogues as demonstrated by Stefanovich and co-workers.⁵⁹ They have shown that negative ionization potentials can be calculated for a number of small multiply charged anions in the gas phase even with IPs as negative as -5 or -6 eV. These states are greatly stabilized in solvent, giving positive IPs in water. We have made and reported similar calculations on dinuclear and polynuclear iron-sulfur clusters over many years.^57,60 What then accounts for the ability to calculate these states, despite their instability toward loss of an electron?⁶¹ The resolution to this paradox lies in the Coulomb repulsive barrier, and in the fact that the continuum orbital into which the excess electron must "leak" has to be represented by a set of very diffuse functions, equivalent to plane waves or augmented plane waves. Even quite good (but not perfect) basis sets cannot represent the continuum electron, and a metastable state confining the excess electron is constructed. The energy and charge distribution in this state are generally very reasonable, because after the solvation energy due to the solvent cage is included, the IPs in solvent are reasonably close to those seen experimentally.^57,59,60c Solvent, counterions, or a protein environment rich in hydrogen bonds create a stabilizing environment while not changing the electron density in a major way.

Multiply charged anions do not, however, constitute the only stability problem for transition metal active sites. DFT calculations on active site cluster models, whether having a net positive or net negative charge, are susceptible either to oscillations and/or to giving an incorrect ground-state electronic structure in many cases unless care is taken and the physical origin of the problem identified. We will discuss two frequent problems that we have found in our recent calculations on FeSOD,⁴² the 2Fe2S Rieske center,⁶² and the FeMo cofactor center of nitrogenase.⁶³ During the self-consistent-field (SCF) process, it is easy for anionic side chains, particularly carboxylates and thiolates, to lose electrons to the metal ion (for Fe(III), in particular) and produce ligand radical states. These can appear to be the ground state even when they are not according to a proper test of total energies of alternative states. This problem often occurs when the cluster model does not contain the H-bonds to the side chains stabilizing the anionic form of the ligands. The cluster model can then be expanded to include the required H-bonds. Alternatively, level-shifting and orbital following procedures can be used to force the "hopefully correct" electronic structure for the SCF ground state. After one has found an SCF solution for one electronic state (corresponding to one specific orbital occupation scheme), the total energy of this state can be compared to those for alternative electronic states (corresponding to different orbital occupancies), and then the "correct ground state" can be determined at least for the specified exchange-correlation potential used. Orbital level shifting is usually turned off as self-consistency is approached to test the stability of the "trial electronic ground state" as compared to others. However, this procedure sometimes leads to renewed oscillations between electronic states. The "correct ground state" may well have some empty one-electron energy levels which lie below the filled energies, giving a "non-Fermi" occupation scheme. It becomes important to recognize when this is occurring, and how to deal with this problem in general.

This second problem of oscillations between filled and empty levels has its origin in the way energy levels are ordinarily filled within the SCF scheme in density functional calculations. It is typical and convenient to use the "aufbau" or building principle, which implies that the one-electron energy levels (derived from the solutions of the Kohn-Sham equations for a given exchange-correlation potential) are ordered by increasing energy, and filled in that order with integer occupation numbers (1 or 0 for occupied or empty spin-orbitals). However, the difference in orbital energies between an occupied and an empty orbital _j (n_j = 0) - _i (n_i = 1) (where n_i, n_j are occupation numbers for the active initial and final orbitals (i,j) of the excitation) can be a poor approximation to the one-electron excitation energy between states if the orbitals have very different spatial extent. A very good approximation to the total energy difference for a one-electron excitation is the Slater transition state energy,⁶⁴ defined as _j (n_j = ¹/₂) - _i (n_i = ¹/₂), with all other occupation numbers unchanged. If two energy levels are close and/or oscillating with changes in occupation numbers, the Slater transition state energy can be computed. This value will report which electronic state has the lower energy with good reliability and is stable against state-to-state oscillation because both orbitals have equal occupation numbers. The major case where this problem occurs is when an empty metal d orbital is close in energy to a filled ligand orbital (or conversely). The _j for that metal level (being the more localized) then drops substantially when the orbital is empty and rises substantially when it is filled, while the more delocalized ligand level _i (n_i) changes much less. The Slater transition state is also a good "halfway point" for trying to converge and compare different SCF states. It can then be followed by level shifting with more confidence that the "correct levels" are being shifted away. Our group prefers this approach to "orbital smearing" using variable occupation numbers, where the physical meaning is much less clear despite its utility.^49b State crossings or near avoided crossings between different electronic states can occur with variations in molecular geometry as well.⁶⁵ For example, a "ligand radical state" with a reduced metal ion may lie close to the "ligand anion state" with an oxidized metal ion.⁶⁶ In "valence tautomeric complexes", these distinctive electronic ground states may each occur under different conditions of temperature and pressure due to their differing entropies and volumes. A proper physical analysis of these problems should include an assessment of spectroscopic properties as well as geometries and energies because the electronic states can be very close to one another.²⁹

(3) Inclusion of an "Average" Polar Environment Using a Continuum Solvation Model, Treating Interaction Energies Either at the Single Point Quantum Level (after the "Gas-Phase" Geometry Optimization on the Quantum Cluster) or with Subsequent Geometry Optimization Including Solvation. One of the most popular models of this type is COSMO (conductor like screening model for two dielectric media, quantum cluster plus solvent),⁶⁷ which has good accuracy for polar solvent media, when the quantum cluster charge and polarity is not too high. The protein and solvent are not explicitly represented; rather an average dielectric constant is chosen for the dielectric medium. COSMO starts from a solvent model with a dielectric constant equal to infinity, and then rescales this back to a finite dielectric using a well-known rescaling formula: f() = ( - 1)/( + x). The variable x is adjustable: for a finite charge in a spherical cavity, x = 0, while for a dipole in a cavity, x = ¹/₂. For large , x becomes much less important, and COSMO is most accurate for large in any event. One great advantage of COSMO is that geometry optimization can be performed with COSMO as a part of the potential, and this option is available in a number of computer programs, including Gaussian and ADF.⁴⁹

(4) Poisson-Boltzmann-Based Models. These models are used to represent the energetic interaction of the active site with the protein and solvent. They can be used as "classical electrostatic models",⁶⁸ but we will discuss their use in combination with quantum models for the active site.^25,60b

The goal is to derive potentials at the active site and corresponding interaction energies with the active site cluster by solving the linearized Poisson-Boltzmann equation, which takes the form:

For the combined DFT-PB method, the simplest model is the single step rigid charge model.¹⁷ Here, all of the generating charges are point charges (either cluster ESP charges or protein charges) and the total interaction energy is calculated as the sum of the reaction field and the protein field energies E_pr = E_react + E_prot.

In the more complete self-consistent-reaction-field (SCRF) method, the full electron density distribution is used in quadratures, and the reaction-field potential produced by the cluster charges is iterated to self-consistency along with the quantum electron density (and its ESP charge representation).^25,52,60b,70 The total interaction energy E_pr can be partitioned into three terms:

One major advantage of Poisson-Boltzmann (PB)-based methods is that a statistical mechanical evaluation of the protonation state of the various titrating residues of the protein versus pH can be carried out prior to a quantum mechanical calculation. Monte Carlo and related methods can be used to reduce the (2)^N scaling of this problem when there are many (N) titrating sites.⁷¹ Another advantage, which also carries over into the combined QM-Poisson-Boltzmann methods, is that the protein field interaction energy with the quantum cluster can be partitioned into terms from the individual protein residues (or into main-chain versus side-chain terms) when the linearized PB equation is used. This greatly facilitates energy analysis and has been used extensively for designing calculations with larger quantum clusters from prior calculations with smaller QM regions plus an electrostatic PB representation of the remaining protein/solvent. The QM-electrostatic region boundary is typically represented using a link atom approach, similar to that often used in QM-molecular mechanics (see below). For both QM-PB and QM-MM, charge transfer across the boundary is not properly treated because this is intrinsically quantum mechanical in nature; Pauli repulsion between QM and PB or MM atoms is either omitted or treated with a simplified classical force field approach. However, combined QM-PB or QM-MM approaches are clearly more rapid than the equivalent fully quantum mechanical representation of the problem. For QM-PB, the dielectric constant for the protein region can represent both electronic polarization and orientational polarization, which roughly reflects protein dynamics, rather than a single structural minimum.

(5) QM-MM Method. An alternative approach to represent the protein and solvent environment is the hybrid quantum mechanical (QM)-molecular mechanical (MM) method.⁷² Just like the continuum dielectric method described previously, this method partitions the whole system into two regions: an active region and the rest of the system. The active region includes a chemically active part such as the catalytic site of the enzyme and is treated with quantum mechanics. The rest is treated with an empirical force field. The whole system is described by a mixed Hamiltonian:⁷³

One advantage of the QM-MM approach is that the computational effort can be focused on the active region of the system where quantum phenomena such as bond-breaking/bond-forming and electron/proton transfer occur, whereas the effects of surrounding environments are taken into account by fast but less accurate approaches. In this context, the geometry constraint from the protein environment to the active site cluster is automatically implemented during the geometry optimization step. The accuracy of the QM-MM calculation relies not only on the level of the quantum mechanical Hamiltonian, but also on the quality of the classical force field. Two of the shortcomings in many classical force fields are the point charge model for the electrostatic potential and the lack of polarization. Development of a polarizable force field⁸¹ and integration of the electrostatic potential with a Poisson-Boltzmann (PB) or General Born (GB) model⁸² are still the subjects of ongoing investigation. Nonetheless, the QM-MM method has been widely applied to study the solvation effect on protein structures⁸³ and the mechanisms of many catalytic reactions in solution phase or in protein environments,⁸⁴ and the results have been reviewed elsewhere.⁸⁵

(6) Car-Parrinello Molecular Dynamics. The fluctuation of protein conformation and side-chain orientation has profound effects on the thermodynamics and kinetics of enzyme catalytic reactions. Traditionally, simulation of protein dynamics has been dominated by the molecular dynamics (MD) method based on the classical force field, which was hard to extend to the quantum region. Parametrization according to quantum mechanical calculations will provide reasonable force field parameters for the system under study.⁸⁶ In the core of the MD method is the calculation of forces acting on the nuclei for the propagation of the nuclear trajectories. The calculation can be easily done with classical molecular mechanics, but is prohibitively expensive with quantum mechanics for the larger systems such as proteins. Car-Parrinello molecular dynamics (CPMD), however, opened an avenue to carry out quantum MD simulation, especially under the density functional theory framework.⁸⁷ This was achieved by a revolutionary approach that solves the electronic wave function and the forces acting on the nuclei simultaneously. The approach is based on the fact that the variational principle is simply a minimization procedure and the parameters in the electronic wave function may be treated as similar to the other dynamical variables nuclear positions, in the Lagrangian. The electronic structure problem and the dynamics of the atoms are hence included in a set of Newtonian equations of motion and can be solved concurrently with the steepest descent method. In this way, the full dynamic time evolution of a structure is computed without resorting to a predefined potential energy surface (PES). The PES is calculated on-the-fly as the nuclear trajectory is generated. The implementation of CPMD was empowered by the pseudopotential approach using plane wave basis sets.⁸⁸

CPMD has been applied to many systems including a wide range of biological systems such as metalloenzymes where the metal-protein interactions are hard to capture by a classical force field and the catalytic procedures involve bond-breaking and bond-forming, as well as proton and electron transfer. These applications have been reviewed thoroughly in several recent reviews.⁸⁹ Combining CPMD with QM-MM methods extends the horizon of CPMD simulation to more realistic complex systems.⁹⁰ Recent work from QM-MM CPMD simulations includes identification of the likely location of Cu²⁺ binding in proteins and the mechanism of cis-trans photoisomerization of the covalently bound chromophore in rhodopsin. The model system in the later study consists of chromophore, protein, and membrane mimetic environment, and the whole system is 24 000 atoms.^89a However, many challenges still have to be met before the CPMD method can become as popular and powerful as the classical MD method.

Protein tyrosine phosphatases (PTPases) are enzymes central in regulating many cellular processes, particularly in response to extacellular signals.⁹¹ PTPases are complementary in their function to protein tyrosine kinases, because in PTPases, a phosphate group is hydrolyzed catalytically from a tyrosine side chain, while in kinases, phosphate is added to the side chain.⁹² A proper level and timing of tyrosine phosphorylation is critical for regulating cell growth, differentiation, metabolism, and progression through the cell cycle, as well as cell-to-cell communication and cell death versus survival.^93,94

PTPases constitute a large and diverse superfamily of enzymes, which include both cytosolic and membrane-bound receptor enzymes.^93,95 The different subfamilies are diverse in sequence, molecular weight, and specificity (dual specificity enzymes can hydrolyze phosphoserine and phosphothreonine as well as phosphotyrosine), but there is considerable experimental evidence that PTPases exhibit a common mechanism with related common structural features.

Although these enzymes operate without metal ion cofactors, tyrosine phosphatases are extremely efficient enzymes, showing maximum catalytic rate constants (k_cat) 10 orders of magnitude greater than the corresponding rate constants for uncatalyzed hydrolysis of phosphate monoesters.⁸¹ While there is general agreement about the mechanism of the uncatalyzed reactions, there is considerable dispute about the mechanism of the enzymatic reaction. We were drawn to the study of PTPases both because of their biological importance and because the proposed catalytic mechanisms based on theoretical and experimental studies displayed sharp contrasts.^25,96 Central to this dispute is the apparent high charge of the free substrate (a dianion) and the question of whether or how this charge persists (or is altered) first in the Michaelis complex, and then in forming the first stable intermediate. It is known that the reaction proceeds through a Tyr-phosphate-cysteine intermediate with a direct P-S bond. If the cysteine is a 1^- anion,⁹³ the reaction takes the form:

As we shall see, the difference in charge between Mechanism A and Mechanisms B and B' leads to very different reaction barriers and pathways. In the following analysis, we will focus on the formation of the covalent phosphocysteine intermediate rather than the final hydrolysis step. Depending on the specific enzyme studied, either of these can be the rate-limiting step, and both are biochemically significant, with the formation of the first intermediate constituting the first irreversible step of the reaction.

In a series of mechanistic studies, Aquist et al.⁹⁷ used the empirical valence bond method to examine the reaction mechanism of PTPases and concluded that mechanisms B and B' are energetically most favorable. Mechanism B (or B') does seem intuitively reasonable because the nucleophilic attack of the CysS^- anion on the phosphorus of the substrate should be easier if the substrate were a monoanion, as in Mechanism B, or after the activation step in Mechanism B', as compared to the dianionic substrate of Mechanism A, where greater charge-charge repulsion would be expected. However, eqs 9 and 10 or eqs 12 and 13 alone do not uncover the actual sequence and timing of the proton transfer from the general acid to the tyrosine anion, or the role of other protein residues at or near the active site.

Gao and co-workers⁹⁸ used an AM1/MNDO semiempirical quantum mechanics treatment of thiolate, substrate, and general acid in the context of the protein environment to explore the reaction pathway, concluding in favor of mechanism A. The protein environment was described by molecular mechanics using the PARAM22 force field in CHARMM, so there was a significant representation of the surroundings. However, this group did not predict the protonation of the tyrosine anion that is needed for completion of the first half of the reaction.

Czyryca and Hengge⁹⁹ used a large quantum cluster of nearly 300 atoms to model the reaction pathway assuming Mechanism A with a PM3 semiempirical Hamiltonian. Their focus was on the structural transformation along the reaction pathway. They predicted a dissociative transition state for the formation of the intermediate. The phosphoryl transfer was predicted to be concerted, but with a "very loose" transition state, and where the bond-breaking and bond-forming processes occur concurrently. Even with such a large cluster model, extensive electrostatic interactions beyond the quantum cluster and long-range solvation effects are very significant due to the large charge on the substrate and charges also on a number of protein residues. Further, a number of protein dipolar groups are aligned near the substrate reaction site to stabilize the substrate and thiolate charges. There are two transition states seen using coordinate driving by incrementally stretching the P-O bond and optimizing the remainder of the structure. These transition states occur at P-O bond lengths of 2.2 and 2.6 Å; for the latter, the calculated P-S bond length is 3.3 Å, which is a reasonable geometry for a very "loose" transition state. While geometrically sensible, the predicted transition state energy barrier is very low, about 3 kcal/mol, which is much smaller than that experimentally observed in a different PTPase (about 14 kcal/mol).¹⁰⁰

Despite their extensive structural diversity, all PTPases share a common active site sequence, consisting of (H/V)CX₅R(S/T) residues, which is called the phosphate binding loop or "P-loop". This loop contains the invariant cysteine nucleophile, a serine or theonine that H-bonds to it, and an arginine nearby. In another region where there is little sequence similarity among PTPases, there is a general acid, which is always Asp, that resides in a similar position with respect to the substrate binding site in all cases.

Dillet et al.⁹⁶ have used Poisson-Boltzmann electrostatics calculations to examine the expected protonation states of the titratable residues in a variety of PTPases. They have found that the nucleophilic Cys is always in the anionic form in both the free unliganded enzyme and the Michaelis complex. The pK_a of the general acid Asp is always shifted up from its reference value in the unliganded enzyme and becomes even more basic in the Michaelis complex. (That is, the Asp tends to be protonated even at higher pH, a typical range being 5.4-9.1.) Figure 1 shows schematically the expected form of the Michaelis complex (adapted from Figure 3 of Zhang⁹⁴) showing how the dipoles of the P-loop can stabilize the bound substrate along with the arginine cation, and showing also the positions of the conserved cysteine and the aspartate general acid. Dillet et al. found that the dipoles of the P-loop and a number (2-7) of positively charged side chains from Arg, His, and Lys make major contributions in stabilizing the Cys anion state, counteracting both the phenyl phosphate substrate dianion and the cost of burying charge within the protein which would favor the neutral Cys form. However, these charges and dipoles act fairly locally so that the Asp general acid pK_a stays high. Experimental support for these results derives from the pH dependence of the k_cat/K_m ratio, which indicates apparent pK_a's for the free enzymes, and from the pH dependence of k_cat, which indicates the apparent pK_a's of the Michaelis complex for those enzymes where the formation of the phosphoenzyme intermediate is rate limiting.^93,94,101 In these cases, the concentration of the Michaelis complex will greatly exceed that of the phosphoenzyme intermediate. Additional support for the protonation state of Asp and its role as a general acid catalyst comes from mutational studies in a low molecular weight phosphatase.

Figure 1 The active site of protein tyrosine phosphatase with bound substrate (adapted from Figure 3 of ref 94).

In recent work, our group²⁵ used a combination of DFT methods (Vosko-Stoll exchange-correlation potential) along with a self-consistent-reaction-field approach for the extended protein-solvent environment to evaluate the comparative reaction pathways of Mechanism A versus B for the low molecular weight (bovine) phosphatase. In this way, the effect of the charges and dielectric media outside the quantum region can be included in the description of the reaction path energetics, and the quantum cluster can polarize in response to the protein and reaction field potentials. The large extended electrostatic effects indicated by Dillet's prior electrostatics calculations can be treated along with bond-breaking and bond-making processes. The quantum cluster used is depicted in Figure 2. The P-loop dipoles reside in the electrostatic "protein" region, while the Asp129, Arg18, Cys12, Ser19 side chains, and model phenyl phosphate substrate form the quantum cluster.

Figure 2 Schematic of the quantum cluster for the active site of bovine phosphatase. The P-S distance is the reaction coordinate.25

The energetic pathway for Mechanism A along the reaction coordinate, taken as the P-S distance as the intermediate is formed, is shown in Figure 3, and the Michaelis complex and the transition state structure are shown in Figure 4.

	Figure 3 Mechanism A energetics along the reaction coordinate. The transition state (TS) and Michaelis complex (MC) points are appropriately labeled. The forward direction of the reaction corresponds to the right-to-left direction along the plot.25
	Figure 4 Michaelis complex (top figure) and transition state (bottom figure) structures for mechanism A. For clarity, the Ser 19 residue has been deleted.25

The general acid Asp129 transfers a proton to the phenoxide group early, so that this process is nearly completed in the transition state. A "loose" transition state geometry is seen with a metaphosphate ion passing between the phenol and the Cys anion. This result agrees with the analysis of observed bridging and nonbridging O kinetic isotope effects, indicating a rather "free" metaphosphate ion, but without a stable intermediate for the reaction. With Mechanism A, the catalytic pathway has not changed the basic mechanism of phosphate monoester hydrolysis, but it has made it much more efficient. (We will also see below that Mechanism B is much poorer at lowering the transition state.)^92,99 The computed barrier height for Mechanism A is about 9.1 kcal/mol, which is somewhat lower than the experimental barrier of about 14 kcal/mol (see Figure 5). This agreement is reasonable given that zero-point vibrational energies and entropic terms were neglected, the size and complexity of this system, and that the experimental energy profile was obtained with p-nitrophenyl phosphate as a substrate, rather than phenyl phosphate. The role of the Arg cation is very interesting; it stays highly charged and shows some mobility, moving in concert with the PO₃^- group as the P-S bond is formed. One can think of the Arg cation as a "molecular guidewire", providing electrostatic guidance for the "free range of motion" of the PO₃^-. The arginine charge and mobility clearly facilitate the reaction, but charge transfer with the PO₃^- is limited. In fact, we may well have underestimated the mobility of the Arg because Arg was truncated at the N()-C() linkage and C() was replaced by a linking H. However, the resulting N-H bond was allowed to rotate about the H pivot, and the C-N-H bond angle was unconstrained. The Arg is also more highly charged (as is the PO₃^- group and the Cys thiolate) after the SCRF procedure than for the "bare cluster" (gas phase) calculation all the way along the reaction path, which is indicative of the positive potential exerted by the protein and reaction field, particularly from the P loop. (There are few charged residues near the active site in this low molecular weight phosphatase, in contrast to other PTPases where several positively charged residues act in concert with the P loop.) The gas-phase cluster calculation yields essentially no barrier at all from the Michaelis complex to the transition state. Further, if the Poisson-Boltzmann energy quadrature is calculated using the gas-phase electronic density and added to the gas-phase energy, the barrier is quite low, only 4 kcal/mol. This result is similar to the low barrier seen by Czyryca and Hengge, who had a large cluster, but no extended protein environment.⁹⁹ The structural implications and energetics of Mechanism B are shown in Figures 6and 7. An associative mechanism is clearly indicated with a very late proton transfer to the leaving phenoxide and an associative distorted trigonal bipyramidal complex in the transition state. With the late proton transfer to the phenoxide, there is no advantage to less charged starting structure. Our predicted barrier height is 22 kcal/mol for Mechanism B, which is 8 kcal/mol higher than experiment and 13 kcal/mol higher than the mechanism A results. A similar energy profile is displayed by Mechanism A if the generalized acid Asp is removed, which is likely related to the late proton transfer in Mechanism B.

	Figure 5 Comparison of different models for environment effects in mechanism A. ·, full-SCRF; , SCRF scheme terminated after one cycle; , electron density as in the gas phase (initial cycle); bare dashed line, initial cycle with interaction calculated with point charges rather than quadrature over density.25
	Figure 6 Michaelis complex (top figure) and transition state (bottom figure) structures for mechanism B. The rest is as in Figure 4.25
	Figure 7 Energetics along the reaction coordinate for mechanism A in the absence of acid catalysis. For comparison, the profiles for mechanisms A and B are also shown. ·, mechanism A without acid catalysis; , mechanism A; and , mechanism B.25

Mechanism A (eqs 9,10) is a more accurate representation of the reaction pathway than is Mechanism B (eqs 12,13,14). The DFT energy barrier from the SCRF calculation is much lower for the modified Mechanism A than for Mechanism B. Consistently, the experimental kinetic isotope effects (KIEs) favor Mechanism A because the KIEs clearly indicate a dissociative transition state involving a metaphosphate (PO₃^-) as does the DFT/electrostatics calculated pathway. Further, the pH profile for formation of the Michaelis complex implies that the substrate is a dianion, as in Mechanism A, and not a monoanion as in Mechanism B or B'. An essential feature of Mechanism A is that the Michaelis complex depicted in eq 9 and the transition state (Figure 4) are not associative because the aspartic acid transfers the proton early to the tyrosine (leaving group). A pentacoordinate Tyr-phosphate thiolate complex never forms.

In summary, in PTPases, we can see how a sophisticated and highly charged and polar protein structure allows the first step of a very efficient phosphate hydrolysis to occur. Our studies of the final hydrolytic step of the phosphocysteine intermediate are also nearing completion. The Arg cation acts in some ways as a metal cation, but it receives a great deal of help from the P loop, the conserved Asp and Cys, and other groups. In the following section, we will see how a simple Mg(2+) cation with its hydration shell performs in catalyzing hydrolysis of phosphodiesters.

Since their discovery some 20 years ago, the study of ribozymes has become enormously active in the fields of molecular biology and biomedical science. Ribozymes are catalytic RNAs that can promote chemical reactions in the absence of proteins. These catalytic RNAs may have played a role in self-replication early in evolution. Some ribozymes are self-splicing introns. These will catalyze the cleavage and removal of the intervening sequence (the intron) from an RNA transcript, and then ligate the message from the flanking exons. This process is extremely important for the regulation and function of messenger and ribosomal RNAs in eukaryotic organisms.¹⁰² Ribozymes can also be engineered to cleave other target RNA molecules, and, consequently, they are now widely accepted as agents capable of inhibiting gene expression. They are therefore, naturally, also associated with very promising links to candidates in gene therapy.

The hammerhead ribozyme is a small RNA molecule that makes up the genome of several plant viruses, viroids, and satellite RNAs.¹⁰³ For the cleavage of the phosphodiester backbone, a divalent metal ion is required as a cofactor to promote activity under physiological conditions.¹⁰⁴ The metal ion typically utilized by the organisms to activate the phosphate group of the phosphodiester linkage leading to hydrolysis is Mg²⁺, although other metal ions have been shown to promote catalysis in vitro.¹⁰⁵

The mechanism of nucleophilic displacement at the phosphorus has been proposed to be of S_N2(P) type, with an in-line geometry, where the approach of the 2'-oxygen to the phosphorus is at a 180 angle relative to the leaving 5'-oxygen,¹⁰⁶ leading to inversion of configuration about the phosphorus (Scheme 1).¹⁰⁷

Scheme 1

It has been shown experimentally that a metal ion ligates to the pro-R_P oxygen of the -O-(PO₂^-)-O- group undergoing reaction, a feature crucial in the metal ion catalysis of phosphodiesters, as cancellation of the negative charge allows for nucleophilic addition to phosphorus (Scheme 2).^{104,107a,b,108} Furthermore, pH titration experiments indicate that a single deprotonation event is required for cleavage.^105c,109,110

Scheme 2

The metal ion can provide a ligated hydroxide to remove the proton from the 2'-hydroxide to generate the nucleophile (general base catalysis), or the nucleophile can be generated by a lyate HO^- (specific base catalysis) and the metal ion has an alternate role in the reaction. Of these metal ion binding sites, at least one appears necessary to achieve proper folding and to stabilize or position other essential functional groups, while one or more additional divalent metal ions participate in the cleavage reaction.^105d This divalent metal ion, or possibly another, can provide the necessary stabilization to the leaving group either by direct interaction with the leaving group oxygen or by orienting a water molecule such that a proton can be donated to this group. The way in which stabilization of the leaving group is achieved is a subject of debate,¹¹¹ and despite the numerous experiments performed on many different hammerhead ribozyme motifs, the precise details of the mechanism of hydrolysis and the number of metal ions essential to the cleavage reaction are not known.

There are a number of crystal structures of the hammerhead ribozyme that contain inhibitor, substrate, or modified substrate complexes,¹¹² and each has provided invaluable structural information and insight into the metal ion binding sites. None of these crystal structures, however, shows the 2'-hydroxyl group poised for an in-line displacement of the 5'-leaving group.

Several studies have used quantum chemical methods to investigate various aspects of the phosphodiesters hydrolysis and to model the mechanism in the hammerhead ribozyme.¹¹³ A transition state has been located depicting the P-O bond cleavage as the rate-determining step along the reaction coordinate, but this transition state was characterized by rotations of P-O bonds, an event that is highly improbable in the hammerhead ribozyme system. Lyne and Karplus suggested that the pK_a values of ionizable groups at the active site of the ribozyme might aid in indicating possible mechanistic pathways and in inferring the nature of the chemically active species.

Incorporation of explicit water molecules has resulted in an improved representation of the computational model in theoretical investigations of acid-catalyzed phosphodiester hydrolysis.¹¹⁴ Using Car-Parinello combined DFT-molecular dynamics calculations, Boero et al. examined reaction mechanisms for a model phosphodiester in a system which also contained 62 water molecules. For this ribose-phosphodiester in water, as a 1^- anion, the authors found a very high reaction barrier (~57 kcal/mol) when the system was driven through the P-O2' bond formation reaction coordinate. When Mg²⁺ was incorporated in the system as Mg²⁺(H₂O)₅, Boero et al. find that the Mg²⁺ ion is ineffective in phosphodiester hydrolysis, although the activation energy for deprotonating the O2' is lowered.

Very recently, DFT calculations were used to examine the effect of general base catalysis (Scheme 2) in the hammerhead ribozyme self-cleavage reaction.¹¹⁵ A neutral chemical model was used, consisting of a mixed phosphodiester of methanol and 3,4-dihydroxytetrahydrofuran (charge = -1) with a [Mg²⁺(HO^-)(H₂O)₄]¹⁺ bound to the pro-R_P oxygen, as shown in Scheme 3. The idea was to establish whether the presence of a single hydrated metal ion would be sufficient to enable catalysis and promote cleavage of the P-OCH₃ bond, the latter being commonly accepted as the rate-determining step in the reaction. This chemical model incorporates the basic features and necessary functional groups to adequately model the proposed hydrolysis reaction.^114,116

Scheme 3

A distinct, although short-lived, trigonal bipyramidal (TBP) pentacoordinated intermediate could be identified along the reaction coordinate. The transition state leading to the intermediate (TS1) was calculated to be +18.6 kcal/mol higher in energy relative to the reacting species and corresponds to the transfer of the proton from O2' to the Mg²⁺-bound OH^- and simultaneous nucleophilic attack of the O2' on the scissile phosphate. The furanose ring spans one apical and one equatorial position in accord with the proposed mechanism for hydrolysis in cyclic phosphate esters.¹¹⁷ In this asymmetric transition state structure, the O2'-P attack distance is 2.06 Å and the O2'-P-OCH₃ angle for nucleophilic displacement is 157.5. The phosphorus and equatorial oxygen atoms are nearly planar (+5.8 phosphorus deviation above plane) in contrast to that observed in the reactant (+24.0 deviation). Optimized geometries of the key structures along the reaction pathway are presented in Figure 8.

The energy of the TBP intermediate was found to be almost identical to the TS1 energy (18.6 kcal/mol above that of the reactant). It is characterized by an O2'-P distance of 1.98 Å, a P-OCH₃ distance of 1.80 Å, and an O4-P-OCH₃ angle of 157.4. The deviation of the phosphorus from the plane defined by the equatorial oxygen atoms is +4.0, indicating inversion of configuration about the phosphorus has not yet occurred. It is interesting here to note that a similar pentacoordinated intermediate, formed in the -phosphoglucomutase reactions, was identified by means of X-ray crystallographic techniques by Allen and co-workers.¹¹⁸

The intermediate is very unstable. A second transition state (TS2) corresponding to a P-OCH₃ bond breaking coupled to the transfer of a proton from Mg²⁺-bound water to the methanol leaving group is calculated to be only 2.2 kcal/mol higher than the intermediate energy (20.8 kcal/mol higher than the reactant). TS2 is an asymmetric transition state with an O2'-P distance of 1.84 Å, a P-OCH₃ bond distance of 1.96 Å, and an O2'-P-OCH₃ angle of displacement is 157.8, slightly larger than in TS1 and the intermediate structure. The O5'-P bond breaking appears to be late, while the first proton transfer is early. The final product is calculated to be 2.3 kcal/mol lower than the starting reactant structure. The structures of the reactant, TS1, the TBP intermediate, TS2, and the final product are shown in Figure 8. Relevant distances are indicated.

Figure 8 Geometries of key structures along the reaction pathway of the hammerhead ribozyme cleavage reaction.

The rate-determining step in the hammerhead reaction has been determined to be the departure of the 5'O^- leaving group.^105c,d,107a This step corresponds to the highest point (TS2) on the potential energy surface in our calculations. The calculated overall energy barrier of ca. 21 kcal/mol is consistent with the reported value of the rate of reaction of approximately 1 min^-1 (~20 kcal/mol) found in vitro for the hammerhead ribozyme.^119,120

The potential energy surface, displayed in Figure 9, suggests further that the lifetime of the intermediate structure is very limited, amounting to a kinetically insignificant^121,122 species along the reaction pathway. From TS1 to TS2, substantial changes in geometry occur with only a small change in energy (~2 kcal/mol). It can be seen that the structure of TS1 corresponds to the transition state for endocyclic cleavage, while the structure of TS2 corresponds to that for exocyclic cleavage. Therefore, this reaction does not appear to proceed by a simple concerted mechanism.

Figure 9 Potential energy curve for the hammerhead ribozyme reactions.

Metal-ligated hydroxide is a requirement to create the nucleophile for initial reaction to take place, while metal-liganded water provides the critical proton to the leaving group, consistent with experimental^{105a,b,107a,b,108} and molecular dynamics studies.^123,124 The energy associated with the highest lying transition state, TS2, shows that one metal ion is sufficient to catalyze the self-cleavage of the hammerhead ribozyme. The possibility of two metal ions participating in the cleavage reaction, however, cannot be eliminated.

Tetrazole chemistry is gaining increasing attention due to its wide-ranging applications. Tetrazoles are, for instance, used in pharmaceuticals as lipophilic spacers and carboxylic acid surrogates,^125,126 in specialty explosives,¹²⁶ photography, and information recording systems,¹²⁷ in addition to being precursors to a variety of nitrogen-containing heterocycles.¹²⁸

The most direct method to form tetrazoles is via the formal [2+3] cycloaddition of azides and nitriles. Depending on the nature of the azide species, evidence in the literature indicates that the mechanism of the reaction is different.

With organic azides, only certain highly activated nitriles are competent dipolarophiles.¹²⁹ In these cases, the reaction is regioselective, and only the 1-alkylated product is observed.¹³⁰ It is accepted that in these cases the reaction proceeds via a traditional [2+3] mechanism (see Scheme 4).

Scheme 4. [2+3] Cycloaddition of Nitriles and Azides

A mechanistically more interesting case is the addition of azide salts and nitriles to give 1H-tetrazoles. It has long been known that simple heating of certain azide salts with a nitrile in solution (typically 100-150 C) produces the corresponding tetrazole in high yield (see Scheme 5). This variant is much more synthetically useful, as the scope of nitriles that are competent reactants in this reaction is very broad, in contrast to the case of organic azides. In addition, a wide variety of metal-azide complexes are competent azide donors.

Scheme 5. Tetrazole Formation by Addition of Azide Salts to Nitriles

Mechanistically, these cases are considerably more complicated. Several possible reaction pathways can be envisioned. Claims have been made for both an anionic two-step mechanism¹³¹ and a concerted [2+3] cycloaddition,¹³² but the data are not conclusive. Moreover, there is evidence against both of these mechanisms. For example, while protic ammonium salts of azide are competent azide donors, tetraalkylammonium salts were not, refuting the strictly anionic two-step mechanism.^132a Also, while virtually all nitriles are engaged by ammonium azide salts at elevated temperature, organic azides only react with the most activated nitriles.¹³³ The fact that these azide salts and organic azides are electronically very similar, yet have significantly different reactivities, indicates that different mechanisms are likely in effect.

Density functional theory calculations using the hybrid functional B3LYP have been performed to study different mechanisms of tetrazole formation, including concerted cycloaddition and stepwise addition of neutral or anionic azide species.¹³⁴ In particular, the question of how the availability of a proton catalyzes this reaction was addressed.

It was shown that the activation barriers for the concerted cycloaddition with both the organic (neutral) and the anionic azide strongly depend on the nature of the substituent of the nitrile, that is, the level of activation of the nitrile (see Figure 10). For example, while CH₃CN has a calculated barrier of 31.6 kcal/mol for the [2+3] cycloaddition to CH₃N₃, the more activated CH₃SO₂CN has a significantly lower barrier of 20.4 kcal/mol.

	Figure 10 Calculated barriers for the [2+3] cycloaddition reaction of nitrile and azide as a function of nitrile substituent. Inserted are two typical transition state structures.
	Figure 11 Optimized structures of TS1 (A), intermediate P (B), and TS2 (C) of the mechanism displayed in Scheme 6. Here, R = Me and R' = H.

As mentioned above, only ammonium salts of azide that contain a proton are competent dipoles; tetrabutylammonium azide does not work, indicating a specific role for the proton. The calculations suggest that, when a proton is available, the reaction proceeds through a stepwise mechanism as shown in Scheme 6. In this mechanism, the azide anion and the proton source attack the nitrile simultaneously to yield a protonated intermediate (intermediate P). The activation of the nitrile by the proton facilitates the attack of the azide on the carbon of the nitrile. From intermediate P, simple 1,5 cyclization (TS2) occurs to give the 1H-tetrazole. In this way, the high barrier associated with the [2+3] cycloaddition is broken up into two steps, each with a smaller barrier, which is a common theme in catalytic processes.

Scheme 6. Suggested Mechanism for Tetrazole Formation by Addition of Azide Salts to Nitriles, Catalyzed by the Presence of Protons

In the calculations, the ammonium salt was modeled as a simple NH₄⁺ species. The optimized structures of TS1, intermediate P, and TS2 are displayed in Figure 11. TS1 was found to be ca. 21 kcal/mol high with acetonitrile as the dipolarophile, intermediate P is ca. 3 kcal/mol higher than the reacting species, while TS2 is only ca. 15 kcal/mol higher than the intermediate. The barriers of this pathway are to be compared to the barriers of the neutral cycloaddition (31.6 kcal/mol) or the anionic cycloaddition (33.8 kcal/mol), see Figure 12.

Figure 12 Comparison of the energies obtained for the various mechanisms for the formation of tetrazole from azide and acetonitrile.

It has been shown that zinc salts are excellent catalysts for the tetrazole reaction.¹³⁵ They work well even in aqueous media, allowing for an environmentally friendly large scale protocol to 1H-tetrazoles.¹³⁶ In these cases, simple heating (80-170 C) of an aqueous reaction mixture of nitrile, sodium azide, and catalytic zinc salt provides the 1H-tetrazole in good yield following acidic work up.

DFT calculations were performed to uncover the catalytic role of the zinc ion.¹³⁷ Several possibilities were considered. First, the azide anion (N₃^-) may be bound to Zn, and the acetonitrile performs the cycloaddition without coordinating to Zn (optimized transition state structure in Figure 13A). This barrier turned out to be rather high, 35.6 kcal/mol, providing no reduction as compared to the uncatalyzed situation (31.6 kcal/mol for cycloaddition of CH₃CN and CH₃N₃ and 33.8 kcal/mol for cycloaddition of CH₃CN and N₃^-).

Figure 13 Optimized transition state structures for different possible scenarios of Zn-catalyzed tetrazole formation. (A) Azide coordinated to Zn, nitrile attacks from outside. (B) Nitrile coordinated to Zn, azide attacks from outside. (C) Both nitrile and azide coordinated to Zn. The calculated barriers are 35.6, 28.7, and 27.3 kcal/mol for A, B, and C, respectively.

Second, the acetonitrile may be bound to the zinc ion, and the azide comes in from outside (transition state in Figure 13B). The calculated barrier for this reaction is 28.7 kcal/mol, representing a lower barrier as compared to the reactions without the Zn. Third, both the azide and the nitrile molecules may be bound to the zinc ion (Figure 13C). The barrier was found to be 27.3 kcal/mol, quite similar to the second case. These results suggest that the critical element of the catalysis is the activation of the nitrile by coordination to the metal ion. Hence, the Zn ion plays a role similar to that of the proton in the acid-catalyzed reaction discussed above.

To test this hypothesis and to further isolate the effect of zinc, we considered the intramolecular cycloaddition of the azidocyanamide shown in Scheme 7. This system was chosen to eliminate the effects of the charged inorganic azide species, as well as covalent coordination of the zinc and azide.

Scheme 7. Intramolecular Tetrazole Formation Studied To Isolate the Zn Effect

In the calculations, the Zn ion reduced the barrier by 5.3 kcal/mol, giving a barrier of 26.3 as compared to 31.6 kcal/mol (Table 1). From kinetic experiments, a reduction of 5.9 kcal/mol was measured, which is in excellent agreement with the theoretical value. This supports the hypothesis that the coordination of the nitrile to the Zn is the major factor responsible for catalysis.

To examine whether the Lewis acidity of the Zn is the factor responsible for the catalysis, a more potent Lewis acid was considered, AlCl₃. It proved to be a more powerful catalyst, with a calculated reduction of the barrier by as much as 12.5 kcal/mol (Table 1). Experimentally, this reduction was measured as 10.2 kcal/mol, again in excellent agreement with the theoretical prediction. Optimized transition state structures for the uncatalyzed, Zn-, and AlCl₃-catalyzed reactions are shown in Figure 14.

Figure 14 Optimized transition state structures for [2+3] intramolecular cycloaddition reactions. (A) uncatalyzed, (B) Zn-catalyzed, and (C) AlCl3-catalyzed. Calculated barriers are 31.6, 26.4, and 19.1 kcal/mol, respectively.

Superoxide dismutase (SOD) catalyzes the dismutation of superoxide anion radicals to hydrogen peroxide and molecular oxygen, and thus protects living cells from toxic oxygen metabolites.¹³⁸ SODs form an important part of biological defenses against toxic oxygen intermediates and radical damage. Although their biological roles are not completely understood, SODs have been shown to prevent inflammation and oxidative damage and are involved in anticancer and antiaging mechanisms.¹³⁹ Four types of SODs are known according to the redox-active metal involved: copper-zinc-, iron-, manganese-, and nickel-containing forms.¹⁴⁰ Mutational defects in CuZnSODs have been associated with familial amyotropic lateral sclerosis (FALS), and defects in MnSOD have been associated with neurodegenerative diseases and cancers.^141,142 CuZnSODs are generally found in the cytosol and peroxisomes of eukaryotic cells (and also in the periplasmic space in some bacteria), FeSODs are found in prokaryotes and plants, and MnSODs are found in prokaryotes and in the mitochondria of higher organisms. Because mitochondria use over 90% of the cell's oxygen, the mitochondrial electron transport chain produces a large quantity of oxygen radicals, and MnSOD is a primary biological defense against radical damage.^139,142,143

All of these enzymes function by dismuting the superoxide anion through the net reaction:

We have focused on the comparative redox properties and associated catalytic cycles of Mn, Fe, and CuZn SODs. First, we used small models containing mainly the first-shell side-chain structures for MnSOD from human and bacterial T. thermophilus enzymes^41,144 and CuZnSOD from bovine erythrocyte.¹⁴⁵ Models have been enlarged to include the second-shell H-bonding partners for FeSOD from E. coli,⁴² MnSODs from T. thermophilus, E. coli, and human mitochondria,⁴² and the MnSOD mutant Q143N from human mitochondria.⁴² We have established a combined density functional (DF)/electrostatics method for the coupled redox potentials (that is, a one-electron redox event coupled with a single protonation or proton transfer) of these systems. Recently, we obtained the coupled redox potentials for Mn- and FeSOD based on analysis of the measured kinetic rate constants of Mn- and FeSODs. This scheme is a valuable and practical alternative to direct electrochemical measurements of redox potentials in SODs (see ref 42 for details). We have examined some aspects of the catalytic reaction pathways for Mn, Fe, and CuZn SODs involving electron and proton transfer based on both experimental data and DFT electrostatic calculations.

MnSODs and FeSODs are normally grouped into one class because of their identical set of coordinating ligands, similarities in sequence and protein structure, particularly for the first- and second-shell residues around the active site, and similar catalytic pathways for superoxide dismutation.

Letting "M" represent Mn or Fe, a general reaction scheme outlining the catalytic dismutation of superoxide ion (O₂^-) via alternating reduction of the M³⁺ and oxidation of the M²⁺ SOD enzyme is given below:^40,42,146

-

-

Scheme 8

The complete Mn- and FeSOD protein is a tetramer, comprising two subunits related by a dyad axis. Each subunit consists of two chains, and each chain contains one metal active center. The resting forms of Mn- and FeSODs have active sites (see Figure 15) organized as approximate trigonal bipyramids with three histidines, one aspartate, and one hydroxyl (or water) ligand forming a five-coordinate metal site with one empty site. Our DF/electrostatics calculations show that the oxidized and reduced Mn- and FeSODs are in the [Mn³⁺- and Fe³⁺(OH^-)] and [Mn²⁺- and Fe²⁺(H₂O)] forms, respectively. In the 5-6-5 coordination reaction scheme proposed by Lah et al.⁴⁰ for FeSOD (and applied also to MnSOD), the reaction in the first and second step proceeds by addition of O₂^- to the empty (sixth) site, with inner sphere electron transfer to or from the metal ion.

Figure 15 The active site of Mn- and FeSODs. Labels are for wild-type human MnSOD.

In the first half cycle of Scheme 8 (eq 16), a superoxide anion O₂^- first binds at the sixth ligand site and shares an electron with the M³⁺ center of M³⁺(OH^-)SOD. It then donates this electron to the metal center and leaves as an O₂ molecule. The rate constants describing the combination and the decomposition of M³⁺(OH^-)-O₂^- are k₁ and k₂, respectively. Bull et al.^147a obtained (see Table 2) k₁ = 1.5 × 10⁹ M^-1 s^-1 and k₂ = 2.5 × 10⁴ s^-1 for the catalysis by T. thermophilus MnSOD, and Hsu et al.¹⁴⁸ reported k₁ = 2 × 10⁹ M^-1 s^-1 and k₂ = 8 × 10⁴ s^-1 for the catalysis by wild-type human MnSOD (hMnSOD). The large values of k₁ are near the diffusion-controlled limit. All kinetic studies for MnSODs show that k₂ is relatively small, with the back reaction (k_-1) competing with product dissociation (k₂). We propose that the superoxide anion O₂^- cannot donate the electron to the M³⁺ center and leave as an O₂ molecule immediately upon binding, because the presence of the negatively charged -OH^- ligand inhibits this electron transfer. After binding with O₂^-, the metal center becomes less positive. Consquently, the M-OH^- distance becomes longer, and the net negative charge on group -OH^- pulls the metal and its ligands toward the side chain of Gln (see Figure 15, Gln143 in hMnSOD). Our DFT geometry optimizations show that the average H-bonding distances (between heavy atoms) of HO^-···N-Gln are 3.08 and 3.44 Å, respectively, in the oxidized resting form Mn³⁺(OH^-) and Fe³⁺(OH^-) SODs.¹⁶⁴ This H-bond is much shorter on average by 0.44 and 0.48 Å, respectively, in the Mn²⁺(OH^-) and Fe²⁺(OH^-) SODs.⁴² Finally, we proposed that the reaction involves successive proton transfer along a H-bonding chain. The proton H₂₂ on the side chain of Gln (143 in hMnSOD) will transfer to the -OH^- ligand, the H on the hydroxyl group of the Tyr34 side chain will go to N₂ of Gln143, and the O(Tyr34) site will extract a proton from its surrounding water molecules. With a water molecule now bound as the fifth ligand, the -O₂^- group will easily transfer an electron to the metal center and leave as an O₂ molecule, giving the M²⁺(H₂O) form of MSOD depicted in the first half of Scheme 8.

The assumption of proton transfer when O₂^- binds to the M³⁺(OH^-) center is also consistent with spectroscopic data on the azide adducts of Mn³⁺SOD.^149,150 The Mn³⁺SOD-azide complex exhibits temperature-dependent absorption and circular dichroism (CD) spectra (thermochromism) changes from a lower temperature six-coordinate form to a higher temperature five-coordinate form. It was proposed that azide binding causes a proton to be transferred either to the hydroxyl group, forming a water molecule, or to be further transferred to reside on the aspartic acid.¹⁵⁰ The five-coordinate form may therefore result either from desolvation of the water molecule or from dissociation of the protonated carboxylate group.

In the second half of Scheme 8 (reaction 17), a superoxide anion O₂^- first binds to the M²⁺ ion of M²⁺(H₂O)SOD. The kinetic rate constant for this step is the same as k₁ (k₃ = k₁) (see Table 2), within the limits of the kinetic simulations. The superoxide then either leaves as O₂^- (rate constant k_-3) or abstracts two protons (probably one from the ligand H₂O and one from the second-shell residue Tyr) and dissociates as a H₂O₂ (rate constant k₄). For both T. thermophilus and the wild-type hMnSODs, k₄ = k₂ from kinetics analysis, and both are relatively small. As in the first half of the cycle, the transfer of protons also controls the rate for generating the product H₂O₂.

An inactive form (X) of MnSOD has been observed in the second half of the catalytic cycle, as depicted in Scheme 8.¹⁴⁷ During our DFT geometry optimizations on the larger Mn³⁺(H₂O)SOD active site clusters,⁴² we found that one of the protons of the water ligand, which originally H-bonded to the atom O₁ of the Asp ligand (see Figure 15), gradually moved closer to O₁ and finally transferred to O_1. After transfer, this proton remained H-bonded to the oxygen atom of the OH^- ligand. The Mn³⁺-O₂(Asp) distance then becomes the longest of the metal-ligand distances in the optimized Mn³⁺(H₂O) cluster.⁴² No such proton transfer was observed during the geometry optimizations for the Fe³⁺(H₂O)SOD active site cluster.⁴² We therefore propose that the inactive Mn²⁺-O₂^- complex is formed after such a proton transfer from the H₂O ligand to the oxygen atom of the Asp ligand. When a superoxide anion binds in the end-on form with the Mn²⁺(H₂O), the metal center will possess partial Mn³⁺ character. Further, if an electron transfers from Mn²⁺ to the -O₂^- group to form O₂^2-, the metal changes to Mn³⁺. The bound solvent then has lower energy as a hydroxyl group than as a water molecule. If one of the protons of the H₂O ligand transfers to the H-bonded oxygen of the Asp ligand (rather than to the O₂^- group), the Asp will probably dissociate, and the O₂^2- group can change to the side-on form.^147a The dead-end complex X of Mn³⁺-O₂^2- is thus formed. A feasible alternative mechanism is H₂O dissociation from the Mn²⁺ center.

The Trp (Trp161 in hMnSOD) indole ring which lies over the active site probably has an important, although indirect, effect on H-bonding in the active site. For both Mn- and FeSOD, the atom O₁ of the Asp ligand is H-bonded with the backbone NH group of its neighboring Trp residue, as well as the water ligand. For instance, in the crystal structure of hMnSOD (PDB code: 1N0J), the distances of O₁(Asp159)···O(water ligand) and O₁(Asp159)···N(Trp161) are 2.901 and 3.085 Å in chain A, and 2.985 and 3.163 Å in chain B, respectively. Here, we see the importance of Trp161 in balancing the interaction between the water ligand and the Asp159 side chain. It has been found in Silverman's group that mutating the Trp161 to other residues,^147b,c especially to Phe (W161F), will increase the product inhibition. The replacement with Phe has only a minor effect on the first half cycle (eq 16) involving reduction of the manganese. However, this mutant exhibits strong product inhibition in reaction with the reduced enzyme form.^147b The inhibited state was very similar to that observed for the inhibited wild-type enzyme.^147b We think one possible reason is that the H-bonding balance of Asp159 with the ligand water and residue-161 is altered in the mutants. As compared to Phe161, the Trp161 side chain will have a stronger -cation interaction with the Mn ion and its first-shell ligands, particularly the bound H₂O. This, in turn, will shift the position of the Trp161 main-chain peptide NH. In mutants lacking a strong -cation interaction, the strength of the interaction between the water ligand and Asp159 may increase. Therefore, once the superoxide binds with the Mn²⁺(H₂O) center, the proton transfer from the water molecule to its H-bonding partner atom O₁ of Asp159 will occur more easily than the wild-type form. It was reported from the X-ray crystal structure that the H-bonds between some of the residues, including the Mn-bound water molecule, Gln143, Tyr34, a water molecule, and His30, are strengthened in W161F, and the distance between the bound water molecule and the Mn center is increased in the mutant.^147b Cabelli et al.^147b have proposed that the positioning of Tyr34 may be critical, because this can serve as a proton donor to the metal-bound solvent (OH^-).

To investigate the function of the residues Tyr34 and Gln143 (see Figure 15), which are in the H-bonding chain, the mutants Y34F hMnSOD (Tyr34 Phe34)¹⁵¹ and Q143N hMnSOD (Gln143 Asn143)¹⁵² were characterized kinetically and spectroscopically. The H-bonding chain is broken in the mutant Y34F. In the active site of Q143N hMnSOD, a new water molecule was found to fill the cavity created by the Gln143 Asn mutation.¹⁵² In chain A of Q143N, the side chain -OH of Tyr34 is pushed away 0.9 Å by this new water molecule and hence no longer forms a direct H-bond with the side chain of Asn143. Both Tyr34 and Asn143 are now H-bonded to the new water molecule, which also H-bonds to the original H₂O or OH^- ligand. The kinetic results suggest that the replacement of Tyr34 Phe does not affect the diffusion-controlled steady-state constant k_cat/K_m, which has a value near 10⁹ M^-1 s^-1 for both wild-type and Y34F hMnSOD.¹⁵¹ However, the Tyr34 Phe replacement does affect the rate of maximal catalysis k_cat, reducing by about 10-fold the steps that determine k_cat.¹⁵¹ Comparing the kinetic rate constants in the forward direction (k₁, k₂, k₃, and k₄) for the wild-type (wt) and the mutants Y34F and Q143N hMnSODs (see Scheme 8 and Table 2), we see that k₁(Y34F) = k₁(Q143N) = k₁(wt) = k₃(Y34F) = k₃(Q143N) = k₃(wt), k₂(Y34F) < k₂(wt), k₄(Y34F) < k₄(wt), and k₂(Q143N) = k₄(Q143N) < k₂(wt) = k₄(wt). These results show that the three enzymes have the same rate for the combination of superoxide anion with the Mn²⁺ or Mn³⁺ center. The main differences of the catalysis by these enzymes are within steps 2 and 4, where the catalysis by Y34F and Q143N is much slower than that by wild-type hMnSOD, which means that the alteration of the proton-transfer pathways reduced the rates for producing the products O₂ and H₂O₂. This further supports our proposal that gated proton transfer will happen in both steps 2 and 4 and the proton transfer will occur prior to or concerted with the electron transfer from the superoxide -O₂^- group to the M³⁺SOD metal center in the first half of the reaction cycle.

Edwards et al.¹⁵³ recently reported that mutation of either His30 or Tyr174B in E. coli MnSOD reduces the activity to 30-40% of that of the wild-type enzyme. These two residues are highly conserved and correspond to His30 and Tyr166B in hMnSOD, His32 and Tyr173B in T. thermophilus MnSOD, and His30 and Tyr163B in E. coli FeSOD. Silverman and co-workers^{151,152,154,155} found that mutation of any residue in the H-bonding network of Gln143···Tyr34···H₂O···His30···Tyr166B in hMnSOD will decrease the catalytic activity. This extended H-bonding network may therefore be the complete proton-transfer pathway.

Although FeSOD and MnSOD proteins are highly homologous, they exert very different redox tuning on the active site metal ion.¹⁵⁶ The reduction potential, E_m, of FeSOD is higher by several hundred millivolts than that in Fe-substituted MnSOD (Fe_sub(Mn)SOD), and much higher than that in MnSOD.¹⁵⁶ Because metal ion reduction is intimately coupled to proton uptake in SOD,^42,157 the degree to which protonation of coordinated solvent is favored in each oxidation state contributes significantly to the observed E_m. Thus, stronger hydrogen bond donation to the metal-ion-bound solvent has been proposed to strongly depress the E_m value of the bound metal ion in the MnSOD protein.^156a,b,158 Very recently, Miller and co-workers examined the active site structures of FeSOD and MnSOD from E. coli and concluded that the E_m's can be tuned by the precise positioning of the conserved second-sphere Gln.^156c Their DFT calculations show that, in the oxidized state, the distance between the Gln amide hydrogen and the oxygen of the Fe³⁺-bound OH^- is decreased from 2.70 Å in FeSOD to 1.94 Å in the Fe_sub(Mn)SOD model, which suggests formation of a stronger hydrogen bond in Fe_sub(Mn)SOD. However, protonation of the solvent ligand upon Fe reduction results in considerable steric interference between the Gln and the H₂O protons in the Fe_sub(Mn)SOD model. This indicates that the closer active site Gln side chain in Fe_sub(Mn)SOD than in FeSOD destabilizes coordinated H₂O versus OH^-, thus strongly favoring the Fe³⁺ state and lowering the E_m value in Fe_sub(Mn)SOD.

CuZnSOD is a homodimeric protein in which the two active sites act independently.¹⁵⁹ Each subunit contains one Cu²⁺ and one Zn²⁺ ion in the oxidized ("resting") state. During catalysis, copper is the redox partner of the superoxide radical, whereas the oxidation state of Zn²⁺ does not change during the dismutation reaction. The interaction between the two subunits in the dimer is probably not important for the catalytic process because there is a large distance between the two remote copper sites (~34 Å).

The Cu²⁺ and Zn²⁺ in each subunit in the oxidized form of CuZnSOD are connected by a histidine imidazolate bridge (His61 in bovine erythrocyte CuZnSOD) (see Figure 16). The role of this histidine bridge and its protonation is critical for the catalytic cycle, and this will be the main focus of our analysis. The Cu²⁺ ion is coordinated by nitrogen atoms of His46, His61, and His118 and the N of His44 (and more distantly by a water molecule) in a square planar geometry with tetrahedral distortion. The Zn²⁺ site has the geometry of a distorted tetrahedron created by three N-coordinated histidines and an aspartate residue: His61, His69, His78, and Asp81. The bridging ligand His61 lies nearly on a straight line between Cu²⁺ and Zn²⁺, which are separated by 6 Å.

Figure 16 Representation of the active site in the oxidized CuZnSOD. Ligand numberings are from bovine erythrocyte CuZnSOD.

The active site Cu²⁺ ion is at the bottom of a deep and narrow (<4 Å) channel in the protein, while the Zn²⁺ ion is completely buried. There is no evidence that the potential substrates can bind to the zinc center; rather, the Zn²⁺ ion appears to play both a structural role and a role in electronic polarization and electrostatic stabilization based on our DFT electrostatics studies.¹⁴⁵ While the Zn²⁺ ion is not directly involved in the catalytic process, it probably participates in the catalysis indirectly by interacting with the Cu center and the His61 bridge. The positively charged copper and zinc atoms in the active site, together with other positively charged residues, form a region of strong positive potential which attracts the superoxide to the copper. The binding of the superoxide to the Cu²⁺ center may be also assisted by a positive arginine residue (Arg141) which is located about 5 Å from the Cu²⁺ ion.¹⁶⁰

Mechanistic studies using pulse radiolysis established that, similar to the mechanism of Mn(Fe)SOD, the copper center is reduced and oxidized by superoxide in two steps.¹⁶¹ Superoxide first reduces the Cu²⁺ ion, yielding dioxygen, and then another molecule of superoxide oxidizes the Cu⁺ ion to produce hydrogen peroxide. The protons required to form the product are taken from the aqueous medium.¹⁶² The catalytic reactions are nearly diffusion limited, occurring at a specific rate of 2 × 10⁹ M^-1 s^-1 at the copper site.¹⁶² The overall mechanism can be summarized in Scheme 9, where E symbolizes the rest of the protein.

Scheme 9

In the first step, the protonation of the imidazolate bridge upon reduction of Cu²⁺ to Cu⁺ is accompanied by the dissociation of the His61-N-Cu bond. This leads to an approximately trigonal planar Cu coordinated to three terminal histidines, with the bridging His61-N-Cu bond being broken. In the reoxidation step, the His61 donates the proton to the substrate upon oxidation of Cu⁺ to Cu²⁺ and rebinds to Cu²⁺.

In our DFT calculations,¹⁴⁵ the His61-N-Cu bond lengthens from