Worm Model for Electron Tunneling in Proteins: Consolidation of the

4424

J. Phys. Chem. B 2001, 105, 4424-4435

Worm Model for Electron Tunneling in Proteins: Consolidation of the Pathway Model and the Dutton Plot Tsutomu Kawatsu, Toshiaki Kakitani,* and Takahisa Yamato Department of Physics, Graduate School of Science, Nagoya UniVersity, Chikusa-ku, Nagoya 464-8602, Japan ReceiVed: October 25, 2000; In Final Form: March 2, 2001

We present a novel method for selecting important electron tunneling pathways in proteins by connecting important interatomic tunneling currents. Then, we constituted an electron tunneling route, called a “worm”, which is formed by averaging over all of the interatomic tunneling currents. The method is applied to six kinds of Ru-modified azurins where the electron transfer takes place from the copper ion to the ruthenium ion linked to the surface of the azurin. We found that the worm is not straight but winds, reflecting the specific role of the microenvironment of the protein structure and ligands in the electron tunneling pathway and rather narrow radius of the worm, ca. 1.5 Å at most. The donor-acceptor distance dependence of the electron-transfer rate ke is examined by taking into account the one-dimensional atomic density in the worm. We found that the Dutton plot, in which the logarithm of ke is in a linear relation with the donor-acceptor distance, which, in turn, depends somewhat on the one-dimensional atom density, is reproduced for a certain level of analysis by the worm model. Thus, we can consolidate the electron tunneling pathway model with the Dutton plot by using the worm model.

1. Introduction The long-range electron transfer between a donor (D) and an acceptor (A) embedded in a protein environment takes place by the superexchange mechanism. The rate of such electron transfer is expressed as

ke )

2π |T |2(FC) p DA

(1)

where TDA is the electron tunneling matrix element and (FC) is the Franck-Condon factor. In the superexchange mechanism, the protein environment plays an important role in the TDA term as a mediator of the virtual state. The protein environment is composed of various arrangements of secondary structures such as R helixes, β strands, and random coils, as well as hydrogenbond and van der Waals contacts, surrounded by lipids or water molecules. Therefore, determining which part of the protein environment is most significantly used in the electron tunneling between a redox pair (donor and acceptor) is an intriguing. Up to the present time, two different models for this problem have been proposed. A model suggested by Beratan and Onuchic is that specific electron tunneling pathways exist in which chemical bonds are predominantly used, hydrogen bonds are subsidiarily used, and through-space routes are rarely used.1-5 According to this model, the electron tunneling pathways depend on the detail of the 3D structure of the protein. The theoretically calculated overall shapes of the electron tunneling pathways in this model are far from straight, but rather branch and assume a significant breadth.1-5 This model is called the “pathway” model1,2,5 or the “tube” model in its advanced version.3,4 Experimental support for this model was provided by Gray and co-workers.6,9 They devised Ru-modified cytochrome c6,7,9 and Ru-modified azurin,8,9 so that each Ru complex could function as an acceptor at anarbitrarily introduced position. They * Corresponding author. E-mail: [email protected].

observed that there is a better correlation between the electrontransfer rate and the effective pathway length than between the electron-transfer rate and the direct distance between the donor and acceptor in the Ru-modified cytochrome c systems.6,7 The other model assumes that the electron tunneling pathway can be approximately expressed by the straight line with a width defined by the sizes of the connecting donor and acceptor. In particular, Dutton and co-workers plotted the logarithm of the rate (log10 ke) as a function of the distance R between donor and acceptor taken from the various biological systems.10 They found that an almost linear relationship holds true for measured optional electron-transfer rates ranging from 1ps to 1s (the slope β ≈ 1.4 Å-1), with a deviation of at most one order of the rate from this straight line.10 Later, they improved the correlation diagram by modifying the slope β so that it could depend on the packing density of atoms F in the protein environment confined by the straight pipe connecting the donor and acceptor molecules.11 Such correlation diagrams between log ke and R for various redox pairs in proteins are called Dutton plots. According to the concept of Dutton plots, the electron transfer rate is determined only by the distance between the redox pair, with specific electron tunneling pathways being insignificant in biological electron transfer. This is a great difference from the pathway model according to whichspecific tunneling paths are expected to play a significant role in biological electron transfer. On the other hand, much theoretical effort has been applied to the determination of the electron tunneling pathways using quantum chemical calculations. Siddarth and Marcus calculated the electron tunneling matrix elements TDA by the extended Hu¨ckel method12 for the important parts of the protein environment, which were selected by the artificial intelligence method,13 and applied them to cytochrome c and myoglobin derivatives.14,15 The advanced method for calculating TDA within the framework of the Hu¨ckel method was devised by Stuchebrukhov et al.16-18 An index for directly searching the important amino

10.1021/jp003918l CCC: $20.00 © 2001 American Chemical Society Published on Web 04/20/2001

