Electron Transfer Contact Maps - American Chemical Society

J. Phys. Chem. B 1997, 101, 1215-1234

1215

Electron Transfer Contact Maps Spiros S. Skourtis*,† and David N. Beratan* Department of Chemistry, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15260 ReceiVed: June 27, 1996; In Final Form: October 29, 1996X

An exact approach is introduced to establish and compare quantitative structure-function relationships in electron transfer proteins. The strategy is based on the idea of an electron transfer contact map and is motivated by reduced descriptions of 3D protein structure that emphasize the folded conformation (ribbon drawings, distance contact maps, and Ramachandran plots). Electron transfer contact maps render the protein-mediated electronic coupling information in a compact and transferable form. The method includes all contributions of multiple interfering tunneling pathways to the tunneling propagation and can be implemented at many levels of electronic structure theory. Contact maps are used to interpret the influence of mutations on electronic propagation in model peptides.

G(Et) ≡ (EtI - H)-1

1. Introduction Protein-mediated electronic interactions facilitate electron transfer (ET) reactions.1 Pathway analysis2-4 and other methods5-12 reveal numerous effects of structure on ET reactions in proteins,13-17 yet two considerable challenges remain. One is to establish reliable valence-orbital effective Hamiltonians that accurately describe electron and hole mediation. Recent progress in this effort appears promising.18 The other is to extractsfrom the millions of elements in a protein’s Hamiltonian matrixsspecific rules that describe how primary sequence composition and folding-induced interactions influence the ET rate. For this latter task it is clear that a reduction of the information related to the intraprotein interactions is essential. The goal of this paper is to develop an exact method for the comparison of electron tunneling propagation in peptides and proteins that exploits information reduction. In protein ET the rate is nonadiabatic, given by19-22

kDfA )

2π |T |2(FC) p DA

(1)

FC is the thermally weighted Franck-Condon factor between donor (D) and acceptor (A) vibronic manifolds. TD A is the effective electronic matrix element between D and A electronic wave functions. If an orthogonal basis set of local orbitals is used in a calculation, it is given by23

TD A )

vA aGa d(Et)vd D ∑ ad

(2)

The general case of a nonorthogonal basis set is discussed in the Appendix. TD A describes how the intervening medium between D and A facilitates tunneling from D to A by providing virtual orbitals for the transferring electron (at energy Et).23-27 These virtual orbitals are oxidized or reduced bridge orbitals that are not thermally accessible to the electron. In eq 2, d(a) are protein-localized orbitals that couple directly to the D(A) wave functions (through matrix elements vd D(vA a)). The D and A wave functions can be chosen in a number of ways (pure ET cofactor centered, mixed cofactor/protein, etc.). The protein’s influence on the tunneling propagation for a specific |D〉, |A〉 pair is reflected in the Ga d(Et) elements, where † Current address: Department of Natural Sciences, University of Cyprus, PO Box 537,1678, Nicosia, Cyprus. X Abstract published in AdVance ACS Abstracts, January 15, 1997.

S1089-5647(96)01924-4 CCC: $14.00

(3)

is the Green function matrix for the isolated (uncoupled from D and A) protein. H is the Hamiltonian of the isolated protein. In general, any matrix element Gi j(Et)28-32 reflects how electron amplitude that enters the protein at orbital j propagates through the protein to orbital i. As such, all relationships between primary sequence, folded conformation, and tunneling mediation are contained in the Green function matrix. We propose a set of new methods to analyze such relationships. Rather than concentrating on a specific DA pair, we emphasize global features of the G(Et) matrix. The methods developed here may be used at any level of quantum theory that is computationally tractable. That is, H could be a simple extended-Hu¨ckel Hamiltonian, a Fock matrix of a more advanced semiempirical technique (e.g., CNDO), or the Fock matrix of an ab initio calculation. The main body of the paper uses the CNDO method to compute H, which serves as a starting point in our analysis. The CNDO approach sets the overlap matrix to unity in the secular equation. In the Appendix, we rederive all formulas in our analysis so that it can be used with electronic structure methods that do not set the overlap matrix equal to the unit matrix in the secular equation. Our goal is to devise tools to address typical questions of structure-function analysis. Are there particular folded protein motifs that propagate tunneling electrons more readily than others? How important is primary structure in an ET protein compared to its secondary and tertiary structure? For example, are proteins with different amino acid sequences but similar folded conformations comparable ET mediators? Do identical secondary structures embedded in different tertiary motifs have similar ET characteristics? Structure-function analysis of any kind interprets the function of a protein in terms of its primary sequence and/or folded conformation. However, the 3D atomic structure of any protein is complex, and simplified representations are necessary for structural interpretation. Examples of reduced representations of protein structure include ribbon drawings, distance contact maps, and Ramachandran plots.33 A common feature of these representations (Figure 1) is that they (partially) suppress structural complexity while revealing certain essential characteristics. Distance contact maps are often used to study secondary and tertiary structural homologies (or differences) between proteins. A distance contact map is a 2D plot of the distances between © 1997 American Chemical Society

1216 J. Phys. Chem. B, Vol. 101, No. 7, 1997

Skourtis and Beratan

Figure 1. Reduced descriptions of protein structure. (a) The left side shows a ribbon drawing of a 20-residue R-helix. The upper triangle in the 2D plot is a contact map for this R-helix. The right side shows a ribbon drawing of a 20-residue antiparallel β-sheet. The lower triangle of the 2D plot is a contact map for the β-sheet. Ribbon drawings highlight the backbone of each peptide, and therefore they emphasize the folded conformation. The contact maps for the R-helix and β-sheet are computed from all backbone atoms of the second to 18th residues. The numbering of atoms i(j) is from the N to C terminus, and adjacent residues are separated by lines. The shade of gray at position i j represents distance between atoms i and j (distance cutoff is 4.5 Å). The topology of the gray dots between nonadjacent residues reflects the folded conformation of each peptide. For example, the gray dots between residues k and k + 4 in the R-helix reflect R-helical hydrogen bond distances. The differences in the hydrogen bond topologies of the R-helix and the β-sheet are clearly seen in the 2D plot. (b) The left side shows the two torsion angles φk, ψk of each residue k that are often used to describe backbone conformation in a Ramachandran plot. The right side is a schematic representation of a Ramachandran plot for (ψ, φ) conformations found in typical secondary-structural motifs (R: right-handed R-helices, β: β-sheets, L: left-handed R-helices33,34).

atoms in a protein (often backbone atoms). If atom i is within a cutoff distance Ri j from atom j, this is represented by a shade of gray or a cross at position i, j of the 2D plot. When two proteins are compared, their contact maps are first calculated from a subset of atoms that is common to both proteins, and then the two contact maps are juxtaposed (Figure 1a). This description of protein structure is reduced because (a) only a few atoms per residue are shown, and (b) a cutoff distance is often chosen as a criterion for proximity between atoms (or residues). The reduction in the number of variables describing the protein and the choice of the unbranched portion of the backbone as a reference make distance contact maps useful for visualizing and comparing global patterns of folding topology. Ramachandran plots35 describe protein backbone torsion angles, usually two torsion angles per residue. These torsion

angles are plotted in 2D scatter plots (Figure 1b). The information reduction in a Ramachandran plot does not in itself make the plot an approximate description of protein conformation. The Ramachandran plot is as accurate and precise as the input torsion angles. 2. A Reduced Description of Protein ET Motivated by the usefulness of reduced representations, we have constructed a backbone-based description of protein electron transfer (Figure 2). This new view allows exact comparisons of electron tunneling propagation in different protein structures, since the backbone chemical constitution is identical in all proteins (barring length). Our starting point is the Hamiltonian H of the protein, represented in a hybrid atomic

Electron Transfer Contact Maps

J. Phys. Chem. B, Vol. 101, No. 7, 1997 1217

Figure 2. (a) Subspace b of hybrid orbitals for residue k. All remaining orbitals comprising the branched portion of the residue are grouped into subspace hr. These involve backbone hanging group (h) and side chain (r) hybrids. The backbone hanging groups that participate in hydrogen bonding are shown with darker lines. These are denoted hy (section 3). (b) Schematic representation of direct Hb b and effective hb b(Et) matrix elements between b orbitals (only a few b orbitals are shown, as thickened dots). The entire unbranched b portion of the backbone is represented by a thin solid line. The effective matrix elements arise from tunneling through the branched (hr) portion of the protein chain (represented by branching, thick solid lines). Both types of matrix elements reflect the folded conformation of the protein. However, only hb b(Et) elements incorporate the effects of side chain chemical structure on the tunneling.

orbital basis set. The basis set is divided into two subspaces (Figure 2). The subspace b contains hybrid orbitals arising from bonds in the unbranched part of the backbone [-N-CRC-]k, and the subspace hr contains all remaining hybrids that constitute the branched portion of the protein chain (Figure 2a). The latter subspace includes hybrids of backbone hanging groups [:N-H, CR-H, Cd:O:]k and side chains [CR-R]k (denoted h and r, respectively). The hr subspace also includes hybrids of disulfide bridges [-CH2-:S:-:S:-CH2-] when they are present. The exact tunneling propagation between hybrid orbitals of b is given by the b b submatrix of the protein’s Green function matrix G(Et) (eq 3). Using the projection method23,36-40 we find

Gb b(Et) )(EtIb b - [Hb b + hb b(Et)])-1

