J. Phys. Chem. 1993,97, 13083-13088
13083
Protein Electron Transport: Single versus Multiple Pathways J. J. Regan,+ S. M. Risser,* D. N. Beratan,’?*and J. N. Onuchict Department of Physics. University of California-Sun Diego, La Jolla. California 92093, and Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania I5260 Received: July 26, 1993; In Final Form: October 15, 1993’
Pathway analysis provides a tool for the design of tailored electron-transfer proteins and is a useful starting point from which to build up multipathway views of electron mediation that include the influence of interference and all of the chemical tunability that is accessible. We present a new theoretical strategy to determine when a single-pathway picture is sufficient or when one must consider multiple paths. A definition of a single pathway in the context of a Green’s function formalism is presented. To illustrate these effects, examples are given for cytochrome c. We also show that full protein Green’s function calculations can be carried out at the tightbinding (extended-Hiickel) level on cytochrome c including all valence orbitals of the protein. Although many questions remain about appropriate parameterization, the simple Pathway prediction that proteins are not structureless isotropic electron-transfer barriers holds as multiple pathways are included in the calculations and backscattering of electron amplitude within and between pathways is added.
I. Introduction The electron-transfer (ET) problem is central in much of chemistry, biology, and physics. Considerable recent progress has been made in understanding how the structure of a macromolecule controls the nuclear and electronic aspects of the transfer process.1-4 Generally speaking, in long-range macromolecule electron transfer, the rate is proportional to a donoracceptor ( E A ) electronic coupling strength ( TDA)and a nuclear factor associated with the motion along the reaction coordinate. In the last five years, we have developed methods that dissect the interactions between the thousands of protein orbitals into pathways, smaller subsets of these orbitals or effective orbitals. This strategy has generated a satisfactory understanding of a wide range of protein electron-transfer rates and, more importantly, provides a design strategyfor new electron-transfer systems. Pathway analysis511 involves performing a search for the dominant set of bonded and nonbonded interactions between two sites in a protein. These methods generate estimates of relative electronic coupling strengths by applying a simple set of renormalized decay parameters to the best path found between the two points. The reason that Pathway methods work so well12-” is that there is a fundamental differencein the nature of throughbond (TB) and through space (TS) electronic coupling. And the Pathway method was the first to successfully partition protein interactions into these two .classes for analyzing long-range electronic coupling. Pathways, therefore, provides a framework for asking questions about the mechanism of protein electronic coupling mediation. Such questions include the following: What roles can structural motif, quantum interference, and protein dynamics play in electron-transfer processes? However, the question of when the single-pathway picture is valid remainsopen. A new theoretical strategy to address this point is the focus of this manuscript. This paper has two goals. The first is to show how Pathway calculations can help guide the execution and interpretation of exact electronic coupling calculations within an effective twostate model of electron transfer. An extended discussion comparing single and multiplepathways in the context of cytochrome cis presented. Thesecond goal is todemonstrate that full-protein Green’s function calculations can be carried out readily at the t University of t University
California.
of Pittsburgh.
*Abstract published in Advance ACS Abstracts, November 15, 1993.
tight-binding (extended-Huckel)level on cytochromec including all valence orbitals of the protein. All of the methods of calculating TDAthat are described here have a common and critical element. All predict that the protein is an inhomogeneous barrier to electron tunneling. As such, TDA is found to be sensitive not only to the donor-acceptor distance but also to the structure of the intervening protein. Our singleand multiple-pathway models all predict considerable scatter of log(le1ectronic couplingl) vs distance, reflecting the fact that the balanceof TB and TS interactions is not homogeneous in proteins. Full proteins can now be analyzed using Green’s function multipathway18techniques on a workstation or larger computer. This battery of new tools for analyzing electronic couplings and dissecting pathway contributions provides a means to design proteins with tailored ET rates (as well as to design proteins that probe weak points of the theory) and to establish a general description of how protein motifs control ET rates. All of the new methods discussed here are implemented with preliminary parameter sets, so the specific predictions are qualitative. The goal is to work on a regime of coupling energies relevant to proteins and to show the importance of various effects on the tunneling matrix elements. Considerable further study of parameters is clearly needed before quantitative estimates can be made.
Why Are Most hoteins Inhomogeneous Barriers to Electron Tunneling? 11.
Much of biological electron transfer falls in the nonadiabatic limit, and upon separating electronic and nuclear motion, a rate proportional to an electronic tunneling matrix element (TDA) and a nuclear Franck-Condon factor (FC) is found (eq 1). In the
last 10 years, significant effort focused on determining the dependence of TDAon the structure of the chemical (see refs 19-27, for example) or biochemicaP.4 bridge between donor and acceptor. While interesting questions remain concerning the separability of the rate equation, particularly for very long range processes, the approximate parameter scale is well defined by the energetics of the bridge and redox states involved.28
0022-3654/93/209713083304.00/0 Q 1993 American Chemical Society
Regan et al.
13084 The Journal of Physical Chemistry, Vol. 97, No. 50, 1993
Electronic coupling through a one-dimensional constant potential barrier (1DSB) drops approximately as
Here EL^ is the energy of the localized state that is decaying through space (TS). For a 10-eV bound state (all distances in
A)
fully assisted in the design of new ET proteins. For a single pathway consisting of a mix of covalent (C), hydrogen bonded (H), and through-space (S)interactions
and the following parameter set is used: eC
T g a exp[-l.7RDA]
(2b) However, when this localized state is coupled through bonds (TB), with an interaction energy y between adjacent bonds,
where Rl is the projection of the bond length along the DA axis and E ~ NisDthe energy of one of the mediating bonds. (For the sake of this discussion, all bonds are taken to have the same energy.) If the energy of the localized (LOC) state is 2 eV from the energy of a chemical bond (and the interaction between neighboring bonds is about 1 eV), the fall of TDAwith distance through the bonded medium is considerably softer:
Ti:
0:
exp[q.7RDA]
(3b)
Although simple, this analysis describes the trade-off between TB and TS mediation mechanisms in much of chemistry.29 The complexity of “real” systems will arise from a range of values for EL^ and orbital symmetry effects on y,but the drastic separation of distance scales for TB and TS interactions pervades the longrange electron-transfer problem. This observation provided the and it will be seen starting point for our Pathway tocarry through to multiple-pathway calculations with both simple and more complicated tight-binding Hamiltonians. It is important to note that the early models for protein electron transfer used 1DSB models with single-exponential decay constants. However, these decays were chosen to model bondmediated coupling by choosing decay constants smaller than 1.7 A-1.30.31 If chemical and biochemicalbridges were simple isotropic lDSB potentials and all tunneling energies were the same, the drop of electron-transfer rates with distance in all systems would be identical. Not only is the average decay expected to vary with structure, but the decay is also predicted to be rather anisotropic.* Indeed these expectations are confirmed experimentallyl2-17 in both chemical and biochemical systems. 111. The Pathway Model: Atom-Centered and Bond-Centered Approach
The Parhway model for calculating TDAin proteins was developed on the basis of earlier studies of electronic coupling in model compounds.23J4 Pathways partitions electronic interactions in the protein into three types: covalent, hydrogen-bonded, and through-space. This division was based on the fact that TB interactions are much longer range than TS interactions.28 As already described, the barrier to tunneling through a bonded medium is considerably lower than tunneling through a vacuum, so the exponential decay of the bond-mediated coupling is slower than through-space. Also, since hydrogen bonds bring lone-pair and bonding orbitals into close proximity, we expect their mediation properties to be substantial. The nature of bonded and nonbonded interactions is obviously a functionof atom type, hybridization,andorientation. However, since the distinction between bonded and nonbonded interactions is so strong, an understanding of coupling pathways arose from noting the necessary balance of these three interaction types along the dominant mediation pathways in proteins. The simple pathway model has survived experimental tests and has success-
= 0.6, fH =: e
exp[-1.7(R - 2Rc)], es = (1/2)tc exp[-1.7(R-Rc)]
(5)
In these expressions, the distances, R, are in angstroms and the decay factors, c, are unitless. The referencecovalent bond distance is chosen as RC = 1.4 The hydrogen-bond decay is treated as two covalent bonds (covalent bond plus lone-pair) from heteroatom to heteroatom, allowing the coupling to be adjusted if the hydrogen-bond length is longer or shorter than the reference length. TS jumps are establishedwhen the jump provides a direct site to site coupling that is greater than some scaling factor times the best covalent pathway coupling already in existence between the two sites. The pathway coupling strength between a given donor and all points in the protein is anisotropic. A plot of this coupling with an exponential fit drawn through it reveals couplings that differ by orders of magnitude at different positions around the donor but at the same distance from it.9 To achieve a more realistic model, Pathways was recently improved to utilize bond-centered states, as opposed to atomcentered states. This permits the inclusion of all bonds to hydrogens, all lone-pair orbitals, and permits the inclusion of an angular dependent orbital coupling (as opposed to a simple distance dependence). Hydrogen bonds in this picture are simple TS jumps between a hydrogen and a nearby “virtual hydrogen” holding a position in a lone-pair orbital. The aforementioned scatter in coupling with distance is not diminished by the inclusion of these hydrogen and lone-pair orbitals (see section V). All of the results described to this point are based on the dominant single (renormalized) Pathway picture.ls This picture will breakdown when multiple pathways become important. When that is the case, a sum of appropriate contributions is needed to estimate TDA. When are these multiple pathways important? What is the definition of a single pathway? To begin to answer these questions, one can expand the concept of a single pathway and introduce the concept of the pathway family in the context of a Green’s function formalism. In some cases, even this renormalized view (based on pathway families) can break down. Coupling between pathway families could become significant and/or multiple families might control TDA.In such cases, no single “route” representation (or family of routes) will successfully describe the coupling. IV. When is a Pathway a Pathway? An appropriate estimate of TDAis made using the Green’s function of the mediating bridge. The Green’s function for a Hamiltonian H is defined as G(E) = ( E - H)-l. One partitions H into three operators H = PA + Hbr + e‘where mAis the Hamiltonian for the space containing only the D and A states, Hbroperates on a space containing only the bridge states, and describes the interaction between these spaces. With this partitioning in mind, an exact expression for GDA(E)can be written as
Protein Electron Transport
The Journal of Physical Chemistry, Vol. 97, No. 50, 1993 13085
where
transition. The terms in the sum are always proportional to tk where k indexes the n states along the path associated with the term and Ek is the coupling between state k and state k + 1 on the path. In this paper, this serves as a definition of “path coupling” (as opposed to full Gbr matrix element coupling). Since all the couplings are less than 1 in absolute value, a path with more steps will generally contribute less to the sum of amplitudes than one with fewer steps-but this is not always the case since the coupling values are not all the same (see below). Thus the path coupling is not strictly proportional to the number of steps. The first term in the expansion corresponds to the route with the largest coupling. Often, this path will also be a path of shortest physical distance through the molecule along a series of atommediated couplings, but not always. Sometimes the coupling between orbitals that do not share an atom (a TS jump) will provide a shortcut that beats a longer TB path. The Pathways model is presumed to work best when only one dominant path, or perhaps a limited family of structurally related paths, is present in the molecule. It is certainly possible that for a given D-A in a protein, a dominant pathway cannot be identified. In such a case there will be many terms (paths) in the GF expansion that have comparable amplitude contributions to the total coupling, and this may indeed lead to a D-A coupling in the protein that simply falls off exponentially with D-A distance. But an interesting situation arises in many proteins when one pathway, or a family of similar pathways (paths that go through the same groups and are of about the same length), makes a dominant contributionto the total amplitudeof the D-A transition. The remaining paths (terms), when included in the sum, only modify the total amplitude by factors that are orders of magnitude less than the total coupling already provided. When a dominant pathway can be identified, the rate can be controlled by changing the protein with amino acid mutations to either improve or degrade the contributions of the pathway. Further, in a protein where no dominant pathway exists, it may be possible to identify a mutation or set of mutations that will introduce a dominant path. The matrix element TDA,being just a single number, provides limited information about the pathways. That is why, despite the availability (within this model) of the exact TDA,a need for calculatingthe first-order terms (paths) in TDAremains important. As described in our earlier work,15J8two different kinds of interferencemay influence TDA:interferencefrom pendant groups off a path or interference from intersecting pathways (loops). The former effect can be eliminated by renormalizing the couplings and site energies of the orbitals in the main path, retaining the physical meaning of the path. Loops, however, do not lend themselves to simple renormalization. When they become important, multiple pathways must be considered.
ll:lf
(m, n = D,A) (7) H,,,is,,the direct D-A matrix element of the two-state subsystem P A , and the sum is over states in the bridge. G b r is the Green’s function for the system consisting only of the bridge states (Le., @ ( E ) = (E - El”‘)-’). Equation 6 is isomorphic to the expression found for a simple two-state system;
Iftheenergies of theD-Astatesaremuch further from theenergies of the bridge states than they are from each other, and if the coupling to the bridge is weak, one can identify a small parameter encapsulating these conditions and expand GDA(E)in terms of it.32.33 If only the zeroth-order term is kept, one finds that GDA(E) in eq 6 is equal to the CDA(E)in eq 8 with y replaced by HddT,(E,,,), where Et,, is the so-called “tunneling energy” of the electron. Et,, is a parameter representing the energy of the donor and acceptor species; it appears below only in the form yc/Etun, since this parameter is central to the Pathways model. At this point the donor-bridgeacceptor system is recognized as an effective two-state system with D-A coupling T D A = HOA(Et,,). This coupling is then used in the rate expression. Gt(Etun)(required in eq 7) is computed by the Greenpath program. It can be computed by inversion of (EtunI-Z P ) . Also, if one has partitioned some orbitals out of the Hamiltonian to examinecertain structural features, the stepwise Green’s function methodla provides a natural way to re-include these states. For example, after examining a best family of pathways, one can use this method to include a layer of surrounding orbitals to see what effect they have, without restarting the calculation with a new Hamiltonian. The matrix elements of Gbr for cytochrome c appearing in this paper should not be confused with TDA. This is an important distinction, since only TDAis meant to be associated with a rate. Letting fl denote coupling from the D-A states to the bridge, and omitting the direct D-A coupling term of eq 7 (it is negligible for long-distance ET in proteins), one has
(9) TDAincorporates the followingparameters: Et,,, the intrabridgestate coupling y (in C) and the D-A to bridge coupling 8. If one imagines a donor coupling by 8 to only one of the bridge states (say state d), and likewise the acceptor (to a state a), then TDA = fl’G:. The matrix element magnitude squared would be, in this case, proportional to a rate. Though this picture is not representative of a true D-A configuration, it is sufficient to consider the coupling provided between two points on the bridge. When iron appears below in a matrix element, it is understood that the iron is nor in the bridge, and that the matrix element is being taken with respect to the bridge states the iron couples to. For rate prediction, one must address the problem of how the donor and acceptor are coupled to the bridge (to which states, and with what coupling). An expansion of Gr(EtUn)in terms of localized bridge states (atomic or bonding orbitals) gives a sequence of terms that can each be associated with a particular pathway for ET. Such pathways can take a more or less direct route between the D and A sites: some may meander through the molecule, and some may backscatter and retrace steps over and over again. The terms in this sum can interfere with one another, constructively or destructively, to change the probability (and thus the rate) of the
V.
Case Study: Cytochrome c
Calculations were performed on the X-ray structure of yeast cytochrome c (Brookhaven PDB entry ~ Y C C This ) . ~ ~structure has native histidines at positions 33 and 39, and four computer mutations were created to obtain structures with histidines at positions 58, 62,66, and 72. These structures reflect six of the cytochrome c ET experiments examined by Gray and coworkers.35J6 Seven internal waters were left in the structure, but the solvent exterior to the molecule was discarded for the calculation. States associated with hydrogen atoms were taken into account, as were states associated with lone-pair orbitals. To provide a connection between the Green’s function calculations and Pathways, the parameter yc/Etu, = -0.48 was chosen so that the decay per C-C bond in a long alkane chain (with only covalent coupling) was equal to cc = -0.6. The coupling for TS
13086 The Journal of Physical Chemistry, Vol. 97,No. 50, 1993
Regan et al.
TABLE I: Pathway Guided Bridge Coupling‘ HisX His33 11.9 A
no. of oaths
no. of plus states nbrs
1 53 25579 191122
15 45 77 113 2144 14 35 59 85 2144 16 37 66 100 2144 17 72 113
His 39 14.6 A
1 24 7787 63436
His 58 15.3 A
1 25 7886 67867
His62 15.9 A
1 174 86333
His66 13.3 A
1 47 18679 142466
His72 10.4 A
1 23 10217 88 134
2144 14 37 43 70 2144 10 29 49 78 2144
38 78 113 155 36 64 94 123 39 67 107 148 41 108 164 36 62 68 112 31 SO 81 112
(i%lcblCa-IV) y/Et,=-0.48 4.56 e-7 1.85 e-7 -1.81 e-6 -1.21 e-6 1.04 e-6 5.03 e-7 4.70 e-7 3.98 e-7 3.45 e-7 2.67 e-7 -3.63 e-9 3.94 e-9 4.73 e-8 5.02 e-8 5.05 e-8 -2.41 e-9 -1.64 e-8 3.95 e-9 9.01 e-9 -3.88 e-8 -2.58 e-7 -1.22 e-7 -1.31 e-7 -1.05 e-7 -1.17 e-5 1.53 e-5 1.78 e-5 1.62 e-5 1.71 e-5
IGb’l IGkd 1.o 0.4 4.0 2.6 2.3 1.o 0.9 0.8 0.7 0.5 1.o 1.1 13.0 13.8 13.9 1.o 6.6 1.6 3.6 1.o 6.6 3.1 3.3 2.1 1.o 1.3 1.5 1.4 1.5
state choice best SO% 10% 5%
all best 50% 10% 5% all best 50% 10% 5% all best 50% 10% 5% all best 50% 10% 5% all best 50% 10% 5% all
0 Gbr matrix elements computed over various subsets of states in the protein for six different acceptors. All the states in the fattened best pathway comprise the first subset (for each acceptor). All the states in the fattened paths with a path coupling of at least 50% of the best-path coupling comprise the second subset, etc. The last line for each acceptor shows the result for including all protein states. Note sign changes in the matrix element. Column six compares all numbers to the best path in each case.
jumps was set to yc exp(-1.7(r - 1,4)), with r being the bond center-to-center distance in angstroms. This choice may overestimateTS coupling because at this point no angular corrections to orbital overlaps have been included. Also, a standard covalent bond reference length (of 1.4 A) may be too long, as is certainly the case with bonds involving hydrogens. However, although further elaboration of these parameters is called for, this set is sufficient to retain the important distinction between TS and TB couplings, as proposed in the pathway model. In Table I, matrix elements of Gbrbetween bonding orbitals associated with select C i s and bridge states directly coupled to the heme iron are presented (this would be appropriate if the donor orbital were considered to be fully localized on the iron). Table1 addressesthequestionposedearlier. Whenisa pathway in the bridge a meaningful concept? How does the contribution of the states in a path family compare to the coupling provided by the entire protein? These numbers are Gbr matrix elements, but over different subsets of the bridge states. The state subsets were chosen by performing path searches. A path search is a depth-first graph search, which yields only “forward” paths, paths that never visit a state more than once (as opposed to paths with backscatter or loops). The “best” path is the path with the highest magnitude path coupling (nrek). The best paths from the iron to the C,-N bonds of the six histidines were found. This path generally consistedof about 14-17 states (except for His72, which has only 10 steps in the path-see “no. of states” column in Table I). To incorporate most of the effects due to pendant groups when computing the Gbr matrix element for a single path, all bonds sharing a common atom with path bonds are included in a Gbr
Figure 1. Backbone ribbon tract of cytochrome c with the heme, His33 and His39, highlighted. The pathways (dark lines) shown provide at least 10% of the best-path coupling from the iron to each His (see Table I). The paths to His39 are relatively confined (except for one weak outlier) compared to the paths meeting the same criterion in the case of His33.
calculation over a pathway. This is called a “fattened” pathway. Each best path was fattened, resulting in sets consisting of 36-41 states (31 for His72-see “plus nbrs” column in Table I). The Gbr matrix element was then calculated over this expanded subset. Successively larger sets of states were selected by computing all paths that provided some percentage of the path coupling of the best path and then fattening these path sets by including any states with direct covalent links (no TS jumps) to path states that were not themselves already included in one of the paths. The percentages used were 50% 10%. and 5%, and these were followed by a calculation which included all orbitals in the protein (see the last line for each HisX). Thisprocessgeneratedlargenumbersofpaths(191 122distinct forward paths have coupling that is greater than or equal to 5% of the best path in the case of His33). Even so, the number of states involved is never higher than 10% of all the states in the protein. The matrix elements for each HisX are sorted by the number of states included in the calculation. One can now ask how these state subsets (reduced proteins) compare to the entire protein. If the best pathway is not dominant, one expects to see both a change in the magnitude of the matrix element as the entire protein is included and sign changes in the matrix elements as interference effects come into play. Although no striking difference was found in the magnitudes shown here, sign oscillations are evident in some cases. An example of this is provided by the two native histidines at 33 and 39, which are shown along with the heme and a backbone trace in Figure 1, The dark lines are iron to C,-N paths from the ”10% of best” level in Table I. The 7787 paths to His39 (meeting the 10% of best criterion) appear in a fairly confined tube (pathway family). No sign oscillation is observed here. On the other hand, the 25 579 paths to His33 (meeting the same criterion) follow a number of distinct routes, providing an opportunity for a greater degree of interference. Figure 2 presents a further investigation of the parameters in the model. A set of plots present bridge coupling from the iron (Le., from states directly coupled to the donor) to all the Cis in the protein. Two values of -yc/Etunare used, and a fitted line is shown for each. Scatter about this line over orders of magnitude at thesamedistances from thedonorisindicativeof theanisotropy in the protein due to its structure. From comparison of these results with the Pathway result, it appears that a -yc/E,, =
Protein Electron Transport -4
xx I
k
The Journal of Physical Chemistry, Vol. 97, No. 50, 1993 13087 100
0 )
#
X
\
*
6
-8
h
*-M
-
.10
.-a a
g
’
-12
2
-14
Distance to alpha-carbon Figure 3. Magnitude of the extended-Hiickel Green’s function matrix elements between the heme Fe (calculated as described in the text) and the sp3C . of each amino acid directed toward the aminoacid’s side chain. The tunneling energy is taken to be 1/4 of the HOMO-LUMO gap energy above the HOMO.
-18
t
.an -” I
8
39
Pathways vs. Green X
0
V
Best path to C, (top line) (FelGICJ, “(/Et,, = -0.556 (FelGICJ, *I/EU, = -0.48 10
12
14
v\ V
V
V
18
18
Distance from 1H:HEM:FE Figwe 2. Magnitude squared of G” matrix elements to C, orbitals, plotted against direct distance from the heme iron. Results for twovalues of y c / E t , are shown, along with a third set of points showing Puthwuys results. In all cases, scatter about the fitted line is indicative of the anisotropy of the protein. Three acceptors addressed in the text (72,33, 39) are highlighted with oversized icons. Note that 72 is below the b t fit line for the pathway case (predicting a weaker than average coupling for the distance it is from the donor) and that this agrees with the yC/Etun = -0.556 matrix element. The parameter y c / E t , is being investigated, and results for the -0.48 case are shown as well. Note that the couplings for these three acceptor sites is similar in the Pufhwuys and y c / E t , = -0.556 case, despite the differences in their distance from the donor. -0.556 is more appropriate than the -0.48 used in the results
above. In the latter case, only the covalent couplings in alkane were used in the calibration. The results for Pathways and Cbrqualitatively agree for the -0.556 case. His72 is below the line as expectedI4 (a major TS coupling dominates the pathway). The values for His72, His33, and His39 are not very different, as was observed experimentally. 14.35.36 Sections IV and V of this paper serve to take the pathway model another step forward. In the model, one can always find a “best pathway”, but such a pathway may or may not be critical to an ET reaction. We have proposed to further define an ET pathway as a small subset of bridge states that provide a bridge coupling (in the context of an effective two-state system) on the same order of magnitude as the entire molecule. This could be called a “reduced protein”, in the sense that the remainder of the molecule does not substantially change bridge coupling when included. Qualitatively,for the idea of a pathway to have meaning, this subset of states should appear somewhat “tubular” in a physical model of the molecule.
VI. Extended-Hlickel Parameterization The results shown so far used a Hamiltonian with very simplified parameters. It is also possible to exactly calculate the Green’s function for a protein using a tight-binding Hamiltonian with more diversity in the coupling parameters. For example, the specific Hamiltonian that we have used for this is the extended Hiickel Hamiltonian, which is a single-particle Hamiltonian that explicitly treats the valence electrons. A minimal basis set of four valence orbitals per heavy atom and one per hydrogen atom
is used. The orbitals are taken to be of Slater form and are nonorthogonal. The Hamiltonian matrix is formed by calculating matrix elements of these basis orbitals. Standard orbital parameters were used for all the elements, with the cutoff for overlap set at 10.0 A. The diagonal matrix elements are the ionizationpotentialsoftheatomicorbitals,andoff-diagonal matrix elements are proportional to the average ionization potentials of the two orbitals and their overlap.3’ Because the basis is nonorthogonal, the Green’s function is not simply the inverse of (EI - H) but The orbital overlap matrix is S. The tunneling energy is determined from the energy of the porphyrin HOMO. For cytochrome c there are approximately 4500 orbitals, which necessitates inverting a 4500 X 4500 matrix. It is feasible to perform such large-scale inversions directly, without eliminating amino acids using a screening meth0d.3~ Similar calculations have been performed with an independent-particlemultiorbital Hamiltonian in a spherical orbital/united atom limit, which decreases the size of the numerical task.39 The Green’s function matrix gives the electronic coupling between any pair of orbitals in the protein, while what is desired is the total coupling between the donor and acceptor sites. To approximate the contribution of the summation in eq 9, we sum the Green’s function matrix elements between the six hybrid orbitals surrounding the Fe and the C,sp3 hybrid orbital pointing in the direction of the first atom in the side chain. Figure 3 is a scatter plot showing predicted matrix elements vs distancein yeast cytochromec using the standard extended-Hiickel parameters and the approximation to the tunneling energy describedabove. As with Pathways, and with theGreen’sfunction calculations in section V, scatter is quite pronounced, showing the importance of the detailed protein environment.
VII. Conclusions The broad distribution of chemical interactions (covalent, hydrogen-bonded, van der Waals, ionic) within proteins makes them rather inhomogeneous barriers to electron tunneling. This is clear in the fact that all reasonable estimates of 7c/Etunand ( 2 m J 3 ~ ~ ) 1 / * lead / h to scatter of the coupling vs distance plots, whether calculated with Green’s function or Pathways methods. An excitingchallengein the field of electron-transfer chemistry it to determine what qualitative differences might arise in the electron-transfer rates as protein secondary structure, tertiary
Regan et al.
13088 The Journal of Physical Chemistry, Vol. 97, No. 50, 1993
structure, D-A redox potentials, and other parameters vary. It will be particularly interesting to see how structural motifs might modify interference effects, dominant pathways, etc. Green's function and Pathway methods provide a means of addressing these issues. Pathways is apparently the first theoretical approach to be used for the successful design of new electron-transfer proteins with prescribed rates.35 The idea of partitioning the decay of the electronic wave function through the protein into two classes (TB and TS) permitted us to simplify substantially the description of the protein environment while capturing the major effects of electroniccontrol. However, severalsimplifyingassumptionswere made in developing Pathways. The most drastic onewas to assume that TDA is controlled by a single pathway or by a few noninterfering paths. Such an assumption appears valid for several systems we have been studying. We can now identify when this is indeed the case and, more importantly, estimate TDA when this assumption fails. In this manuscript we have presented a strategy that can be used to address prior assumptions. Although the results should be viewed as preliminary, the theoretical framework is described and initial comparisons are presented. Using a Green's function technique, we have been able to probe the validity of the pathway concept. Sections IV and V of the manuscript were dedicated to this question. The broadened definition of a pathway as a pathway family shows why our early Pathways approach works at all. Complicating effects arising from backscattering and pendant groups can be included while keeping the pathway picture intact, i.e., TDAcan be computed from a product of decay factors (see eq 4). Again, the parameters used in Pathways should be thought of as renormalized parameters that includesuch effects. Figure 2 shows good correlation between the Green's function results and Pathways results for cytochrome C.
The main goal of this manuscript is to describe the new set of theoretical tools that we have developed to attack the problem of electron-transferin proteins. The parameters used here, orbital couplings and Et"",are preliminary-chosen to give the average decay of coupling with direct distance in cytochrome c, consistent with the Pathways results. However, since TB and TS coupling interactions are so qualitatively different, we continue to believe that a simple set of parameters, that does not differentiate among covalent orbitals, will be sufficient in many cases. Input from quantum chemistry calculations on smaller systems will permit tests of this hypothesis and will clearly play a major role in the future development of parameters. All of our results, whether based on single-, restricted-, or multiple-pathway methods, reflect the overwhelming anisotropy of the protein medium. This anisotropy cannot reasonably be ignored in making estimates of TDA, despite the appealing simplicity of isotropic barrier models.40
Acknowledgment. The authors thank the National Science Foundation (Grant No. MCB-9018768 and an NYI award), the Department of Energy's Advanced Industrial Concepts Division, and National Institutes of Health (PHS Training Grant
GM08326-03), and The Arnold and Mable Beckman Foundation for generous support of this research. Computational resources provided by Cray Research Inc. and Autodesk Inc. (HyperChem) were invaluable.
References and Notes (1) Newton, M. D. Chem. Reu. 1991, 91, 767. (2) Mikkelson, K. V.; Ratner, M. A. Chem. Reu. 1987,87, 113. (3) Structure and Bonding: Long Range Electron Transfer in Biology;
Springer-Verlag: New York, 1991. (4) Metal Ions in Biological Systems; Sigel, H., Sigel, A,, Eds.; Marcel Dekker Press: New York, 1991; Vol. 27. (5) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488. (6) Beratan, D. N.; Onuchic, J. N. Photosynth. Res. 1989, 22, 173. (7) Onuchic, J. N.; Beratan, D. N. J . Chem. Phys. 1990, 92, 722. (8) Beratan, D. N.; Betts, J. N.; Onuchic, J. N. Science 1991,252,1285. (9) Beratan, D. N.; Betts, J. N.; Onuchic, J. N. J . Phys. Chem. 1992, 96, 2852. (10) Betts, J. N. Beratan, D. N.; Onuchic, J. N. J . Am. Chem. Soc. 1992, 114,4043. ( 1 1) Regan, J. J.; Betts, J. N.; Beratan, D. N.;Onuchic. J. N. InPrinceton
Lectures in Biophysics; Bialek, W., Ed.; World Scientific Publishing Co.: Singapore, 1992; pp 175-195. (12) Beratan, D. N.; Onuchic, J. N.; Betts, J. N.; Bowler, B.; Gray, H. B. J. Am. Chem. SOC.1990,112, 7915. (13) Beratan, D. N.; Onuchic,J. N.;Gray, H. B. In Metallom in, Bio/ogical Systems; Sigel, H., Sigel, A., Eds.; Marcel Dekker Press: New York, 1991; pp 97-127. (14) Onuchic, J. N.; Beratan, D. N.; Winkler, J. R.; Gray, H. B. Science 1992, 258, 1740. (15) Onuchic, J. N.; Beratan, D. N.; Winkler, J. R.; Gray, H. B. Annu. Rev. Biophys. Biomol. Strucr. 1992, 21, 349-377. (16) Casimiro, D. R.; Beratan, D. N.; Onuchic, J. N.; Winkler, J. R.;
Gray, H. B. In ACS Books: Mechanistic Bioinorganic Chemistry; Thorp, H., Pecoraro, V., Eds.; ACS Press: Washington, DC, 1993; in press. (17) Beratan, D. N.; Onuchic, J. N. Aduances in Chemistry Series 228; Bolton, J. R., Mataga, N., McLendon, G., Eds.; ACS Press: Washington, DC, 1991; p 72. (18) Onuchic, J. N.; Andrade, P. C. P.; Beratan, D. N. J. Chem. Phys. 1991, 95, 1131. (19) Closs, G. L.; Miller, J. R. Science 1988, 240, 440. (20) Gust, D.; Moore, T. A. Acc. Chem. Res. 1993, 26, 198. (21) Wasielewski, M. R. Chem. Rev. 1992, 92,435. (22) Larsson, S.J . Am. Chem. Soc. 1981, 103, 4034. (23) Beratan, D. N.; Hopfield, J. J. J . Am. Chem. Soc. 1984,106, 1584. (24) Onuchic, J. N.; Beratan, D. N. J . Am. Chem. SOC.1987,109,6771. (25) Jordan, K. D.; Paddon-Row, M. N. Chem. Rev. 1992, 92, 395. (26) Liang, C. X.;Newton, M. D. J . Phys. Chem. 1992, 96, 2855. (27) Naleway, C. A.; Curtiss, L. A.; Miller, J. R. J. Phys. Chem. 1991, 95, 8434. (28) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1985, 83, 5325. (29) Hoffmann, R.Acc. Chem. Res. 171, 4, 1. (30) Hopfield, J. J. Proc. Natl. Acad. Sci. U.S.A. 1974, 71, 3640. (31) Jortner, J. J . Chem. Phys. 1976, 64, 4860. (32) Skourtis, S. S.;Onuchic, J. N. Chem. Phys. Lett. 1993, 209, 171. (33) Skourtis, S.S.;Beratan, D. N.; Onuchic, J. N. Chem. Phys. 1993, 176, 501. (34) Louie, G. V.; Brayer, G. D. J . Mol. Biol. 1990, 214, 527. (35) Wuttke, D. S.; Bjerrum, M. J.; Winkler, J. R.; Gray, H. B. Science 1992, 256, 1007. (36) Casimiro, D. R.; Wong, L. L.; Colon, J. L.; Zewert. T. E.;Richards, J. H.; Chang, I. J.; Winkler, J. R.; Gray, H. B.J. Am. Chem.Soc. 1993,/15, 1485. (37) Yatw, K. Huckel Molecular Orbital Theory;Academic Press: New York, 1978. (38) Siddarth, P.; Marcus, R. A. J . Phys. Chem. 1993, 97,2400. (39) Gruschus, J. M.; Kuki, A. J . Phys. Chem. 1993, 97, 5581. (40) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.;Dutton, P. L. Nature 1992, 355, 6363.