Worm Model for Electron Tunneling in Proteins

J. Phys. Chem. B, Vol. 105, No. 19, 2001 4425

acids for electron transfer was implemented by Okada et al.19 Stuchebrukhov and co-workers presented a theory for the quantum mechanical tunneling currents flowing through individual atoms and searched the important atoms and important amino acids in the electron tunneling process.20-23 Recently, we proposed a novel method for determining the electron tunneling pathways based on the calculated interatomic tunneling currents.24 The theory was applied to the electron tunneling pathways in Ru-modified azurins. In that study, we found that the through-space routes play a more important role in the tunneling pathway than previously expected.24 In the present paper, we develop the previous theory using the interatomic tunneling current so as to establish a method for describing an effective tunneling region in a more conceptual way. We show that the shape of such a tunneling region resembles that of a worm. Manipulating ke by the worm model, we show that the relation between log ke and R is quite similar to that obtained by the Dutton group. We show how we can consolidate the Dutton plot and the pathway model by applying the worm model. 2. Theoretical Method We start from theoretical formulation for tunneling currents, devised in previous work.24 Molecular orbitals in the initial and final states ψi and ψf are given as

ψi ) CiDφD +

∑ν Ciνφν

(2)

ψf ) CfAφA +

∑ν Cfνφν

(3)

where φD and φA are orbitals in the donor and acceptor, respectively; φν is the νth atomic orbital in the protein, including ligands; and Ciν and Cfν are molecular orbital coefficients in the initial and final states, respectively. Because the donor and the acceptor weakly couple with the protein environment, |CiD| ≈ 1 and |CfA| ≈ 1 hold. All of the coefficients are determined by solving the secular equation in the extended Hu¨ckel theory. The tunneling current from atomic orbital µ to atomic orbital ν, Jµν, is written as20

Figure 1. (a) Model of the normalized interatomic tunneling currents K (large arrows) through the atoms. The open circles are the mediator atoms (a, b, ...). The filled circles are the donor and acceptor atoms. Sj is the jth perpendicular plane to the line connecting the donor and acceptor. Rjab represents the vector of the point at which the tunneling current Kab crosses to the plane Sj. (b) Worm model of the tunneling pathway. The worm is made by connecting each circle with the neighboring circles. The center and radius of the circle on the jth plane are Rj and σj, where Rj and σj are the barycenter and standard deviation, respectively, of the tunneling currents on the jth perpendicular plane to the line connecting donor and acceptor.

1 Jµν ) -Jνµ ) (Cfµ Ciν - CfνCiµ)(1/G)µν p

The interatomic tunneling current Jab from atom a to atom b is defined as

(4)

where

Jab ≡ (1/G)µν ≡ ESµν - Hµν

∑ ∑Jµν µ∈a ν∈b

(10)

(5) The total tunneling matrix element TDA is given by16-18

E is the energy of the initial and final states, Sµν is the overlap integral, and Hµν is the matrix element of the electronic Hamiltonian. The coefficients C are expressed as24

Ciν ) -CiD

∑µ D†µGµν

(6)

∑µ GνµAµ

(7)

Cfν ) -CfA

TDA ) CiD CfA

D†µGµνAν ∑ µν

(11)

TDA is expressed as20

TDA ) p

∑ ∑ Jab

(12)

i b∉Ωi a∈ΩD D

where

D†ν ≡ ESDν - HDν

(8)

Aν ≡ ESνA - HνA

(9)

where ΩjD represents the region of the donor side as the protein is divided by the plane Sj (see Figure 1a). Equation 12 indicates that the tunneling matrix element TDA is expressed by a sum of the current Jab across the arbitrary plane located between the donor and acceptor. The normalized tunneling current Kab is

4426 J. Phys. Chem. B, Vol. 105, No. 19, 2001

Kawatsu et al.

Figure 2. Electron tunneling currents (arrows) and worm (shaded region) in the His83 derivative. Important atoms that take part in the electron tunneling current with |Kab| > 0.1 are drawn by spheres.25 Only the interatomic currents that have |Kab| > 0.2 are drawn. The width of the arrow is drawn in proportion to the magnitude of the current.

defined as

Jab Kab ≡ p TDA

(13)

Because of the conservation of the total current across an arbitrary plane, we obtain the following sum rule for Kab (see Figure 1a)24

∑ ∑

a∈ΩiD

Kab ) 1

(14)

b∉ΩiD

where Kab is positive when the current flows from the atom in the donor region to the atom in the acceptor region, and it is

negative when the current flows from the atom in the acceptor region to the atom in the donor region. We express the position of barycenter Rj of all the tunneling currents passing through the plane Sj as

∑ ∑

Rj )

a∈ΩiD

|Kab|Rjab

b∉ΩiD

∑ ∑

(15) |Kab|

a∈ΩiD b∉ΩiD