(4)

Hb b is the b b submatrix of the protein’s Hamiltonian H. It contains all direct interactions between b hybrids (Figure 2b).

The diagonal matrix elements of Hb b are the site-energies of these orbitals. The off-diagonal elements of Hb b describe all covalent and through-space interactions between the b hybrids. The other term, (hr) hb b(Et) ) Hb hrGhr hr(Et)Hhr b, where (hr) -1 Ghr (5) hr(Et) ≡ (EtIhr hr - Hhr hr)

represents effective Hamiltonian matrix elements between b hybrids. Hb hr (Hhr b) is the b hr submatrix of H that contains all direct interactions between the b hybrids and all remaining hybrids belonging to the branched portion of the protein (hr). (hr) Ghr hr(Et) is the Green function matrix of the entire branched portion of the protein; it describes the tunneling propagation through the hr hybrid orbitals. Therefore, hb b(Et) represents the interactions between b hybrids that are mediated only by tunneling through backbone hanging groups and/or side chains (with no intermediate propagation via b hybrids). These



Figure 3. (a) Reference R-helical structure (denoted I) of primary sequence GLY8. (b) Final R-helical structure (denoted II) created from GLY8 by the double mutation GLY1 f GLN1, GLY8 f GLN8. The ends of the GLN side chains are close to each other and to GLY4, GLY5. The positions of the backbone atoms in both structures are identical. The b portion of the backbone (Figure 2a) is represented by a tube.

interactions arise from all possible tunneling pathways through the branched portions of the protein, including hydrogen bond pathways [N-H‚‚‚:O:dC], side chain pathways [CR-R‚‚‚R′CR], and their combinations (Figure 2b). The above description of ET is reduced because it only considers electronic tunneling between certain backbone hybrid orbitals. It is exact because the effect on this tunneling of all remaining orbitals is taken into account fully (via eqs 4, 5). Matrix elements of Hb b and hb b(Et) between orbitals in nonadjacent amino acids (along the protein chain) reflect the folded conformation of the protein. Such matrix elements are denoted long-range, where range means separation along the protein backbone. In addition, hb b(Et) contains the effect of primary structure on tunneling. To interpret this primary structure effect, imagine superimposing the backbones of two proteins that have identical backbone lengths and folded conformations, but different side chains. In this case, differences in the tunneling propagation (Gb b(Et)) among the two proteins arise predominantly from the different electronic structures and branchings of the side chains (not from their relative positions in space). In this situation, the Hb b matrices of the two proteins are nearly identical (differences arise from polarization effects due to the side chains, as shown in section 4). The hb b(Et) matrices are not identical, reflecting the differences in primary structures.

The interpretation of [Hb b + hb b(Et)] is aided by colored 2D plots. The x (y) axis of such a plot represents b hybrid orbital number (measured from the N to the C terminus, with residue identity indicated). The value (or absolute value) of Hi j (or hi j(Et)) is represented by a color at the x ) i, y ) j position on the plot. The 2D representation of these interactions is analogous to the more traditional distance contact map in the sense that it represents “electronic” proximity between different residues rather than “spatial” proximity. This new kind of map is called an electron transfer contact map. The usefulness of this representation is that it shows the global effects of protein conformation and electronic structure on the electron tunneling propagation. In the following sections we give examples of contact map analysis by comparing tunneling propagation for two R-helical peptides (Figure 3). In section 3, we use ET contact maps to interpret typical secondary-structural interactions in an R-helix. Although secondary-structural effects present an important application of the method, we do not give a detailed analysis of these effects in this paper. Rather, we emphasize two other issues (section 4): (i) definition of the perturbation associated with a mutation (in this case the double mutation that relates the two R-helices in Figure 3); (ii) interpretation of the interferences between the new pathways that arise upon mutation and the initial unchanged propagation. We also demonstrate



that drastically different pathway regimes can coexist for a fixed tunneling energy, even in a small peptide. These regimes range from the single to the multiple pathway limit. For all calculations, the peptide Hamiltonian is constructed at the CNDO/S SCF level.41,42 The hybrid orbital basis of each peptide is derived from the converged SCF density matrix, via Weinhold’s natural hybrid method.43 A modification of Weinhold’s BONDO program44 is used that allows computation on larger systems, incorporation of the CNDO/S parametrization, and output of the atomic orbital to hybrid orbital transformation matrix. Figures 4-6 in the following sections are interpolated (smoothed) plots of the tunneling interactions in order to assist visualization. Figures 9-11 are not interpolated. 3. Hydrogen Bond Topology in an r-Helix In this section we discuss how the proposed method can be used to analyze the topologies of electronic interactions in typical secondary-structural motifs. In particular, we concentrate on the model R-helical peptide GLY8 (Figure 3a) and show that for this peptide the effective interactions, hb b(Et), reflect the typical hydrogen bond topology of an R-helix. ET model systems of this general structure are the subject of current interest.45,46 Figure 4a is a 2D plot of the direct interactions Hb b for this peptide (eq 4), where the x (y) axis represents the b orbital number (numbering of b orbitals is from the N to the C terminus, as in Figure 2a). The color of the 2D plot at x ) i, y ) j corresponds to the value of |Hi j| in atomic units (au), as shown in the color legend. Figure 4a shows that the dominant matrix elements of Hb b, shown in green, are short-range (at most, between nearest neighbor amino acids). The strongest elements of Hb b involve covalent interactions between pairs of adjacent orbitals [N-CR-C-]k (magnitude j1.0 au), and the weakest involve interactions between pairs of b orbitals belonging to nearest neighbor amino acids (magnitude j10-2 au). Figure 4b is a similar absolute value plot of the effective interactions hb b(Et) (eqs 4, 5). Et ) -0.2 au (-5.4 eV), a tunneling energy value that is approximately at the middle of the peptide’s HOMO-LUMO gap and is typical of biologically relevant values. The dominant matrix elements in this plot are also short-range (intra-amino acid). Their maximum strengths are an order of magnitude smaller (j10-1 au) than the maximum strengths of the Hb b elements (Figure 4a). These interactions arise from backscattering propagation through the hanging groups and side chains. This backscattering leads to renormalization of the b hybrid orbital couplings and energies. The dominant longer range effective interactions connect amino acids k to k + 2 (magnitude ∼10-3 au), k to k + 3, and k to k + 4 (both magnitudes ∼10-3 au) . The hb b(Et) matrix represents the sum of hanging group [:NH, CR-H, Cd:O:] and side chain [CR-R] mediated interactions between b orbitals. Since hb b(Et) contains the Green function (hr) Ghr hr(Et) for the entire branched portion of the peptide (eq 5), it also captures the tunneling propagation through all typical R-helix hydrogen bonds [N-H‚‚‚:O:dC]. We define the matrix -1 hb(hy) Hhy b b (Et) ≡ Hb hy(EtIhy hy - Hhy hy)

(6)

In eq 6 Hhy hy is the submatrix of the total Hamiltonian H that involves all N-H and Cd:O: hybrid orbitals of the peptide. The N-H and Cd:O: groups of each residue are shown with darker lines in Figure 2a. As defined, hb(hy) b (Et) describes the effective interactions between b hybrids that are mediated solely by tunneling propagation through these hydrogen bond groups, i.e., (EtIhy hy - Hhy hy)-1. The remaining effective interactions,

hb b(Et) - hb(hy) b (Et), arise from all combinations of tunneling propagation through the hybrid orbitals of the backbone and side chain groups [:N, CR-Hback, CR-Hsidec] and from combinations involving both these orbitals and the orbitals of the hydrogen bond groups [N-H, Cd:O:]. Figure 4c is a plot of hi j(Et) - h(hy) i j (Et). A comparison between parts c and b of Figure 4 shows clearly that the topology of the longer range matrix elements of hb b(Et) (residues k to k + 3 and k to k + 4) arises predominantly from the “typical” hydrogen bond mediation (as defined in eq 6). The other backbone hanging groups and the side chains contribute to shorter range effective interactions (jk to k + 2). The influence of the direct and effective interactions (Figure 4a,b) on the overall tunneling propagation between b orbitals is shown in Figure 4d. This is a plot of the absolute value of Gb b(Et), where Gb b(Et) ) (EtIb b - [Hb b + hb b(Et)])-1 (eq 4). The dominance of the hydrogen bonds is clearly seen in Figure 4d. The effective interactions in Figure 4b between residues (k, k + 3) and (k, k + 4) give rise to peaks in Gb b(Et) that involve tunneling propagation between these same residues. A detailed analysis of such hydrogen bond effects on Gb b(Et) is not the emphasis of this paper. 4. Exact Treatment of Mutations and Analysis of Interferences Having discussed tunneling propagation briefly for the GLY8 structure (denoted I in Figure 3a), we now consider the effect of the double mutation introduced to GLY8 to create GLN1 GLY6 GLN8 (denoted II in Figure 3b). Our goal is to interpret the tunneling propagation between b orbitals in structure II (Gb(II)b (Et)) in terms of the original propagation between b orbitals in structure I (Gb(I)b(Et)) and the changes in electronic interactions between these orbitals that arise from the I f II mutation. The latter changes are described by a “mutation” matrix that is defined and interpreted in this section. In the following calculations we minimize the effects of interatomic distance changes by keeping the positions of all common atoms of GLY8 and GLN1 GLY6 GLN8 frozen. This restriction is not essential. In general, a mutation affects backbone conformation and the mutation matrix accounts for such conformational effects. We use a common natural-hybridorbital43 basis set to describe the b portion of both structures (the natural b hybrids of GLY8). The tunneling propagation for the two structures is given by (eqs 4, 5) (I) Gb(I)b(Et) ) (EtIb b - [Hb(I)b + hbb (Et)])-1

Gb(II)b (Et) ) (EtIb b - [Hb(II)b + hb(II)b (Et)])-1 (7) where b denotes their common natural-hybrid basis set and Et ) -0.2 au (-5.4 eV). To measure the effect of the I f II mutation on the tunneling propagation, we compute the fractional change in the propagation for all b orbital pairs (i j). The resulting matrix Rb b(Et),

Ri j ≡

(I) |G(II) i j | - |Gi j |

|G(I) ij|

(8)

is the importance value matrix.32 It is used to identify regions of the peptide that are most strongly affected by the mutation. If the mutation enhances the magnitude of the tunneling (I) propagation from orbital j to i, then |G(II) i j | > |Gi j | and Ri j > 0. In the opposite limit, where the magnitude is reduced, (I) |G(II) i j | < |Gi j | and Ri j < 0. When there is no change, Ri j ) 0.


Figure 4 (Continued on page XXX.)




Figure 4. (a) Absolute values of direct Hamiltonian matrix elements Hb b between b orbitals of GLY8. These matrix elements are schematically shown in the drawing on the top left of the 2D plot in part a. The x(y) axis numbers of the 2D plot represent b orbital numbers i(j), and the coloring at position i, j corresponds to the value of |Hi j| in atomic units (color bar). The ordering of the b hybrid orbitals is from the N to the C terminus. The grid lines separate adjacent amino acids, and the ordering of b orbitals within each amino acid is shown in Figure 2a. In this plot the maximum absolute value is =1 au, and all values between 0.06 and 1 au are shown as green. (b) Plot similar to part a for the absolute values of all effective matrix elements hb b between b orbitals i, j (i.e. |hi j(Et)|). (c) Plot of |hi j(Et) - hi(hy) j (Et)| showing all tunneling interactions in hi j(Et) that are not mediated solely by the typical hydrogen bonds [N-H‚‚‚:O:dC]. Comparison of b and c shows that the matrix elements of hb b(Et) connecting residues k to k + 3(4) arise from typical hydrogen bond mediation, defined by h(hy) b b (Et) in eq 6. This mediation is schematically shown in the drawing above part b. (d) Similar plot showing the overall tunneling propagation between b orbitals (Gb b(Et) in eq 4). Only absolute values |Gi j(Et)| are shown. There are peaks in the tunneling propagation between residues k to k + 3 and k to k + 4, arising from the hydrogen bonds.

Only (i j) orbital pairs that are affected by the mutation (Ri j * 0) are of interest to us. Figure 5 is a plot of the negative and positive matrix elements of Rb b. At the residue level, the highest density of changes is seen for propagation between residues 1 and 3-4, 1 and 7-8, and between 3-5 and 7-8. These gross features of Rb b seem to be related to structural changes that arise from the mutation (Figure 3). The mutations GLY1 f GLN1 and GLY8 f GLN8 cause the edges of the GLN side chains to come into contact with each other and with GLY3, GLY4, GLY5. It is, therefore, not surprising that the tunneling propagation between residues 1 and 3-4, 1 and 7-8, and 3-5 and 7-8 is changed significantly by the mutation. Another feature of the Rb b(Et) matrix is that there can exist both positive and negative Rb b(Et) elements for orbitals in a given residue pair (e.g., residues 1 to 7-8, 1 to 3-4). This means that the tunneling propagation

between two amino acids can simultaneously be enhanced and reduced (in different parts of the amino acids) as a result of a mutation. For example, consider the longest range propagation (between residues 1 and 8). There is considerable enhancement in the tunneling between orbitals 1 and 48 (R48 1 ) 9.82, (I) |G(II) 48 1| ) 10.82|G48 1|). On the other hand, there is considerable reduction in the tunneling between orbitals 1 and 46 ((R46 1 (I) ) -0.98, |G(II) 46 1| ) 0.02|G46 1|)). Both 48 and 46 belong to the terminal (eighth) residue, indicating that the influence of the mutation on the tunneling is inhomogeneous at the amino acid level. Among the b orbital pairs in Figure 9 with |Ri j(Et)| * 0, we analyze three in detail, (48 1), (39 1), and (48 10). These pairs (indicated by arrows on the Rb b matrix of Figure 5 and drawn schematically below it) are chosen because they show large changes in tunneling propagation and also because they involve



Figure 5. Importance value matrix Rb b for the mutation I f II (Figure 3). The upper half shows all b hybrid pairs i j with Ri j < 0 (reduction of the total tunneling propagation between i and j). The lower half shows all i j pairs with Ri j > 0 (enhancement of the total tunneling propagation between i and j). The x(y) axis represents b hybrid orbital numbers measured from the N to the C terminus (the b hybrids of GLY8 are used as a common basis). The hybrids within each residue are ordered as in Figure 2a. Rather than labeling residues along the x(y) axes as in the previous plots, we label them along the diagonal. For the residues that have been altered by the mutation, the names of both the original and final residue are shown. This representation relates primary sequence changes to the tunneling propagation changes. The three hybrid orbital pairs (i j) that are analyzed in detail are shown by arrows on the Rb b plot. Below the plot is a schematic representation of the position of these pairs along the b portion of the backbone and of their ranges |i - j|.



Figure 6. (a) Matrix |Vi j(Et)| describing the total perturbation arising from the mutation I f II (Figure 3). The largest matrix elements between nonadjacent residues couple the first residue to residues 3-5, and residues 3-5 to residues 7-8. There are also large matrix elements at longer ranges (not visible in this plot) that arise from the combined tunneling propagation through the side chains of both GLNs. These are shown in a separate plot labeled “Longest range interactions”, at the bottom right corner of part a. The magnitude scale in this plot is different from the color bar of part a. It shows that the longest range matrix elements couple residues 2-3 to 7-8, and 1 to 7-8. The top right corner of part a is a (I) schematic representation of the dominant interactions in Vb b(Et). (b) Matrix |Hi(II) j - Hi j | describing the perturbation of the direct Hamiltonian matrix elements induced by the mutation. In the present example the are no conformational changes of the backbone atoms, and these matrix (I) elements only describe local hybrid orbital distortions. (c) Perturbation matrix |hi(II) j (Et) - hi j (Et)| for the effective interactions between b orbitals. Comparison with parts a and b shows that all long-range mediation induced by the mutation is described by this matrix.



orbitals distant from each other, at the edges of the R-helix. The tunneling propagation between 1 and 48 represents the longest range propagation among b orbitals (|i - j| ) 47) and, as discussed above, it shows 1 order of magnitude enhancement. Both (39 1) and (48 10) have identical ranges, (|i - j| ) 38). For (39 1), the mutation reduces the tunneling propagation by (I) 1 order of magnitude (R39 1 ) -0.84, |G(II) 39 1| ) 0.16|G39 1|). For (48 10), the propagation is reduced by 3 orders of magnitude (I) -3 (|G(II) 48 10| ) 10 |G48 10|). The drastically different effects of the mutation on propagation within these three pairs are interpreted with the mutation matrix. 4.1. The Mutation Matrix. At the reduced Hamiltonian level (i.e., Hb b + hb b(Et) in eqs 4, 7), the perturbation caused by the mutation is described by the difference in the reduced Hamiltonians of the two peptides,

Vb b(Et) ≡ [Hb(II)b - Hb(I)b] + [hb(II)b (Et) - hb(II)b (Et)]

(9)

Gb(II)b (Et) ) Gb(I)b(Et) + Gb(I)b(Et) Tb b(Et) Gb(I)b(Et) (10) where

Tb b(Et) ) (Ib b - Vb b(Et) Gb(I)b(Et))-1 Vb b(Et) Tb b is the scattering matrix associated with Vb b. The term Gb(I)b Tb b Gb(I)b describes the influence of the mutation on Gb(I)b to infinite order in Vb b. Namely, for each (i j) pair of b hybrids, we have (I) G(II) i j ) Gi j +

(I) G(I) ∑ i R T R β Gβ j Rβ

(11)

where

TRβ )

VRβ +

1st order in Vbb

VRmG(I) ∑ mnVnβ + all higher orders mn

(12)

}

}

2nd order in Vbb

(I) (R, β, m, and n are also b hybrids). Each term G(I) i R T R β Gβ j represents the sum of an infinite set of GLN-mediated pathways:48

}

(I) (I) (I) G(I) iR TRβGβj ) GiR VRβGβj + 1st order pathways

(I) (I) G(I) ∑ iR VRmGmnVnβGβj + all higher order pathways mn

(13)

}