where Rjab is the position where the current Kab passes through the plane Sj (see Figure 1a). Equation 15 is weighted by the



Figure 3. Electron tunneling currents (arrows) and worm (shaded region) in the His124 derivative. The others are as in Figure 2.

absolute value |Kab| rather than Kab because the tunneling current circulates in the broad area of the protein and so the global center of such a current is a property described by considering equally the positive and negative current. This is explained in more detail in the following sections. The standard deviation from the barycenter is given as follows:

σj )

(∑

a∈ΩiD

∑

b∉ΩiD

|Rjab - Rj|2|Kab|

)

1/2

(16)

The above considerations apply to an arbitrary plane Sj. In the following, we position a set of planes {Sj} such that they vertically cut a straight line connecting the donor and acceptor,

as shown in Figure 1b. The central line connecting the barycentric positions Rj in the neighboring planes constitutes a barycentric line of the tunneling area. The outline of the tunneling route is obtained by connecting the circles with centers at the positions Rj and with radii of σj, as shown in Figure 1b. Because the outline of the tunneling region resembles a worm, we call this representation a worm model. In actual calculations, we position parallel planes at intervals of 0.001 Å vertically along the straight line connecting the donor and acceptor. However, it happens that the barycentral line is slightly bumpy because the center of the atomic position is a delta function in the Cartesian coordinate representation. Thus, we take an average for the position of a barycenter Rj over the positions of neighboring planes. In the actual calculation, we


Kawatsu et al.

take average Rj’s over 1001 successive planes. The average is calculated as j+500

R hj )

Rk

∑ k)j-5001001

(17)

where the summation over k is made as far as the plane exists. When less than 1001 successive planes exist at positions near the donor or acceptor, we average the Rj values over those small numbers of planes. The barycentric line length LDA is defined as

LDA )

∑j |Rh j+1 - Rh j|

(18)

The straightness of this tunneling pathway can be measured by the index η as follows

η)

LDA RDA

(19)

where RDA is the direct distance between the donor and the acceptor. The straightness is higher when η is closer to 1. 3. Applications to the Ru-Modified Azurins We apply our method to some ruthenium-modified blue copper protein azurins. Experimentally, six surface positions of azurin were substituted individually by histidine (His) to which the Ru complex is coordinated.6,7 This yields six electrontransfer reactions that sample six different parts of the azurin protein from the native Cu+ to the substituted light-generated Ru3+. The electronic calculation is made for the whole protein on the basis of extended Hu¨ckel theory as before.24 We assume the following theoretical model: The donor electron is located in the 4s orbital of the copper, and the acceptor hole is the 5s orbital of the ruthenium, and the initial and final energies for the tunneling electron are -10.7 eV. The protein structure is taken from the crystallographic data for azurin,26-28 and the structure of the Ru complex is taken from the crystallographic data for Ru(bpy)2(im)2(BF4)2.29 The combined structure of the protein and Ru complex is obtained by minimizing the energy using the program PRESTO30 and the pertinent parameters.31-33 The extended Hu¨ckel calculation is made by using the QCPE34 and ITPACK2C35 programs. First, we calculate the interatomic tunneling currents. In Figures 2 and 3, we show the maps of the interatomic currents for |Kab| > 0.2, as well as the atoms that bear interatomic currents for |Kab| > 0.1 for the His83 and His124 derivatives, respectively. (Here, we use the terminology “His83 derivative” for the Ru-modified azurin in which the Ru complex is coordinated to His83). Clear illustrations of the interatomic currents for these derivatives have been presented before.24 The breadth of the arrow is chosen so as to be proportional to the magnitude of the current. In Figures 2 and 3, we see that many arrows are connected in a complicated fashion, forming a broad area for the tunneling route. By connecting the arrows, we obtained many chains of sequential atoms starting from donor and arriving at acceptor.24 Each chain of atoms might correspond to the chain of atoms in the pathway model proposed by Beratan and Onuchic.1,2,5 However, the chain of atoms formed by connecting the interatomic currents is more delicate and informative than the previous pathway model because the current is essentially a vector that retains a property of direction.

Figure 4. Plot of Qj as a function of the distance R of the plane Sj from the donor for His83 (solid line) and His124 (broken line). The interval for the each plane is 0.001 Å.

The vector property of the interatomic currents is most evidently expressed in the quantity Qj defined by

Qj )

∑ ∑

|Kab|

(20)

i b∉Ωi a∈ΩD D