where the reference Hamiltonian is that of GLY8: Hb(I)b + hb(I)b(Et). The mutation matrix is Vb b(Et). The absolute value of Vb b(Et) is shown in Figure 6a. The term [Hb(II)b - Hb(I)b] in Vb b(Et) describes changes in the direct matrix elements between the b hybrids (Figure 2b left). In general, these changes are caused both by conformational changes of the backbone and by local distortions of the hybrids’ electron densities/hybridizations. The first effect dominates for long-distance ET, as it reflects differences in secondary/tertiary structure. Figure 6b shows the absolute value of [Hb(II)b - Hb(I)b]. Conformational effects are not seen in this example because all backbone atoms in I and II were held fixed. Therefore, only differences in Hb b matrix elements arising from local hybrid orbital distortions appear. These are intra-amino acid effects and are most prominent for the first and last residues (whose side chains were altered) and for the central residues (located near the edges of the GLN side chains in structure II, Figure 3b). The matrix term [hb(II)b (Et) - hb(I)b(Et)] in Vb b(Et) describes changes in the effective Hamiltonian matrix elements (Figure 2b right). As already discussed, these effective elements (hb b(Et) matrix) describe interactions between b hybrids that are mediated by tunneling propagation through the branched part of the protein chain (backbone hanging groups and side chains). [hb(II)b (Et) - hb(I)b(Et)] arises from the secondary/tertiary structural modifications caused by a mutation, because such modifications shift the positions of the hanging groups and side chains, and from changes of the side chains themselves (primarystructure effects). When there are minimal secondary/tertiary structural modifications, [hb(II)b (Et) - hb(I)b(Et)] describes purely primary-structural effects of the mutation. These effects arise from differences in the electronic structure and packing of the side chains. Figure 6c shows the absolute value of [hb(II)b (Et) - hb(I)b(Et)] for the mutation GLY8 f GLN1 GLY6 GLN8. As there is no conformational relaxation of the backbone and only the first and last side chains are modified, the effective couplings in Figure 6c describe new interactions between b hybrids that arise from tunneling propagation through these modified side chains. A comparison of part c with parts a and b of Figure 6 shows that most of the structure of the mutation matrix Vb b(Et) is contained in [hb(II)b (Et) - hb(I)b(Et)]. The largest long-range matrix elements of [hb(II)b (Et) - hb(I)b(Et)] are between residues 1 and 3-5 and between residues 3-5 and 7-8. These interactions arise because the side chains of GLN1 and GLN8 are close to residues GLY3, GLY4, and GLY5 (Figure 3b). Interactions between residues 1 and 3-5 arise predominantly from tunneling propagation through the GLN1 side chain. Similarly, interac-

tions between residues 3-5 and 7-8 arise predominantly from tunneling propagation through the GLN8 side chain. Also shown in Figure 6 are the longest range matrix elements of the mutation matrix. There is a set of matrix elements between residues 2-3 and 7-8 and another set between residues 1 and 7-8. These sets describe interactions that arise from tunneling propagation through both GLN side chains. Such propagation is possible because the edges of the side chains are close to each other (Figure 3b). A comparison of the long-range matrix elements in Figure 6a and Figure 5 shows that many of the peaks in the mutation and importance value matrices occur for the same residue pairs. This means that the mutation matrix provides a useful intermediate description of the mutation. Its relation to the the importance value matrix (i.e. the change in tunneling propagation) is described in the following section. 4.2. Infinite Order Perturbation Theory in the Mutation Matrix. To understand fully changes in the tunneling propagation caused by the mutation (i.e., Gb(I)b(Et) f Gb(II)b (Et)), we must relate the mutation matrix (Vb b(Et)) and the original propagation (Gb(I)b(Et)) to the final propagation (Gb(II)b (Et)). In this section we discuss the exact relations among these three quantities. In general47

2nd order pathways

All pathways in this set share the common feature that GLN mediation begins at hybrid orbital β and ends at hybrid orbital R. Figure 7 is a diagrammatic representation of the terms in eq 13. Equations 10-13 are exact. To understand the contributions of the different pathway sets, Tb b in eq 10 is computed and all terms in ∑R β of eq 11 are ordered from maximum to minimum (I) absolute value.48 The G(I) i R TR β Gβ j elements are then added (I) (I) (I) iteratively in this order to Gi j until G(I) i j + ∑R βGi R TR β Gβ j (II) converges to Gi j . This is one of many possible ways of computing eq 11. It emphasizes the relative magnitudes of the different pathway sets and their competition with the original (I) (I) propagation G(I) i j . Once all the pathway sets Gi R TR β Gβ j



Figure 7. Diagrammatic representation of the mutation matrix VR β, the R-helical (original) propagation Gi(I)j , the first- and second-order GLN(I) (I) (I) mediated pathways Gi(I)R VR β Gβ(I)j and Gi(I)RVR mG(I) m nVn βGβ j, and the pathway set Gi R TR β Gβ j. For both first- and second-order pathways, GLN mediation begins at b orbital β and ends at b orbital R. The first-order pathway involves no intermediate traversals between the GLN side chains and the b orbitals of the backbone. The second-order pathway involves two intermediate traversals (arrows at orbitals m and n). Gi(I)R TR β Gβ(I)j is the sum of all such pathways (all possible intermediate traversals). In the diagram of Gi(I)R TR β Gβ(I)j the possibility of multiple intermediate traversals (eq 13) is represented by a double arrow.

required for convergence are identified, the interferences between them (and with the original propagation prior to the (I) (I) mutation) are analyzed. If |G(I) i R TR β Gβ j|/|Gi j | . 1, the mutation creates a set of GLN-mediated pathways (connecting β to R) whose sum exceeds the original propagation G(I) i j . On (I) (I) the other hand, if |G(I) T G |/|G | ≈ 1, the propagation Rβ iR βj ij through this pathway set is of similar magnitude to the original propagation and it may interfere with it. In this latter situation (I) one needs to consider the relative signs of G(I) i R TR β Gβ j and (I) (I) (I) (I) Gi j . When Gi R TR β Gβ j/Gi j > 0, there is constructive interference between the original and GLN-mediated propaga(I) (I) tion. When G(I) i R TR β Gβ j/Gi j < 0, there is destructive interference with the original propagation, and this may lead to substantial reduction of the overall tunneling propagation, (I) G(II) i j . The relative magnitudes and signs of the Gi R TR β (I) (I) Gβ j/Gi j elements also reflect the interferences between the pathway sets themselves. For a given (i j) pair, this type of interference analysis is rationalized in terms of 2D plots of the (I) (I) matrix G(I) i R TR β Gβ j/Gi j , where x(y) ) R(β). The advantage of this representation (discussed in Figure 8) is that it shows the topology of all relevant pathway sets. Examples are given for the (i j) pairs (48 1), (39 1), and (48 10), shown in Figure 5. For the orbital pair (48 1) involving the eighth (last) and first amino acids, tunneling propagation is enhanced by 1 order of (I) magnitude, with R48 1 ) 9.82 (G(II) 48 1 ) 10.82|G48 1|). The sum (I) (I) (I) G48 1 + ∑R βG48 R TR β Gβ 1 nearly convergences to G(II) 48 1 after

the addition of the 20 dominant sets of GLN-mediated pathways, (I) (I) (I) (II) i.e., |G(II) 48 1 - (G48 1 + ∑RβG48 R TR β Gβ 1)| < 0.1|G48 1| if more than these 20 sets are included in the sum. Figure 9a is a (I) histogram of these 20 pathway sets G(I) 48 R TR β Gβ 1, ordered from maximum to minimum absolute value. Figure 9b shows (I) (I) G(I) 48 1 + ∑R βG48 R TR β Gβ 1, as the sets in Figure 9a are added iteratively in the same order. Superimposed on both plots are (II) the values of G(I) 48 1 and G48 1. Figure 9c is a plot of the negative part (upper half) and positive part (lower half) of the (I) (I) matrix G(I) 48 R TR β Gβ 1/G48 1. This matrix is built from the 20 sets in Figure 9a,b. All sets in the negative part have the opposite sign of the original (purely R-helical) propagation G(I) 48 1 and interfere destructively with it. The sets in the positive part interfere constructively with the original propagation but interfere destructively with the remaining sets (in the negative part). (I) Figure 9a,c shows that the dominant set, G(I) 48 47 T47 2 G2 1, is more than 1 order of magnitude larger than and of opposite (I) (I) (I) sign of the original propagation G(I) 48 1 (G48 47 T47 2 G2 1/G48 1 ) -15). The dominant set accounts for most of the propagation, (I) (I) (II) with G(I) 48 1 + G48 47 T47 2 G2 1 ) 1.3G48 1 (Figure 9b). The diagrammatic representation of this set is shown in Figure 9c. The remaining dominant pathway sets are of varying sign. Among them, the three strongest have magnitudes 3-4 times (I) larger than |G(I) 48 1| and have topologies similar to G48 47 T47 2 (I) G2 1, i.e., GLN side chain mediation begins in the first residue and ends in the last residue. As an example, the diagram of