It was already shown in eq 14 that the sum of Kab across the plane Sj is always unity. Therefore, the sum of |Kab| in eq 20, namely Qj, must be larger than unity. Large values of Qj indicate that significant backward currents (namely, negative value of Kab) flow through the plane Sj, as well as forward currents. We plot Qj for all of the planes at intervals of 0.001 Å for the His83 (solid line) and His124 (broken line) derivatives in Figure 4. We find that Qj for His83 is as much as 8 around the central region between the donor and acceptor. We also find that Qj for His124 is about 3 in the central region and is much larger in the region near the donor. These results are consistent with the maps of the interatomic currents in Figures 2 and 3, which show that complicated circular currents flow around the central region and the donor site region of the His83 derivative. Similarly, complicated circular currents flow around the donor region, and almost unidirectional currents flow around the central region in the His 124 derivative. Here, it should be mentioned that the circular current flowing through the plane Sj is as much as 3.5 times larger than the net current in the case of Qj ) 8. Judging from the above results, we can say that the description of the electron tunneling route by interatomic currents actually belongs to the same category as the pathway model and that the description by interatomic current involves proper information of the phase factor for each quantum mechanical electron pathway. Examining Figure 2 in more detail, we find that the overall tunneling current sometimes jumps from one β strand to another β strand. This arises from the structural property of the His83 derivative in which some parts of β strands run perpendicularly to the line connecting donor and acceptor. The overall pathways that we obtained with our method do not deviate much from those predicted by the tube model.3,4 However, our overall pathways are more compact than those of the tube model. This is shown in Figure 3, which suggests that one β strand plays the role of the dominant pathway while another β strand, which runs parallel to the first, plays the role of a subsidiary pathway in the central region of the overall process for the His124 derivative. In both derivatives, the pathways in the region around



Figure 5. Pictures of the worms in the His83 and His124 derivatives drawn on the background of the whole protein structure25 of azurin.

the donor and acceptor are very complicated, allowing for many through-space tunnelings. As we have seen above, the overall shape of the tunneling route of the His83 derivative is considerably different from that of the His124 derivative. The overall shapes of the tunneling routes in other four derivatives (His107, His109, His122, and His126) are similar to that of His124. In particular it is evident that only one β strand is used for electron tunneling in His107, forming an unusually narrow tunneling route. To more easily grasp a conceptual image for the electron tunneling route, we go from the description of a bunch of pathways to the description of a statistically averaged path

yielding the worm model. We calculate the line of the barycenter of the currents and the standard deviation, which becomes the radius of the worm on each cutting plane by the method in section 2. The worm so obtained is drawn by the shaded area in Figures 2 and 3. We also draw the worm on the background of the protein structure in Figure 5 for the His83 and His124 derivatives and in Figure 6 for the His107, His109, His122 and His126 derivatives. We see that these worms are rather narrow and winding. We average the standard deviation σj over all of the planes. In Table 1, we list the thus-calculated average value σ. We find that σ is rather small, ranging from 1.16 Å (His107) to 1.37 Å


Kawatsu et al.

Figure 6. Pictures of the worms in the His122, His109, His126, and His107 derivatives drawn on the background of the whole protein structure25 of azurin.

TABLE 1: Calculated Values of the Average Standard Deviation σ from the Barycenter, Barycentric Line Length LDA, and Straightness η number

derivative

σ (Å)

RDAa (Å)

LDA (Å)

η

1 2 3 4 5 6

His83 His107 His109 His122 His124 His126

1.37 1.16 1.35 1.22 1.28 1.32

17.1 26.0 16.7 12.7 19.7 21.6

20.3 29.8 20.2 16.9 24.7 25.4

1.18 1.14 1.20 1.32 1.25 1.17

a

RDA is the distance between donor and acceptor

(His83). This σ value is read as an average radius of the worm. In the same table, we also list the barycentric line length LDA and the direct distance between the donor and acceptor RDA, as well as the straightness η of the worm. We find that LDA ranges from 16.9 Å (His122) to 29.8 Å (His107) and that η ranges from 1.14 (His107) to 1.32 (His122). The finding that η is not much larger than unity, indicating that the worm deviates only slightly from linear. We also find a correlation that η is small when RDA is large, except in the case of His124. Next, we investigate the distribution of two kinds of currents Kab and |Kab|. For this purpose, we define two kinds of the



Figure 7. Plot of X(pσ) and Y(pσ) as functions of pσ for the His83 derivative. The solid line expresses Y(pσ) ) ∑|Kab|. The broken line expresses X(pσ) ) ∑Kab. σ is the standard deviation for |Kab|.

average current within a worm whose radius has been modified to pσ as follows

X(pσ) )

1

N

∑( ∑ N j)1 ab

Kab)

(21)

|Kab|)

(22)

(|Rj-Riab| σ. Even under such a situation, we use σ as the radius of the worm for convenience of description. We should again state the fact that some amount of the current |Kab| flows through the region outside of it, the manner of which varies among the derivatives. 4. Dependence of ke on the Worm Properties We go on to investigate how the electron-transfer rate ke depends on the distance RDA and the other various properties of the worm. We can calculate the electron-transfer rate kcal e using eq 1. For this purpose, we assume that (FC) is 10-4 cm in all cases. We evaluate TDA by the above extended Hu¨ckel calculations. Calculated values of kcal e roughly agree with the experimental values.17,18 In the following theoretical analysis, however, we use kcal e instead of the experimental value. This is a theoretically consistent way to investigate the correlation between the structures of the Ru-modified azurins determined by calculations and kcal obtained by using the calculated e structures. As we have already shown, we obtained worms for the six derivatives and the corresponding values of kcal e . We now try to establish the best-fit function for the six values of kcal e using the properties of the worms. Typical parameters of the worm are RDA and LDA. In addition to these, the atom density in the worm region can be another parameter. Specifically, we count the number of atoms nw whose centers are on the inside of the worm. We list those numbers for the six derivatives for the cases including and excluding hydrogen atoms in Table 2. Using these nw values, we can consider two kinds of one-dimensional atom


Kawatsu et al.

TABLE 3: Portion P of the Net Tunneling Current Passing through the Worm for the Six Derivatives number

derivative

P

1 2 3 4 5 6


0.81 0.71 0.75 0.81 0.81 0.16

densities in the worm, FL and FR, as follows:

FL )

nw LDA

(23)

FR )

nw RDA

(24)

It should be noted that we did not divide nw by the volume of the worm in eqs 23 and 24. We obtain atom densities along the one-dimensional coordinates LDA and RDA. Then, we obtain larger atom densities for larger worm cross section. Such atom densities must be pertinent in general conduction phenomena including electron tunneling. The calculated values of the atom densities for the six derivatives are listed in Table 2. Another property of the worm is the portion P of the net tunneling current through the worm, which is obtained from

∫0σX(r) dr P) ∞ ∫0 X(r) dr

(25)

where X(r) is the function defined in eq 22. The calculated values of P for the six derivatives are listed in Table 3. We find that P is almost always around 0.8, except for the case of the His126 derivative (P ) 0.16). We recognize the fact that P2kcal represents the rate due to the tunneling net current e through the worm. Now, let us consider the following five fitting functions

) a1 - b1LDA

(26)

log10 kRe ) a2 - b2RDA

(27)

log10 kLe

) a3 - (c3 - d3FL)(LDA - e3)

(28)

log10 kRF e ) a4 - (c4 - d4FR)(RDA - e4)

(29)

log10 kLF e

log10 kRFP ) a5 - (c5 - d5FR)(RDA - e5) - 2 log10 P e

(30)

where ai, bi, ci, di, and ei are parameters to be optimized. These parameters are determined by minimizing the following values 6

DI2 )

2 I [log10 kcal ∑ e (i) - log10 ke(i)] i)1

(31)

where i is the number of the derivative, as shown in Table 2, and the index I represents L, R, LF, RF, and RFP. DI2 is the variance, and DI is the standard deviation of the six points of log10 kcal e from the fitting function.

RF RFP Figure 9. Correlations of kLe , kRe , kLF with kcal e , ke , and ke e . The diagonal line indicates that the correlation is perfect.

TABLE 4: Standard Deviation of the Six Values of log10 kcal e from the Fitting Function log10 kIe I

DI (log10 unit)

L R LF RF RFP

0.860 0.770 0.514 0.360 0.266

Performing the minimization of eq 31, we obtain the leastsquares fitting functions as follows:

log10 kLe ) 13.0 - 0.443(LDA - 6.10)

(32)

log10 kRe ) 13.0 - 0.458(RDA - 2.78)

(33)

log10 kLF e ) 10.3 - (1.38-1.14FL)(LDA - 14.1)

(34)

log10 kRF e ) 10.5 - (1.19 - 0.780FR)(RDA - 9.19)

(35)

log10 kRFP ) 8.28 - (2.19 - 1.90FR)(RDA - 12.7) e 2 log10 P (36) for the six In Figure 9, we plot log10 kIe vs log kcal e derivatives. In Table 4, we list the standard deviation DI. We RF LF R see that the fitting improves in the order kRFP e , ke , ke , ke , and L ke . Above all, it is a remarkable result that the direct distance RDA is a better measure than the worm length LDA. The density also works rather well. (In particular, all of the filled circles and open circles of His83, His109, and His126, which deviated largely from the straight line in Figure 9, decreased the deviation considerably by considering the role of the atom density.) Finally, it should also be noted that the standard deviation DI becomes smaller when the portion of the net tunneling current P is introduced. Now, we notice that the fitting function log10 kRF e in eq 35 looks like the empirical equation that was recently proposed by the Dutton group in the investigation of the natural engineering principles of electron tunneling in biological oxidation-reduction processes11

log10 kem e ) 13.0 - (1.2 - 0.8F)(RDA - 3.6)

(37)



where we neglected the term due to the (FC) factor. In eq 37, F is the packing density. They obtain F by calculating the packing fraction of the atoms inside the pipe, which is determined by connecting all of the atoms of the donor with all of the atoms of the acceptor by straight lines. To make a comparison easier, we rewrite eq 35 as follows