Figure 8. Explanation of the 2D representation of the ratios Gi(I)R TR β Gβ(I)j/Gi(I)j for a given (i j) pair. Suppose that j is a b orbital of residue 1 and i is a b orbital of residue 8. Consider an R, β pair where β belongs to residue 2 and R to residue 7. The top of the figure shows diagrams for Gi(I)R TR β Gβ(I)j and Gi(I)j . The only difference from Figure 7 is that the residue regions of the b orbital space are shown explicitly as boxes. The bottom of the figure shows the 2D representation of Gi(I)R TR β Gβ(I)j/Gi(I)j . Residue identity and (i j) orbital position are shown along the diagonal. The x(y) axes represent b orbital numbers R(β). The color at position x ) R, y ) β corresponds to the value of the ratio. Grid lines separate adjacent residues. For a given (i j) pair, all ratios Gi(I)R TR β Gβ(I)j/Gi(I)j of the contributing pathway sets are shown on the same 2D plot. This allows the visualization of multiple-pathway sets. For example, all pathway sets for which GLN mediation begins at residue 2 and ends at residue 7 will appear within the box of the 2D plot containing the dot. (I) the second largest pathway set (G(I) 48 45 T45 2 G2 1) is shown in Figure 9c. Although the stepwise addition of these remaining (I) pathway sets causes fluctuations of G(I) 48 1 + ∑R βG48 R TR β (I) (II) Gβ 1 about G48 1 (Figure 9b), these fluctuations are relatively small ( β. Therefore, to compress information, the upper half of this matrix (I) (I) (I) (I) (I) (I) shows all negative G48 R TR β Gβ 1/G48 1, and the lower half shows all positive G48 R TR β Gβ 1/G48 1. The grid lines separate adjacent residues, and residue identity is written along the diagonal together with the the (i j) orbital numbers (written on the residues they belong to). Mutated residues (I) (I) have the names of the original and final amino acids. Also shown in part c are diagrams (Figures 7, 8) for the ratios G(I) 48 47 T47 2 G2 1/G48 1 ) -15 (I) (I) (I) (I) (I) (I) (I) and G48 45 T45 2 G2 1/G48 1 ) 4. G48 47 T47 2 G2 1 and G48 45 T45 2 G2 1 are the two strongest pathway sets. They have very similar topologies (with GLN side chain mediation starting at the first amino acid and ending at the eighth amino acid), but opposite signs.



(I) (I) Figure 10. (a, b) Plots (similar to Figure 9a,b) for the orbital pair (i j) ) (39 1). The 40 dominant G39 R TR β Gβ 1 needed for convergence to (II) (I) (II) 0.1G39 1 are shown. The arrow joining G39 1 to G39 1 in parts a and b indicates a reduction of the tunneling propagation. (c) Ratio plot of the 40 (I) (I) (I) (I) (I) G39 R TR β Gβ 1 in parts a and b (similar to Figure 9c). Since R > β for all 40 G39 R TR β Gβ 1, the upper half of part c shows the negative G39 R TR β (I) (I) (I) (I) (I) (I) (I) G(I) /G and the lower half shows the positive G T G /G . Also shown in part c are diagrams for the ratios G T G /G Rβ 40 2 β1 39 1 39 R β1 39 1 39 40 21 39 1 ) (I) (I) (I) (I) (I) (I) (I) (I) (I) (I) -1.13, G39 42 T42 2 G2 1/G39 1 ) 0.7, and G39 21 T21 2 G2 1/G39 1 ) 0.09. G39 40 T40 2 G2 1 and G39 42 T42 2 G2 1 are the two strongest pathway sets. They have similar topologies, with GLN side chain mediation starting at the first amino acid and ending at the seventh amino acid. They strongly interfere with each other and with the R-helical tunneling propagation G(I) 39 1. The cancellations arising from the interference between such sets and (I) make much weaker sets of very different topology significant to the propagation (a, b). The strongest of these sets is G(I) G(I) 39 1 39 21 T21 2 G2 1 (diagram shown). For this set GLN side chain mediation starts at the first amino acid and ends at the fourth amino acid.



seventh. A diagram of the fifth strongest set, G(I) 39 21 T21 2 G(I) , is shown in Figure 10c. We conclude that, in contrast 21 with the previous case, the large reduction in the tunneling propagation arises from destructive interference between the original (R-helical) propagation G(I) 39 1 and several GLN-mediated pathway sets of different topologies. Physically, it is the cancellations (initial fluctuations of the sum in Figure 10b) between the R-helical propagation and a group of pathway sets for which GLN mediation begins at the first amino acid and ends at the seventh to eighth amino acids that make much weaker sets significant for the propagation. For these weaker sets, GLN mediation begins at the first amino acid and ends at the third to fourth amino acids, or it begins at the fourth to fifth amino acids and ends at the seventh. Figure 11 shows the pathway set analysis for the orbital pair (48 10) that involves the eighth (last) and second amino acids. The mutation reduces the tunneling propagation from orbital 10 to 48 by 3 orders of magnitude. The enormous destructive interference between the R-helical and GLN side chain mediation cannot be attributed to a small fraction of pathway sets (as in the previous examples). Rather, 1109 pathway sets (out of (I) (I) a total 2304) are needed to get |G(II) 48 10 - (G48 10 + ∑R βG48 R (I) (II) TR β Gβ 10)| < 0.1|G48 10|. Figure 11a shows these sets, and (I) Figure 11b is a plot of their partial sum, G(I) 48 10 + ∑R βG48 R TR β G(I) β 10. The small 2D plot of Figure 11c shows all R, β (I) pairs (black dots) of the 1109 G(I) 48 R TR β Gβ 10. The larger 2D (I) (I) plot of Figure 11c shows the ratios G48 R TR β G(I) β 10/G48 10. (I) (I) Only the dominant G48 R TR β Gβ 10 are visible, and none of them are much stronger than the purely R-helical propagation G(I) 48 10 (maximum ratios = (1.2). Even for the most dominant pathway sets, the purely R-helical propagation is not separable from the GLN-side chain mediation. For example, consider the three groups of pathway sets (1), (2), and (3) in Figure 11c. Each group contains pathway sets for which GLN side chain mediation begins and ends in common amino acids (diagrams on the left of Figure 11c). However, the beginning and ending amino acid pairs differ among these groups, spanning the whole peptide. All three groups contain pathway sets of comparable magnitudes that destructively interfere with each other. The strong cancellations between the large number of dominant pathway sets make several hundred much weaker sets relevant to the propagation (Figure 11a,b). 4.3. First-Order Perturbation Theory in the Mutation Matrix. To this point, our analysis has been in terms of pathway sets. (I) For a given (i j) pair, the set G(I) i R TR β Gβ j represents the sum of an infinite number of GLN-mediated pathways G(I) i R VR m (I) G(I) ... V G , all sharing the common feature that GLN nβ mn βj mediation begins at b orbital β and ends at b orbital R (eq 13 and Figure 7). There can be many intermediate traversals between the side chain and b portions of the peptide (Figure 7). The most direct pathway belonging to the set G(I) i R TR β (I) Gβ(I)j is given by G(I) i R VR β Gβ j. This pathway makes no intermediate traversals. For the three pairs of orbitals considered previously, we find that

G(I) 48 1 +

(I) (II) G(I) ∑ 48 R VR β Gβ 1 ) 0.9G48 1 Rβ

G(I) 39 1 +

(I) (II) G(I) ∑ 39 R VR β Gβ 1 ) 1.2G39 1 Rβ

G(I) 48 10 +

(14)

(I) (II) G(I) ∑ 48 R VR β Gβ 10 ) 23G48 10 Rβ

where the sum is over all 2304 matrix elements of Vb b (the

(I) summation over all G(I) i R TR β Gβ j always converges to the (I) exact Gi j ). Therefore, first-order perturbation theory in the mutation matrix can account for most of the change in the tunneling propagation between 1 and 48 and between 1 and 39. First-order perturbation theory overestimates G(II) 48 10 by 1 (II) order of magnitude, but considering that |G48 10| ) 10-3 |G(I) 48 10|, it still predicts 2 out of the 3 orders of magnitude reduction in the tunneling propagation. An analysis of the (I) dominant G(I) i R VR β Gβ j elements for the three orbital pairs shows that they have the same topologies as the dominant (I) G(I) i R TR β Gβ j elements (i.e., the same R, β pairs in Figures 9c and 10c and in the large plot of Figure 11c). They also have (I) (I) (I) similar values (1.2 > (G(I) i R VR β Gβ j)/(Gi R TR β Gβ j) > 0.8). In general, we find that 80% of the 2304 TR β matrix elements are well approximated by the corresponding VR β matrix elements with errors less than 20% (Figure 12). This is because each VR β is infinite order in the interactions within and between the modified side chains, as it involves their full Green functions. For these reasons, the initial rate of convergence of G(I) i j + ∑R β (I) (I) (I) G(I) V G is similar to that of G + ∑ G T Gβ(I)j (i.e., R β R β R β iR βj ij iR similar to Figures 9b, 10b, and 11b). The increasing failure of first-order perturbation theory for the orbital pairs (39 1) and (48 10) (eq 14) arises because the direct GLN-mediated (I) pathways G(I) i R VR β Gβ j interfere with each other or with the original R-helical propagation G(I) i j . The resulting cancellations make weaker higher order pathways with intermediate traversals relevant. In summary, both infinite-order and first-order analyses indicate that drastically different types of mediation are in effect for the three pairs of orbitals shown in Figure 5. For (48 1), the mutation creates new shortcuts for the tunneling electron through the two GLN side chains. One of these shortcuts dominates, having a magnitude ∼15 times larger than the purely R-helical mediation (Figure 9). It surpasses the strength of the R-helical mediation, and it is the source of the order of magnitude increase in the tunneling propagation between 1 and 48. For this shortcut, GLN mediation begins in the first residue (where orbital 1 is located) and ends in the last residue (where orbital 48 is located). The other contributing shortcuts are of similar topology (i.e., first to last residue), but they are weaker than the dominant one (by one-third or less). The contributions of these weak shortcuts cause small fluctuations about the value of tunneling propagation that is determined by the dominant shortcut. Altogether, multiple intermediate traversals of the electron between the GLN side chains and the R-helix are insignificant; a small fraction of all possible pathway sets can accurately reproduce the final tunneling propagation (20/2304). For (39 1) no shortcuts are created. Rather a few alternatiVe routes are created through the GLN side chains that compete in magnitude with the purely R-helical propagation (Figure 10). There is substantial destructive interference between the purely R-helical propagation from 1 to 39 and a set of routes that enter the GLN side chains in the vicinity of the first residue (where orbital 1 is located) and leave them in the vicinity of the seventh residue (where orbital 39 is located). Qualitatively, the order of magnitude reduction in the tunneling propagation arises from cancellations between this latter set of routes and the R-helical propagation. These cancellations make several weaker (by 1 order of magnitude) GLN-mediated routes significant. Such routes enter the GLN side chains in the vicinity of the first residue but leave them in the vicinity of the central residues (third and fourth) rather than the seventh. In this case we begin to see further cross-linking of the electron propagation between the GLN side chains and the remaining R-helix (Figure 10). Cross propagation between the GLN side chains and the R-helix