log10 kRF e ) 17.2 - 4.36FR - (1.19 - 0.780FR)(RDA - 3.60) (38) ≈ 13.4 - (1.19 - 0.780FR)(RDA - 3.60)

(39)

where the near equality in eq 39 indicates that we made the approximation of replacing 4.36FR in eq 39 with the average value 4.36 × 0.86. Equation 39 agrees very well with eq 37. The physical meaning of F is not necessarily the same as our atom density FR, but its value is similar. In particular, F takes a value between 0.5 and 0.9 with the most frequent value being 0.75, while FR takes a value between 0.7 and 1.0 with an average of 0.86. Therefore, we provided a theoretical foundation to the empirical eq 37 of Dutton group by our microscopic calculations based on the worm model. 5. Discussion Before going into a detailed discussion, we summarize our strategy and the general concepts applied to our study. To begin, we theoretically calculated all of the interatomic currents for the six derivatives. Within these calculations, we obtained all of the information regarding all of the microscopic tunneling processes. The problem was then reducing the vast amount of microscopic information into an intuitively understandable electron tunneling route. For this purpose, we selected and included normalized interatomic currents whose amplitudes are larger than 0.1. The number of those currents is not prohibitively large so as to prevent their being written down for the protein structure on the atomic level. Examples are given in Figures 2 and 3. This is an intermediate stage of abstracting the microscopic information. We emphasized that the dominant electron tunneling pathways are formed by connecting the current vectors with large amplitudes. At this level, a group of the tunneling pathways corresponds to the pathway model proposed by Beratan and Onuchic, and Gray’s group.3,4 Also, we pointed out that our pathways retain the information of the electron tunneling direction pertinent to the current and that many of the currents are circulating. The result in Figure 4 that the value of Qj can be larger than 8 in some regions indicates the large scale of current circulation. In the next step of abstracting the information, we calculated the barycenter and the standard deviation of all the currents on the plane that was placed perpendicularly to the straight line connecting the donor and acceptor. By connecting the barycenters, we obtained the central line of the tunneling route. By connecting the circles with radii equal to the standard deviations on many planes, we obtained the worm. The worm was considerably winding (the straightness of the central line of worm was in the range 1.14-1.32 for the six derivatives). The average radius of the worm was 1.16-1.37 Å for the six derivatives, which is similar to the sizes of the radii of the donor and acceptor ions. Using this worm model, we obtained some fitting functions for the six values of log10 kcal e . One of the , agreed rather well with the empirifitting functions, log10 kRF e cal function proposed by Dutton group. Although there is a qualitative difference between the atom density and the packing density, we could successfully reproduce the empirical relation

Figure 10. Schematic representation of the worm as a stack of cylinders.

using the somewhat advanced level of the worm model, which is the model abstracted from the pathway model. In this sense, we can say that the pathway model and the Dutton plot (empirical relationship of ke with RDA and F) are consolidated theoretically. However, it should be noted that the region in which F is calculated is different from the region in which FR is calculated. In this respect, it can be said that the Dutton plot does not reflect the real electron tunneling pathway but rather superficially reproduces a good correlation between ke and RDA with the mediation of F. In the previous section, we found that RDA is a better factor for fitting kcal F than LDA. It is an intriguing question why the worm length LDA is not a better factor than RDA. As mentioned in section 2, we positioned parallel planes at intervals of 0.001 Å vertically with respect to the straight line connecting the donor and acceptor. We calculated the radius of the worm on each plane using eq 16. Then, we can define a cell that is a cylinder with the radius of the worm and a height of 0.001 Å. The worm is expressed as the stack of these cells, as drawn schematically in Figure 10. The tunneling current through the cell is conserved among the cells. Furthermore, the radius of the cylindrical cell does not change very much. Therefore, the “resistance” to electron tunneling is essentially the same among the cells. The total resisance to electron tunneling is expressed by an exponential function of the product of the number of cells and the average resistance of the cells. The number of cells is proportional to RDA, namely, RDA/0.001. The resistance of the cell should be a function of the number of atoms in the cell per unit length, called the one-dimensional atom density in the worm. The one-dimensional atom density in the cell is, on the average, equal to the number of atoms in the worm divided by RDA, as expressed by FR in eq 24. The decay function of the tunneling current is expressed by such a total resistance to electron tunneling. Then, the logarithm of the tunneling rate should be expressed by a form such as eq 29. Thus, eq 29 is one of the natural formulas in expressing the present analysis of electron tunneling current. On the other hand, if we try to express the total resistance using LDA, we should initially put the planes vertical with respect to the center line of the worm, which must be obtained after determining barrycenters on the planes. That is, the planes must be determined selfconsistently, which is very complicated and laborious. After that, the radius of the worm, the dispersion on the plane, is determined. The worm produced by such process might be somewhat different from the worm that we obtained in the present paper. In eq 23, we assumed that the worm is the same between the two