(I) (I) Figure 11. (a, b) Plots similar to Figures 9a,b and 10a,b for the orbital pair (i j) ) (48 10). The 1109 dominant pathway sets G48 R TR β Gβ 10 needed for convergence to 0.1G(II) are shown. The ticks in the common x axis of parts a and b separate groups of 10 pathway sets. The arrows 48 10 (I) (I) indicate a very large reduction in the tunneling propagation. (c) The small 2D plot on the left shows all R, β pairs of the 1109 G48 R TR β Gβ 10 in (I) (I) (I) parts a and b. Each R, β pair corresponds to a black dot. The large 2D plot on the right shows all 1109 ratios G48 R TR β Gβ 10/G48 10 (similar to parts 9c and 10c). Only the largest ratios are visible, all of which have R g β. The three diagrams shown describe three groups of dominant pathway sets. All pathway sets in a group have the feature that GLN mediation begins and ends in common residues. For this reason these diagrams do not display specific R, β pairs but rather residue pairs (boxes with outgoing and incoming arrows). Diagram (1) describes a group of pathway sets for which GLN mediation begins at GLY2 and ends at GLN8. In diagram (2) GLN mediation starts at GLY5 and ends at GLN8. In (3) GLN mediation starts at GLY2 and ends at GLY5. The large number of dominant pathway sets interfere destructively (part b). The resulting reduction of the propagation makes several hundred very weak pathway sets significant. The topologies of these sets indicate that GLN and R-helical mediation are not separable.



Figure 12. For 80% of a total 2304 TR β matrix elements, VR β is a good approximation to TR β (1.2 > VR β/TR β > 0.8).

(in the middle of the structure) begins to make a noticeable contribution. A larger fraction of pathway sets is needed to reproduce the final propagation (40/2304). For the third pair (48 10), the 3 orders of magnitude reduction in the tunneling propagation arises from 1109 out of a total 2304 GLN-mediated pathways sets (Figure 11). There is no simple interpretation of this enormous destructive interference. The key observation is that the mutation created a large number of alternative routes through the GLN side chains (a far larger number than in the previous case). These routes compete in magnitude with the purely R-helical propagation and with each other. They also have varying topologies. For many of them, GLN mediation does not begin and end at the residues where the bridge donor and acceptor orbitals are located (orbitals 10 and 48). The large cancellations among these routes and the R-helical propagation reduce the tunneling mediation so much as to make hundreds of very weak routes significant. This represents the extreme opposite of the first case, as the R-helical and GLN mediations are not separable. It is interesting to note that the three drastically different kinds of behavior are observed for (i j) pairs of similar ranges and locations along the backbone, as shown in Figure 5. For real ET systems we do not expect dramatic cancellations between pathways (as in Figure 11). Time dependent fluctuations of the backbone and side chains will certainly wash out many interferences. However, when the fluctuations do not significantly alter the local folded conformation of the system (e.g., in the interior of a protein), partial cancellations between dominant pathways may survive. Different aspects of the effects of bridge fluctuations on ET have been investigated,51-53 and recently bridge fluctuations were incorporated in TAD computations.54 However, their specific effects on pathway interferences remain unknown. 5. Discussion We have presented a method to compare electron tunneling interactions in different proteins. To rationalize the origin of electronic propagation differences in distinct protein structures and to provide a systematic means of interpreting the influence of structure on electron transfer, we have introduced electron transfer contact maps. These maps are motivated by distance contact maps used to analyze protein conformations and by other reduced descriptions of proteins, such as Ramachandran plots and ribbon drawings. Electron transfer contact maps use all orbitals of the unbranched protein backbone ([-N-CR-C-]k) as the reference orbital set. This reduction allows exact comparisons of electronic propagation in proteins that have different amino acid sequences and folded conformations. It also allows the detailed interpretation of how mutations influence electron transfer.

Let us now summarize this strategy of protein analysis. Suppose two proteins are to be compared. First, the projection method is used to construct reduced Hamiltonians for the unbranched backbone of each protein. To accomplish this, a common hybrid orbital basis is chosen. If the proteins have different lengths, the projections are on the unbranched portion of the backbone that is common to the two proteins. This projection involves no approximation because the influence of the remaining backbone in each protein will appear in the protein’s reduced Hamiltonian as a set of effective matrix elements. In the next step, one of the proteins is used as a reference and the mutation matrix associated with the differences (arising from sequence or conformation) between the two protein structures is computed. The mutation matrix is the difference between the reduced Hamiltonians of the two proteins. The differences in tunneling propagations (Green function matrices) between the proteins is interpreted in terms of this mutation matrix and the tunneling propagation of the reference structure. If necessary, this analysis can be done via infinite order perturbation theory in the mutation matrix by constructing its scattering matrix (perturbations arising from large changes in protein structure will be strong). The reduced Hamiltonian matrices, the mutation matrix, and the subsequent scattering matrix analysis of the tunneling are visualized in terms of 2D plots that keep track of the common backbone orbital and residue identities. We call such plots electron transfer contact maps. These maps show the topologies and strengths of tunneling interactions associated with changes in primary sequence and/or folded conformation. In cases where the two proteins have very different structure (e.g., R-helix versus β-sheet), it is more useful to compare the reduced Hamiltonians themselves, rather than to compute differences in reduced Hamiltonians. This comparison is made by plotting them against each other in upper and lower triangular matrices, with the residue identities of the common backbone orbitals indicated along the diagonal separating these matrices (i.e., similar to Figures 1a, 5, and 9c-11c). The scattering matrix analysis of the reduced Hamiltonians is also done separately for each protein and compared in the same way. For example, rather than computing the mutation matrix for GLY8 and GLN1 GLY6 GLN8, we could have analyzed the reduced Hamiltonian and resulting pathway sets in each structure separately and then compared them to each other. This would have given the same physical picture as the mutation matrix analysis. The proposed method is readily generalized to probe larger portions of the proteins under comparison. If two proteins have a large fraction of identical residues, the atoms of these residues may also be included in the reduced set of atoms together with the backbone. For example, in GLY8 and GLN1 GLY6 GLN8,

1232 J. Phys. Chem. B, Vol. 101, No. 7, 1997 only the first and last side chains are different. Rather than projecting on the unbranched backbone of each peptide, we could have projected on the set of all backbone orbitals and all hydrogen side chains of the common GLY2-GLY6 section of the peptides. We would then probe a larger set of tunneling propagation changes. The disadvantage of a large projection space is that the reduced Hamiltonians are not as readily interpretable (for example, if the projection space includes several large side chains, the reference structure is a highly branched structure). On the other hand, a smaller projection space (e.g., only the CR’s) may constitute a very small portion of each protein under comparison. Another obvious generalization of the method is in the definition of the mutation matrix. A “mutation” does not necessarily represent modifications of the protein primary sequence. It may also describe insertions of water molecules or cofactors in a protein. In addition to studying global changes of the tunneling propagation (i.e., large sections of the G(Et) matrix), one can also apply this method to the study of donor-acceptor matrix elements. In such cases, the G(Et) elements of interest are those that enter TD A ) ∑a dvA a Ga d(Et) vd D (eq 2). Suppose that a set of mutations is introduced to a protein I (connected to D (A)) to produce protein II. The change in TD A is analyzed in terms of the individual elements vA a G(I) a d(Et) vd D and vA a′ (I) (II) G(I) (E ) v of T and T , as in section 4. The only new t d′ D a′ d′ DA DA feature of the analysis is the study of interferences between the elements vA a Ga d(Et) vd D in each TD A expression. Destructive interferences between these elements may make each TD A smaller than any single vA a Ga d(Et) vd D contribution. Also, very large donors and acceptors that couple to different parts of each protein can wash out differences in tunneling propagation that are seen through an analysis of G(Et). The interpretation of interferences between vA a Ga d(Et) vd D elements can be performed in terms of partial sums, as in Figures 9a-11a. Although our analysis of the G(Et) matrix is nonperturbative, the formula TD A ) ∑a dvA a Ga d(Et) vd D is perturbative in vA a vd D. However, as long as the two-state limit is valid and the electron transfer is long-distance, this formula can reproduce the exact energy splitting by variation of Et.37 For example, ˜D ≡ consider a symmetric system with ED ) EA. We define E ˜ A ) ∑a1 a2vA a1 Ga1 a2(EA) ED + ∑d1 d2vD d1 Gd1 d2(ED) vd2 D and E va2 A, where d1, d2 (a1, a2) are orbitals of the bridge that couple directly to |D〉 (|A〉). Rather than setting Et ) ED ) EA in ˜D ) E ˜ A.55 This takes into account Ga d(Et), we can choose Et ) E the effect of backscattering between the donor (acceptor) and the bridge that can be significant even for long-distance ET. 6. Conclusions A new strategy for structure-function analysis in protein electron transfer was described based on the idea of electron transfer contact maps. A simple application of this strategy was presented for two small peptides using CNDO semiempirical Hamiltonians. It was shown that completely different tunneling mediation regimes can coexist even in such small systems. These regimes range from the single-pathway to the multiple-pathway limits. A mechanism for producing multipathway behavior in the presence of a few dominant single pathways was identified. When a few dominant pathways of similar magnitude but opposite sign destructively interfere, their sum is comparable in magnitude to the contribution from a multitude of weaker pathways. In this situation, the collective contribution of the weaker pathways to the matrix element is important. Recent large-scale computations (∼3500 orbitals) using the CNDO/S method56 suggest that contact map analysis is computationally feasible (at least for the portions of proteins relevant to ET). Electron transfer contact maps describe global features of