Kawatsu et al.

TABLE 5: Number of Heavy Atoms (C, N, O, S) and Light Atoms (H) within the Worm heavy atoms number i 1 2 3 4 5 6 total

light atoms

derivative

nwa

ncb

nc/nw

nw

ncc

nc/nw


8 22 13 7 16 12

8 22 11 7 16 12

1.00 1.00 0.85 1.00 1.00 1.00

5 2 4 2 3 5

3 (0) 2 (1) 1 (0) 1 (0) 2 (1) 4 (1)

0.60 1.00 0.25 0.50 0.67 0.80

78

76

0.97

21

13

0.62

a

nw is the number of atoms that belong to the worm. b nc is the number of atoms that contribute to the current |Kab| > 0.1 within the worm. c The number in parentheses is the number of hydrogen atoms forming hydrogen bond.

processes in calculating FL. Therefore, some error must be involved in such an analysis using LDA, namely, in eq 28. This might be the reason that the analysis using RDA is better than that using LDA. This fact also indicates that the detail of the form of the worm changes somewhat depending on the analytical method. However, we think that the present method is simple and that an essential part of the electron tunneling pathway is sufficiently involved. Next, we discuss the role of hydrogen atoms in the electron tunneling route. In Table 5, we list the number nw of heavy atoms (C, N, O, S) and light atoms (H) within the worm for the six derivatives. We also list the number nc of heavy atoms and light atoms that contribute to current |Kab| that are larger than 0.1 within the worm. For the heavy atoms, nc/nw is 1.0, except for the His109 derivative, indicating that almost all of the heavy atoms in the worm work to form an effective route for the tunneling current. For the light atoms, nc/nw ranges from 0.25 to 1.00, with an average value of 0.62. This indicates that the hydrogen atoms also work actively to form an effective route in the worm. The hydrogen atoms that work as hydrogen bonds represent relatively few cases, as shown by the numbers in parentheses in Table 5. These results indicate that the role of hydrogen atoms should be taken into account explicitly in the calculations, as well as the role of heavy atoms. Finally, let us reexamine the limitations of the present calculations. First, the basis of our theoretical calculations was the extended Hu¨ckel theory. This theory is easy to use, but it contains the disadvantage that the electron correlation effect is not correctly taken into account. Therefore, it is a general consensus that the results calculated using the extended Hu¨ckel theory are considered to be qualitative at best. In this sense, most of our calculated values are only qualitatively correct. Even if this is the case, how could we derive so many conclusions? The answer is that we consistently made a mostly comparative study among the six derivatives within the theoretical framework, to establish the method of defining tunneling region. For this purpose, we calculated the interatomic currents and proposed the worm as an abstraction. We defined the worm radius σ by the standard deviation on each plane. As we have seen in Figures 7 and 8, this definition of σ still reserves considerable ambiguity. This ambiguity comes not only from the qualitative nature of the extended Hu¨ckel calculations but also from the ambiguous concept of the tunneling route up to the present time. Rather than treating various kinds of redox proteins, we studied rather similar kinds of redox proteins, namely, the six Ru derivatives of azurins. Then, we succeeded in consolidating the pathway model and Dutton plot, both of which could be reproduced by our calculations at different levels of abstraction and approximations. Taking advantage of the experience in analyzing the

tunneling route of the six derivatives, we go on to different types of redox proteins. This becomes our future area of study. In the present study, we used redox systems with the same donor and the same kind of acceptors. This is the reason that we obtained a very simple relation between ke and RDA using worm model. When we applied the present theory to other systems such as cytochrome c oxidase, we could not obtain such a nice relationship as obtained in the present calculation.36 This would indicate that the property of the local structure of the redox site must participate significantly in the electron tunneling pathway and modify the relationship between ke and RDAthrough the worm model. This problem is for our future investigations. 6. Conclusion In this paper, we showed that the average interatomic tunneling currents reduce to a worm that winds with a width of a few angstroms. The manner of the winding of the worm reflects the places through which electron tunneling is either easy or difficult, depending on the microscopic environment of the protein. In particular, the protein environment is not homogeneous for the tunneling current. Even under such a situation, we showed that the electron transfer rate is rather well expressed by an exponentially decreasing function of the donoracceptor distance RDA, with the decay constant being somewhat dependent on the atom density in the worm. This result is consistent with the recently obtained Dutton plot. We elucidated the mechanism by which RDA becomes a better factor in the correlation with the electron-transfer rate than the worm length LDA. In such a way, we succeeded in consolidating the electron tunneling pathways model with the Dutton plot by use of the worm model. Acknowledgment. This work was supported by the Grantin-Aid of Scientific Research (B) to T.K. and Scientific Research on Priority Areas “Biological Machinery” to T.Y. from the Japanese Ministry of Education, Science, Sports and Culture. References and Notes (1) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488. (2) Onuchic, J. N.; Beratan, D. N.; Winkler, J. R.; Gray, H. B. Annu. ReV. Biophys. Biomol. Struct. 1992, 21, 349. (3) Regan, J. J.; Bilio, A. J. D.; Langen, R.; Skov, L. K.; Winkler, J. R.; Gray, H. B.; Onuchic, J. N. Chem. Biol. 1995, 2, 489. (4) Regan, J. J.; Onuchic, J. N. Electron Transfer: From Isolated Molecules to Biomolecules, Part Two; Jortner, J., Bixon, M., Eds.; John Wiley & Sons: New York, 1999; p 467. (5) de Andrade, P. C. P.; Onuchic, J. N. J. Chem. Phys. 1998, 108, 4292. (6) Wuttke, D. S.; Bjerrum, M. J.; Winkler, J. R.; Gray, H. B. Science 1992, 256, 1007. (7) Casimiro, D. R.; Richards, J. H.; Winkler, J. R.; Gray, H. B. J. Phys. Chem. 1993, 97, 13073. (8) Langen, R.; Chang, I.-J.; Germanas, J. P.; Richard, J. H.; Winkler, J. R.; Gray, H. B. Science 1995, 268, 1733. (9) Gray, H. B.; Winkler, J. R. Annu. ReV. Biochem. 1996, 65, 537. (10) Moser, C. C.; Keske, J. M.; Warnucke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 335, 796. (11) Page, C. C.; Moser, C. C.; Chen, X.; Dutton, P. L. Nature 1999, 402, 47. (12) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1990, 94, 8430. (13) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 2440. (14) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 6111. (15) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1993, 97, 13078. (16) Stuchebrukhov, A. A. Chem. Phys. Lett. 1994, 225, 55. (17) Stuchebrukhov, A. A.; Marcus, R. A. J. Phys. Chem. 1995, 99, 7581. (18) Gehlen, J. N.; Daizadeh, I.; Stuchebrukhov, A. A.; Marcus, R. A. Inorg. Chim. Acta 1996, 243, 271. (19) Okada, A.; Kakitani, T.; Inoue, J. J. Phys. Chem. 1995, 99, 2946. (20) Stuchebrukhov, A. A. J. Chem. Phys. 1996, 105, 10819.

Worm Model for Electron Tunneling in Proteins (21) Stuchebrukhov, A. A. J. Chem. Phys. 1997, 107, 6495. (22) Stuchebrukhov, A. A. J. Chem. Phys. 1996, 104, 8424. (23) Cheung, M. S.; Daizadeh, I.; Stuchebrukhov, A. A.; Heelis, P. F. Biophys. J. 1999, 76, 1241. (24) Kawatsu, T.; Kakitani, T.; Yamato, T. Inorg. Chim. Acta 2000, 300-302, 862. (25) Humphrey, W. F.; Dalke, A.; Schulten, K. J. J. Mol. Graphics 1996, 14, 33. (26) Bonander, N.; Va¨nngård, T.; Tsai, L.-C.; Langer, V.; Nar, H.; Sjo¨lin, L. Proteins 1997, 27, 385. (27) Adman, E. T.; Jensen, L. H. Isr. J. Chem. 1981, 21, 8. (28) Adman, E. T.; Stenkamp, R. E.; Sieker, L. C.; Jensen, L. H. J. Mol. Biol. 1978, 124, 35. (29) Reddy, K. B.; Cho, M.-o. P.; Wishart, J. F.; Emge, T. J.; Isied, S. S. Inorg. Chem. 1996, 35, 7241.

J. Phys. Chem. B, Vol. 105, No. 19, 2001 4435 (30) Morikami, K.; Nakai, T.; Kidera, A.; Saito, M.; Nakamura, H. PRESTO (A vectorized molecular mechanics program for biopolymers). Comput. Chem. 1992, 16, 243. (31) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (32) Geremia, S.; Calligaris, M. J. Chem. Soc., Dalton Trans. 1997, 1541. (33) Ungar, L. W.; Scherer, N. F.; Voth, G. A. Biophys. J. 1997, 72, 5. (34) Howell, J.; Rossi, A.; Wallace, D.; Haraki, K.; Hoffmann, R. FORTICON8 (Extend Hu¨kel method program). QCPE 1977, 11, 344. (35) Kincaid, D. R.; Respess, J. R.; Young, D. M.; Grimes, R. G. ITPACK, version 2C; a Fortran package for solving large sparse linear system by adaptive accelerated iterative methods; Center for Numerical Analysis: University of Texas, Austin, Texas, 1999. (36) Kawatsu, T.; Kakitani, T.; Yamato, T. unpublished results.

Worm Model for Electron Tunneling in Proteins: Consolidation of the

Recommend Documents