Skourtis and Beratan tunneling mediation, and they should be useful for studying the transferability of tunneling interactions arising from primary, secondary, and tertiary structure. Acknowledgment. We thank Michael Grabe for assistance with the programming, and Professors S. Priyadarshy and F. Weinhold for guidance with the CNDO/S method and the BONDO program. We thank NIH (GM48043) and NSF (CHE9257093) for support of this research. Appendix Here we derive generalizations of all formulas used in the main text. The generalized case developed here applies to a nonorthogonal basis set. For a nonorthogonal basis set the secular equation of the entire system (D, A, and protein) is

HC B ) ESC B

(15)

H is the total Hamiltonian, C B is an eigenvector of eigenvalue E, and S is the total overlap matrix. All derivations that follow involve partitioning the Hilbert space of the full system into two subspaces. Therefore, we first consider the general problem of partitioning.36 Suppose M is B is a column vector a square matrix (M-1 is its inverse) and C in the same vector space. The vector space is partitioned into two subspaces p and q (p + q ) I), such that

(

)

(

)

M M [M-1]p p [M-1]p q , M ) M p p Mp q , M-1 ) [M-1]q p [M-1]q q qp qq

()

C Bp C B) C Bq

(16)

B) where Mp p, Mp q, etc., are submatrices. If the equation MC 0 is written in its partitioned form and then solved for the p subspace, one finds

(Mp p - Mp q(Mq q)-1 Mq p)C Bp ) 0

(17)

Similarly, writing M M-1 ) M-1 M ) I in its partitioned form and solving for [M-1]p p, leads to49

[M-1]p p ) (Mp p - Mp q[Mq q]-1 Mq p)-1

(18)

Derivation of TD A. First we present the derivation of TD A for a nonorthogonal basis set. The aim is to approximate the two eigenenergies of the entire system (D, A, and protein) that arise from the eigenstates with most dominant |D〉 and |A〉 character. This may be done by an iterative procedure23,36 that diagonalizes a reduced 2 × 2 Hamiltonian in the D, A subspace. The limitations of the two-state reduction are discussed in refs 37 and 38. The secular equation (eq 15) for the entire system (D, A, and protein) is written (ES - H)C B ) 0. p is taken to be the subspace of the |D〉 and |A〉 states (denoted DA), and q is taken to be the subspace of all protein states (denoted pr). These states may be localized atomic or bonding/antibonding protein orbitals or any other convenient basis set of protein orbitals. Setting M ) (ES - H), p f DA, and q f pr in eq 17, we find

HDA DA(E)C B DA ) SDA DAC B DA

(19)

where

HDA DA(E) ) HDA DA + (ESDA pr - HDA pr)(ESpr pr Hpr pr)-1(ESpr DA - Hpr DA) HDA DA(E) is a 2 × 2 reduced (or effective) Hamiltonian in the



D, A subspace. Its off-diagonal matrix element is

HA D(E) ) HA D +

Vb b(Et) ≡ Et[S b(I)b - S b(II)b ] + [H b(II)b - H b(I)b] +

(ESA a - HA a)[(ESpr pr ∑ ad

[H˜ b(II)b (Et) - H˜ b(I)b(Et)] (25)

Hpr pr)-1]a d(ESd D - Hd D) (20) The first step in the iterative procedure is to set E in HDA DA(E) equal to the zeroth-order D and A energies (E ) E(1) ) ED, B DA ) SDA DAC B DA. This EA) and then to diagonalize HDA DA(E(1))C gives two eigenenergies E+(E(1)), E-(E(1)) that are a first approximation to the two relevant eigenenergies of the entire system. In the following step one can set E ) E(2) ) (E+(E(1)) + E-(E(1)))/2, and similar iterations are repeated until E+, Edo not change with subsequent iterations. We call the final E(n) the tunneling energy, Et. The converged reduced Hamiltonian is HDA DA(Et), where HDA DA(E) is given in eq 19. TD A is taken as the off-diagonal (A D) element of HDA DA(Et). For long-distance ET where HA D is negligible, we find from eq 20

TD A ≡ HA D(Et) )

G b(II)b (Et) ) G b(I)b(Et) + G b(I)b(Et)tb b(Et)G b(I)b(Et) (26) where

tb b(Et) ) (Ib b - Vb b(Et)G b(I)b(Et))-1Vb b(Et) tb b(Et) is the scattering matrix associated with the mutation. A similar analysis (starting with eq 18) can be used to obtained a generalization of the hydrogen bond coupling matrix hb(hy) b (Et) (eq 6). It is given by

H˜ b(hy) b (Et) ≡ (Et Sb hy - Hb hy)(Et Shy hy Hhy hy)-1(Et Shy b - Hhy b) (27)

(EtSA a - HA a)[(EtSpr pr ∑ ad Hpr pr)-1]a d(EtSd D - Hd D) (21)

where a, d are protein orbitals that couple directly to the D and A states. For S ) I, eq 21 becomes eq 2, where in eq 2 vA a ≡ HA a, vd D ≡ Hd D, and H ≡ Hpr pr. The effect of the protein on TD A is contained in the matrix (EtSpr pr - Hpr pr)-1. For S ) I this matrix becomes G(Et) in eq 3, where H ≡ Hpr pr. (EtSpr pr - Hpr pr)-1 is not the Green function of the isolated (uncoupled from D and A) protein. The Green function of the protein is given by Spr pr(EtSpr pr - Hpr pr)-1Spr pr. Derivation of the Reduced Hamiltonians, Mutation Matrix, and Scattering Matrix. Equation 21 shows that the effect of the protein on TD A is contained in (EtSpr pr - Hpr pr)-1. Therefore, the derivation of the reduced Hamiltonians, the mutation matrix, and the scattering matrix must start from partitioning (EtSpr pr - Hpr pr)-1. We partition the set pr of all protein orbitals into the subspace of atomic hybrids for the unbranched backbone of the protein (b) and the subspace of hybrids for the remaining, branched portion of the protein (hr). Setting M ) (EtSpr pr - Hpr pr), p f b, and q f hr in eq 18, we find

Gb b(Et) ≡ [(EtSpr pr - Hpr pr)-1]b b ) (EtSb b - [Hb b + H˜b b(Et)])-1 (22) where

H˜b b(Et) ) (EtSb hr - Hb hr)(EtShr hr - Hhr hr)-1(EtShr b Hhr b) (23) Equations 22 and 23 become eqs 4 and 5 for S ) I where, in eqs 4 and 5, H ≡ Hpr pr. Suppose two proteins (denoted I and II) are compared. The total tunneling propagation between b orbitals in each of the two proteins is given by

G b(I)b(Et) ) (EtS b(I)b - [H b(I)b + H˜ b(I)b(Et)])-1

such that

(24)

G b(II)b (Et) ) (EtS b(II)b - [H b(II)b + H˜ b(II)b (Et)])-1 where H˜ (I) and H˜ (II) have the form of eq 23 with S ) S (I), H ) H (I) and S ) S (II), H ) H (II), respectively. The mutation matrix is

A more general approach to partitioning, beginning with the Green function of the entire system S(EtS - H)-1S rather than its secular equation (eq 15), was recently discussed in ref 50. The two approaches are formally equivalent only when S ) I. For a review of other aspects of nonorthogonality in ET see ref 24. References and Notes (1) Barbara, P. F.; Meyer, T. J.; Ratner, M. A. J. Phys. Chem. 1996, 100, 13148. Bertini, I.; Gray, H. B.; Lipard, S. J.; Valentine, J. S. Bioinorganic Chemistry; University Science Books: Mill Valley, CA, 1994. Skourtis, S. S.; Beratan, D. N. AdV. Chem. Phys., in press. (2) Beratan, D. N.; Onuchic, J. N. In Protein Electron Transfer; Bendall, D. S., Ed.; BIOS Scientific: Oxford, 1996; Chapter 2. Curry, W. B.; Grabe M. D.; Kurnikov, I. V.; Skourtis, S. S.; Beratan, D. N.; Regan, J. J.; Aquino, A. J. A.; Beroza, P.; Onuchic, J. N. J. Bioenerg. Biomembr. 1995, 27, 285. (3) Beratan, D. N.; Onuchic, J. N.; Winkler, J. R.; Gray, H. B. Science 1992, 258, 1740. Betts, J. N.; Beratan, D. N.; Onuchic, J. N. J. Am. Chem. Soc. 1992, 114, 4043. Beratan, D. N.; Betts, J. N.; Onuchic, J. N. J. Phys. Chem. 1992, 96, 2852. Beratan, D. N.; Betts, J. N.; Onuchic, J. N. Science 1991 252, 1285. Onuchic, J. N.; de Andrade, P. C. P.; Beratan, D. N. J. Chem. Phys. 1991, 95, 1131. Beratan, D. N.; Onuchic, J. N.; Betts, J. N.; Bowler, B.; Gray, H. B. J. Am. Chem. Soc. 1990, 112, 7915. Onuchic, J. N.; Beratan, D. N. J. Chem. Phys. 1990, 92, 722. Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488. (4) Regan, J. J.; Di Bilio, A. J.; Langen, R.; Skov, L. K.; Winkler, J. R.; Gray, H. B.; Onuchic, J. N. Chem. Biol. 1995, 2, 489. Regan, J. J.; Risser, S. M.; Beratan, D. N.; Onuchic, J. N. J. Phys. Chem. 1993, 97, 13083. (5) Gruschus, J. M.; Kuki, A. J. Phys. Chem. 1993, 97, 5581. Kuki, A.; Wolynes, P. G. Science 1987, 1647. (6) Gehlen, J. N.; Daizadeh, I.; Stuchebrukhov, A. A.; Marcus, R. A. Inorg. Chim. Acta 1996, 234, 271. Stuchebrukhov, A. A.; Marcus, R. A. J. Phys. Chem. 1995, 99, 7581. Siddarth P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 6112. Siddarth P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 2400. (7) Stuchebrukhov, A. A. J. Chem. Phys. 1996, 104, 8424. (8) Nakagawa, H.; Koyama, Y.; Okada, T. J. Biochem. 1994, 115, 891. (9) Okada, A.; Kakitani, T.; Inoue, J. J. Phys. Chem. 1995, 99, 2946. (10) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 355, 796. (11) Lopez Castillo, J. M.; Filali Mouhim, A.; Plante, I. L.; Jay Gerin, J. P. Chem. Phys. Lett. 1995, 239, 223. Lopez Castillo, J. M.; Filali Mouhim, A.; Plante, I. L.; Jay Gerin, J. P. J. Phys. Chem. 1995, 99, 6864. (12) Bicout, D. J.; Field, M. J. J. Phys. Chem. 1995, 99, 12661. (13) Langen, R.; Chang, I.-J.; Germanas, J. P.; Richards, J. H.; Winkler, J. R.; Gray, H. B. Science 1995, 268, 1733. Karpishin, T. B.; Grinstaff, M. W.; Komar-Panicucci, S.; McLendon, G.; Gray, H. B. Structure 1994, 2, 415. Wuttke, D. S.; Bjerrum, M. J.; Winkler, J. R. Gray, H. B. Science 1992, 256, 1007. (14) Onuchic, J. N.; Beratan, D. N.; Winkler, J. R.; Gray, H. B. Annu. ReV. Biophys. Biomol. Struct. 1992, 21, 349. (15) Nocek, J. M.; Zhou, J. S.; De Forest, S., Priyadarshy, S.; Beratan, D. N.; Onuchic, J. N.; Hoffman, B. M. Chem. ReV. 1996, 96, 2459. (16) Moreira, I.; Sun, J.; Cho, M. O. K.; Wishart, J. F.; Isied, S. S. J. Am. Chem. Soc. 1995, 116, 8396. (17) Ullmann, G. M.; Kostic, N. M. J. Am. Chem. Soc. 1995, 116, 8396.

1234 J. Phys. Chem. B, Vol. 101, No. 7, 1997 (18) Kurnikov, I. V.; Beratan, D. N. J. Chem. Phys. 1996, 105, 9561. (19) Marcus, R. A. J. Chem. Phys. 1956, 24, 979. Marcus, R. A. J. Chem. Phys. 1965, 43, 679. (20) Hush, N. S. Trans. Faraday Soc. 1961, 57, 155. Hush, N. S. Electrochim. Acta 1968, 13, 1005. (21) Hopfield, J. J. Proc. Natl. Acad. Sci. (U.S.A.) 1974, 71, 3640. (22) Bixon, M.; Jortner, J. J. Chem. Phys. 1968, 48, 715. Jortner, J.; J. Chem. Phys. 1976, 64, 4860. (23) Larsson, S. J. Am. Chem. Soc. 1981, 103, 4034. Larsson, S. J. Chem. Soc., Faraday Trans. 2 1983, 79, 1375. (24) Newton, M. D. Chem. ReV. 1991, 91, 767. (25) Liang, C.; Newton, M. D. J. Phys. Chem. 1992, 96, 2855. (26) Jordan, K. D.; Paddon-Row, M. N. Chem. ReV. 1992, 92, 395; J. Phys. Chem. 1992, 96, 1188. (27) Curtiss, L. A.; Naleway, C. A.; Miller, J. R. Chem. Phys. 1993, 176, 387. (28) da Gama, A. A. S. J. Theor. Biol. 1990, 142, 251. (29) Ratner, M. A. J. Phys. Chem. 1990, 94, 4887. (30) Goldman, C. Phys. ReV. A 1991, 43, 4500. (31) Malinsky, J.; Magarshak, Y. J. Phys. Chem. 1992, 92, 2849. (32) Skourtis, S. S.; Regan, J. J.; Onuchic, J. N. J. Phys. Chem. 1994, 98, 3379. (33) Creighton, T. E. Proteins. Structures and Molecular Properties; W. H. Freeman: New York, 1984. (34) Branden, C.; Tooze, J. Introduction to Protein Structure; Garland Publishing Inc: New York, London, 1991. (35) Ramachandran, G. N.; Saisekharan, V. AdV. Protein Chem. 1968, 23, 283. (36) Lo¨wdin, PO. J. Math. Phys. 1962, 3, 969; Lo¨wdin, PO. J. Mol. Spectrosc. 1964, 13, 326. (37) Skourtis, S. S.; Mukamel, S. S. Chem. Phys. 1995, 197, 367. Skourtis, S. S.; Beratan, D. N.; Onuchic, J. N. Chem. Phys. 1993, 176, 501. Skourtis, S. S.; Onuchic, J. N. Chem. Phys. Lett. 1993, 209, 171. (38) Mujica, V.; Kemp, M.; Ratner, M. A. J. Chem. Phys. 1994, 8, 6849. Mujica, V.; Kemp, M.; Ratner, M. A. J. Chem. Phys. 1994, 8, 6856.

Skourtis and Beratan (39) Coutinho Neto, M. D.; da Gama, A. A. D. Chem. Phys. 1996, 203, 43. Coutinho Neto, M. D.; da Gama, A. A. D. J. Mol. Struct. (THEOCHEM) 1995, 330, 437. (40) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford University Press: New York, 1995. (41) Pople, J. A.; Beveridge, D. L. Approximate Molecular Orbital Theory; McGraw-Hill Book Co.: New York, 1970. (42) Del Bene, J.; Jaffe, H. H. J. Chem. Phys. 1968, 48, 1807. (43) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899. (44) Weinhold, F. QCPE Program No. 408. (45) Dahiyat, B. I.; Meade, T. J.; Mayo, S. L. Inorg. Chim. Acta 1996, 234, 207. (46) Rabanel, F.; Gibney, B. R.; DeGrado, W. F.; Moser, C. C.; Dutton, P. L. Inorg. Chim. Acta 1996, 234, 213. (47) Economou, E. N. Green's Functions in Quantum Physics, 2nd ed.; Springer-Verlag: Berlin, 1990. (48) Skourtis, S. S.; Onuchic, J. N.; Beratan, D. N. Inorg. Chim. Acta 1996, 234, 167. (49) Ayres, F., Jr. Theory and Problems of Matrices; Schaum Publishing Co.: New York, 1962. (50) Priyadarshy, S.; Skourtis, S. S.; Risser, S. M.; Beratan, D. N. J. Chem. Phys. 1996, 104, 9473. (51) Onuchic, J. N.; da Gamma, A. A. S. Theor. Chim. Acta 1986, 69, 89. (52) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488. (53) Tang, J. J. Chem. Phys. 1993, 98, 6263. (54) Wolfgang, J.; Risser, S. M.; Priyadarshy, S.; Beratan, D. N. J. Phys. Chem., in press. (55) Davydov, A. S.; Gadidei, Y. B. Phys. Status Solidi B 1985, 132, 189. (56) Priyadarshy, S.; Risser, S. M.; Beratan, D. N. Int. J. Quantum. Chem., Quantum Biol. 1996, 60, 65. Priyadarshy, S.; Risser, S. M.; Beratan, D. N. J. Phys. Chem. 1996, 100, 17678.

Electron Transfer Contact Maps - American Chemical Society

Recommend Documents