Average Electron Tunneling Route of the Electron Transfer in Protein

9948

J. Phys. Chem. B 2008, 112, 9948–9958

Average Electron Tunneling Route of the Electron Transfer in Protein Media Hirotaka Nishioka† and Toshiaki Kakitani*,‡ Graduate School of EnVironmental and Human Sciences, and Department of General Education, Faculty of Science and Technology, Meijo UniVersity, Tempaku-ku, Nagoya 468-8502, Japan ReceiVed: NoVember 08, 2007; ReVised Manuscript ReceiVed: April 09, 2008

We present a new theoretical method to determine and visualize the average tunneling route of the electron transfer (ET) in protein media. In this, we properly took into account the fluctuation of the tunneling currents and the quantum-interference effect. The route was correlated with the electronic factor 〈TDA2〉 in the case of ET by the elastic tunneling mechanism. We expanded 〈TDA2〉 by the interatomic tunneling currents Jab’s. Incorporating the quantum-interference effect into the mean-square interatomic tunneling currents, denoted as 〈Jãb2〉, we could express 〈TDA2〉 as a sum of p2〈Jãb2〉. Drawing the distribution of 〈Jãb2〉 on the protein structure, we obtain the 〈Jãb2〉 map which visually represents which parts of bonds and spaces most significantly contribute to 〈TDA2〉. We applied this method to the ET from the bacteriopheophytin anion to the primary quinone in the bacterial photosynthetic reaction center of Rhodobacter sphaeroides. We obtained Jab’s by a combined method of molecular dynamics simulations and quantum chemical calculations. In calculating 〈Jãb2〉, we found that much destructive interference works among the interatomic tunneling currents even after taking the average. We drew the 〈Jãb2〉 map by a pipe model where atoms a and b are connected by a pipe with width proportional to the magnitude of 〈Jãb2〉. We found that two groups of 〈Jãb2〉’s, which are mutually coupled with high correlation in each group, have broad pipes and form the average tunneling routes, called Trp route and Met route. Each of the two average tunneling routes is composed of a few major pathways in the Pathways model which are fused at considerable part to each other. We also analyzed the average tunneling route for the ET by the inelastic tunneling mechanism. 1. Introduction In biological systems, long-range electron transfer (ET) reactions between protein cofactors take place via a superexchange mechanism where the electron tunnels with use of virtual states provided by protein environments.1,2 The protein environment has the heterogeneous and intrinsic structure. It has been a significant subject to elucidate how these protein structures contribute to regulate the long-range electron tunneling route and the ET rate.3–12 A standard theory for the ET rate kDA gives the following formula3,13

kDA(-∆G) )

2π |T |2(FC(-∆G)) p DA

(1)

where the term (FC(-∆G)) is called the thermally averaged Franck-Condon factor which depends on the energy gap -∆G and |TDA|2 is called the electronic factor for the electron tunneling. In the classical approximation of the nuclear motion, the (FC(-∆G)) was expressed by the Marcus energy gap law.3 In many biological ET systems, it is known that the protein environment is adjusted to optimize the energy gap law.2 For the electronic factor, there has been much debate as to how the electron tunneling route and TDA are regulated by the protein environment.2,4,6,7,11,14–23 Many theoretical approaches to determine the electron tunneling route for a fixed protein conformation have been done for these 20 years.4,6,14 Those studies involve the Pathways * Corresponding author. Phone: +81-52-838-2394. Fax: +81-52-8321170. E-mail: [email protected]. † Graduate School of Environmental and Human Sciences. ‡ Department of General Education, Faculty of Science and Technology.

model,15–17 tube model,18,19 tunneling current analysis,7,11,20 worm model,21 and the electronic contact map22 as well as a simpler model of the square-barrier tunneling2 and the packing density model.23,24 On the other hand, theoretical studies by the combined analysis of molecular dynamics (MD) simulations for the protein fluctuation and the quantum chemical (QC) calculations for the TDA revealed that the value of TDA fluctuates significantly according to the thermal fluctuation of protein conformation.25–31 In more detail, it was found that the interatomic tunneling currents change drastically according to the fluctuation of protein conformation.20,21,28,29 Then, the tunneling routes which are obtained by the connection of the interatomic tunneling currents are altered significantly by the fluctuation of the protein conformation. Since the tunneling currents are vectorial quantities, the summation of the tunneling currents which give rise to TDA7,21 is suffered from the quantum-interference among the interatomic tunneling currents. Here, let us define a route as follows: The route is a group of the interatomic tunneling currents which have large amplitude and the fluctuations of which have high correlation among them. This group should satisfy the condition that the combination of the interatomic tunneling currents in the group spans from donor to acceptor. Under such situation, it is desired to take proper average for the fluctuation of the interatomic tunneling currents by taking into account the varying quantum-interference effect and to obtain reasonably averaged tunneling routes. Recently, much theoretical interest was directed to the nonCondon effect for the ET rate under the fluctuating protein conformation.25,32–37 Among them, a non-Condon theory of the ET rate which satisfies the detailed balance condition was presented.37 According to this theory, the ET rate is expressed

10.1021/jp710689s CCC: $40.75  2008 American Chemical Society Published on Web 07/16/2008

Average Electron Tunneling Route

J. Phys. Chem. B, Vol. 112, No. 32, 2008 9949

as the sum of two terms; kDAel due to the elastic tunneling mechanism and kDAinel due to the inelastic tunneling mechanism as follows:37 el inel kDA(-∆G) ) kDA (-∆G) + kDA (-∆G) el (-∆G) ) kDA inel (-∆G) ) σDA2 kDA

1 p2

2π 〈T 2 〉 (FC(-∆G)) p DA

∫-∞∞ dt(A(t) - 1) exp(iεt/p)

∑ Ciνφν

(5)

∑ Cfνφν

(6)

ψi ) CDiφD +

(2)

ν

(3)

ψf ) CAfφA +

ν

∫-∞∞ dε(FC(-∆G - ε))

2 × 1 + exp(-ε/kBT)

Adopting the one electron and effective reduced two states approximation, we express the molecular orbitals ψi and ψf in the initial and final diabatic states as

(4)

where 〈 〉 indicates the statistical average. The factor A(t) is the normalized autocorrelation function of TDA(t). The factor σDA2 is the variance of TDA; σDA2 ) 〈TDA2〉 - 〈TDA〉2. Equation 3 was also derived in the other non-Condon theories under the assumption of slow modulation of the molecular environment.34–36 It was shown that the inelastic tunneling mechanism works dominantly for the energy gap deeply in the inverted region in the case of fast modulation of the molecular environment.25,32,37 Equations 2 and 3 represent that 〈TDA2〉 is the electronic factor of kDAel and that σDA2 is the formal electronic factor of kDAinel. The integral of eq 4 denotes how the thermally averaged Franck-Condon factor (FC(-∆G - ε)) is modulated by the non-Condon effect of the electronic factor. In general, conformational fluctuations of proteins are composed of fluctuations in broad time scales ranging from femtoseconds to seconds. It is impossible to cover such broad time range by the MD simulation. The time range which is covered by the ordinary MD simulations at present will be some nanoseconds. Therefore, we cannot expect that the time average over around 1 ns will correspond to the thermal average in the equilibrium state. This short time average will be related to the average for the fluctuation of the protein environment centered at a global potential energy minimum so far as we start from the protein conformation obtained by the X-ray crystallographic study. In this paper, we present a method in which the electronic factor is decomposed into the sum of the mean-square interatomic tunneling currents incorporating the quantum-interference effect among them. This method is applied to the ET from the bacteriopheophytin anion (BPhe-) to the primary quinone (QA) in the bacterial photosynthetic reaction center (BPRC). We derived average tunneling routes by analyzing the magnitude of the mean-square interatomic tunneling currents and considering the quantum-interference effect among the interatomic tunneling currents. The constitution of this paper is as follows: In section 2, we present a theory for expressing the electronic factor of the ET rate (eq 1) with use of the fluctuating interatomic tunneling currents. The thermal fluctuation of the protein conformation is produced by the MD simulation. Then, in section 3, we describe the method of MD simulations for BPRC and the method of calculation of the electronic structure. In section 4, we present numerically calculated results of the average electron tunneling routes for the elastic and inelastic tunneling mechanisms. Section 5 is devoted to the discussion. In section 6, we make our conclusion. 2. Theory 2.1. Interatomic Tunneling Currents. First, we briefly review the theoretical formula of the interatomic tunneling current and the electron tunneling matrix element.7,11,20,28

where φD and φA are the donor and acceptor molecular orbitals, respectively. The function φν represents the atomic orbital (AO) provided by the mediator. The interatomic tunneling current Jab is written as7,11,20,28

Jab )

∑

Jµν

(7)

µ∈a,ν∈b

1 Jµν ) (Cµi Cfν - Cµf Ciν)(Hµν - EtSµν) p

(8)

where H is the one-electron Hamiltonian, S is the overlap integral, and Et is the tunneling energy. The tunneling matrix element TDA is expressed as7,20

TDA ) p

∑

Jab

(9)

j ,b∉Ωj a∈ΩD D

where ΩDj is a donor side space separated from acceptor side by an arbitrary surface Sj. The renormalization of ψi and ψf is made and its effect is taken into account in calculating the tunneling matrix element and the interatomic tunneling currents by the method described in refs 39 and 40. The sign of Jab is taken positively when the tunneling current flows from donor to acceptor. A consistent set of signs is obtained for all of the Jab’s which contribute to eq 9, and this then leads to an overall sign for TDA. 2.2. Expression of the Mean-Square Electron Tunneling Matrix Element 〈TDA2〉 by Means of the InterferenceCorrected Mean-Square Interatomic Tunneling Currents. The ET rate due to the elastic tunneling mechanism is expressed as eq 3. The electronic factor is 〈TDA2〉. In the following, we decompose 〈TDA2〉 into the mean-square interatomic tunneling currents corrected by the quantum-interference effect. Using eq 9, we rewrite 〈TDA2〉 as sum of two terms: the diagonal term 〈Dj〉 and the nondiagonal term 〈NDj〉 of the product of the two interatomic tunneling currents as follows:

〈

〈TDA2 〉 ) p2

∑

Jab


∑

Jcd

j ,d∉Ωj c∈ΩD D

〉 ) 〈Dj 〉 + 〈NDj 〉 (10)

where

Dj ≡ p2

∑

Jab2

(11)


NDj ≡ p2

∑

∑

JabJcd

(12)

j ,b∉Ωj c(ab*cd)∈Ωj ,d∉Ωj a∈ΩD D D D

Here the term Dj indicates the sum of square values of Jab passing through an arbitrary surface Sj that surrounds the donor region ΩDj. The term NDj indicates the sum of JabJcd(ab*cd), in which both Jab and Jcd pass through the surface Sj. The nondiagonal terms represent the magnitude of the quantuminterference among tunneling currents. The time averages of Dj and of NDj on each surface Sj are calculated as follows:

9950 J. Phys. Chem. B, Vol. 112, No. 32, 2008

1 N

〈Dj 〉 )

1 N

〈NDj 〉 )

Nishioka and Kakitani

N

∑

Dj(tk)

(13)

k

N

∑ NDj(tk)

(14)

k

where tk represents MD time step. Here, we define a parameter δj representing the degree of quantum-interference as follows

δj )

〈NDj 〉 〈Dj 〉

(15)

Then, 〈TDA2〉 can be written as follows

〈TDA2 〉 ) 〈Dj 〉 (1 + δj)

(16)

where -1 < δj < 1 holds true. When δj is close to 1, the quantum-interference is totally constructive among the interatomic tunneling currents passing the surface Sj. When δj is close to -1, the quantum-interference is totally destructive. When δj is zero, the quantum-interference is totally absent. Here, we used the words “quantum-interference is constructive or destructive” to discriminate it from the words “constructive or destructive interference” used in the coherence parameter.27 Equation 16 represents that 〈TDA2〉 is expressed as 〈Dj〉 or the sum of the mean-square interatomic tunneling currents 〈Jab2〉 time p2 passing through an arbitrary surface Sj if we can neglect the quantum-interference effect among the interatomic tunneling currents. However, when the quantum-interference effect is not small, we rewrite eq 16 as

〈TDA2 〉 ) p2

∑

〈Jãb2 〉

(17)


where

〈Jãb2 〉 ) 〈Jab2 〉 (1 + δj) (a ∈ ΩDj, b ∉ ΩDj)

(18)

Equation 17 represents that 〈TDA2〉/p2 is strictly expressed as the sum of 〈Jãb2〉, which may be called the mean-square interatomic tunneling current corrected by the quantum interference. Similarly, the factor (1 + δj) in eq 18 may be called the quantum-interference correction factor for the mean-square interatomic tunneling currents. 2.3. Expression of the Variance σDA2 of the Electron Tunneling Matrix Element by Means of the Variance σ˜ ab2 of the Interatomic Tunneling Currents Corrected by the Quantum Interference. Using the interatomic tunneling currents, the variance of TDA is

σDA2 ≡ 〈TDA2 〉 -〈TDA〉2 ) p2

∑ 〈Jab2 〉 + p2∑ ∑ 〈JabJcd 〉 -p2(∑ 〈Jab〉)2 a,b

a,b c,d

a,b

)

∑ σãb2

(19)

a,b

where

ηj )

σ ãb2 ) σab2(1 + ηj)

(20)

σab2 ) p2(〈Jab2 〉 -〈Jab 〉2)

(21)

p2∑a,b∑c,d(〈JabJcd 〉 -〈Jab 〉 〈Jcd 〉) ∑a,bσab2

(22)

with a condition that a ∈ ΩDj, b ∉ ΩDj, c ∈ ΩDj, d ∉ ΩDj, and ab

* cd. Equation 19 represents that σDA2 is expressed as the sum of σ˜ ab2 which may be called the variance of the interatomic tunneling currents corrected by the quantum interference. Similarly, the factor (1 + ηj) in eq 20 may be called the quantum-interference correction factor for the variance of the interatomic tunneling currents. 3. Method 3.1. Molecular Dynamics. The procedures for molecular dynamics (MD) simulations reported in this paper were essentially same as those described in refs 29, 37, and 38. The initial configuration of the reaction center of Rhodobacter sphaeroides was obtained from the Protein Data Bank, entry code 1AIJ.41 The 262 water molecules incorporated in the structure 1AIJ are taken into account after optimization of their positions. The MD simulations are performed for the whole reaction center including 828 amino acids, all cofactors, 262 water molecules, and 4 detergent molecules of lauryldimethylamineoxide (LDAO) which are found to be attached to the RC in the structure 1AIJ. The phytil chain of the primary quinone is included up to C31 which was identified in the structure 1AIJ. The MD program Presto42 was used with the Amber force field.43 We adopted a dielectric constant of 2 and approximated Coulombic potentials over 12 Å by the PPPC method.44 The harmonic restriction was imposed on the heavy atoms of the residues on the surface of the protein and three LDAO molecules on the surface. The same restriction was imposed on oxygen atoms of water molecules. The Shake algorithm45 was used for bond stretchings of hydrogen atoms. One point which we newly adopt in the present study is that we restrict the position of the hydrophobic phytil chain in the region C21 to C31 since the carbon atoms of this region of the phytil chain are located outside the protein surface and no LDAO molecule is attached to them.41 (One of the three LDAO molecules on the surface covers the phytil chain up to C21 to protect it from exposure to the surface of the protein.) We first obtained the energy minimum conformation of protein. After such treatment, the system was gradually elevated to 300 K. An integration time step of 1 fs was employed. After 600 ps run for equilibration, we generated a trajectory for 1 ns. We collected the conformations of 1 000 001 at every 1 fs, which were used to calculate TDA and Jab. 3.2. Electronic Structure Calculations. The procedures for QC calculations reported in this paper were essentially the same as described in refs 37 and 38 except one point which we describe below. The chemical structures of donor and acceptor which we adopted in calculating φD and φA were shown in ref 29. In these calculations, the long hydrocarbon chains of native bacteriopheophytin a and the native ubiquinone 10 are truncated to hydrogen atoms. We calculated TDA and Jab for a pruned system which consists of BPhe, QA, and the three amino acids TrpM252, MetM218, and HisM219 at each MD time step. In this ET system, these three amino acids were already shown to be necessary and sufficient to calculate TDA in the electron transfer from BPhe- to QA.27,29 The coordination of donor, acceptor, and the three amino acids is shown in Figure 1. The electronic structures of the pruned protein were solved at the extended Hu¨ckel level. We referred to the Forticon8 program46 for the extended Hu¨ckel calculations. The parameters used in the extended Hu¨ckel calculations were taken from the values in ref 50. We adopted the tunneling energy Et, as -9.5 eV.27,29 We used a ITPACK 2C package47 for solving {Ciµ} and {Cfµ}. The calculation of TDA for each protein conformation was made by the same method as the previous one.20,21,28,29,40



Figure 2. Calculated time-dependence of TDA for Et ) -9.5 eV (red line) and for Et ) -10.0 eV (green line).

Figure 1. General view of the coordination of donor (bacteriopheophytin), acceptor (quinone), and the three amino acids (TrpM252, MetM218, and HisM219) which we investigate in this paper. The molecular structure of donor and acceptor is represented by a ball and stick model. The molecular structure of the three amino acids is represented by a stick model. Blue balls, red balls, deep blue balls, and white balls indicate carbon atoms, oxygen atoms, nitrogen atoms, and hydrogen atoms, respectively. The gold stick in the MetM218 represents the sulfur atom.

The electronic structures of donor and acceptor were calculated by the PM3 method48 in the Gaussian package.49 The lowest unoccupied molecular orbital of the neutral BPhe is chosen as φD, and the lowest unoccupied molecular orbital of the neutral QA is chosen as φA. To save the computational time, we fix the atomic orbital coefficients in φD and φA which are solved in the isolated state of BPhe and QA, respectively. In both the PM3 method and the extended Hu¨ckel theory, the Slater type atomic orbitals are used. Even though the overlap integral is neglected in the PM3 method in solving eigen values and eigen functions, this treatment does not indicate that the atomic base is chosen as the orthogonal one but that it is just an approximation and the effect of approximation is reduced by properly choosing the Fock matrix element. Under this situation, we made connection of the electronic structures of the amino acids with donor (acceptor) within the framework of the extended Hu¨ckel theory. 37 One point which we newly adopt in the present study is that we take into account the variable hyperconjugation effect of the methyl group with the π-conjugated part of quinone in the course of rotation of the methyl group. (The methyl group which is obtained by truncating the phytil chain is not concerned with this problem because no truncation is made in the MD simulations and the rotation does not occur.) Namely, the molecular orbital (MO) coefficients of hydrogen atoms in the methyl group change in the course of methyl rotation. As an approximate method to incorporate this variable hyperconjugation effect, we calculated the MO coefficient of each hydrogen atom as a function of rotational angle θ of the methyl group in the isolated state of quinone. We obtain the MO coefficient of one hydrogen atom is well fitted to the following function

0.07284 × cos(θ - 90o)

(23)

where the origin of θ is chosen at the position when one of the C-H bond is parallel to π-conjugated plane of quinone. The MO’s for the other hydrogen atoms are obtained by shifting θ by ( 120°. Such treatment is crucial because we found that the

methyl group plays a significant role in the route through which the electron tunnels dominantly. The variable hyperconjugation effect of the methyl group was also applied to the methyl groups in BPhe. 4. Results 4.1. Dependence of TDA on the Tunneling Energy. In the treatment of the one electron and the effective two states model, the tunneling energy Et must be empirically determined. By the calculation of the extended Hu¨ckel theory, we find that the energy gap between HOMO and LUMO in our three amino acids system is between -11.19 and -8.04 eV. Then, the plausible value of Et was estimated to be around -9.5 to -10 eV which is middle of the energy gap.27,29 In Figure 2, we plot an example of the time dependence of TDA calculated by using Et ) -9.5 and -10 eV. We observe that TDA for Et ) -10.0 eV is only slightly smaller than TDA for Et ) -9.5 eV. In this paper, we tentatively adopt Et ) -9.5 eV. We compare the results of Et ) -10.0 eV with those of Et ) -9.5 eV if necessary. 4.2. Map of the Mean-Square Interatomic Tunneling Currents. First, we calculated the mean-square interatomic tunneling currents 〈Jab2〉 using the tunneling energy Et ) -9.5 eV. In Figure 3, we mapped 〈Jab2〉 by the red pipe whose width is proportional to the magnitude of 〈Jab2〉. In forming the map, we used an instantaneous conformation of the protein, since the methyl group rotates many times during the simulation time and the protein conformation other than the methyl group is not largely altered during the simulation time. We find that there is a series of the broad pipes connected to each other from the carbonyl oxygen atom at C131 in the ring V of BPhe through the π-conjugated rings of Trp to the two π-conjugated atoms C1 and C6 in QA. We also find that there is another series of the broad pipes connected to each other from C12 in the ring III of BPhe, through a part of the Met residue, to the three hydrogen atoms of the methyl group attached to C5 of QA. The residual pipes are extremely narrow except the two pipes with appreciable width connecting the two series of broad pipes. We notice that the pipe which connects the C12 atom in the ring III of BPhe to the closest hydrogen atom of the Met residue is broad. 4.3. Mean-Square Interatomic Tunneling Currents 〈Jãb2〉 Corrected by the Quantum Interference. In order to consider the quantum-interference effect, we calculate the factor (1 + δj) in eqs 15 and 16. For this purpose, we provide the surfaces


Figure 3. Calculated map of the mean-square interatomic tunneling current 〈Jab2〉. All of the values of the mean-square interatomic tunneling currents 〈Jab2〉 are drawn by the red pipe whose width is proportional to the value of 〈Jab2〉. Those pipes which satisfy p2〈Jab2〉/〈TDA2〉 < 0.001 are not appreciably seen. Protein conformation together with that of donor and acceptor are drawn by the ball and stick model. The blue, white, red, and gold balls indicate the carbon, hydrogen, oxygen, and sulfur atoms, respectively. The orange and green balls indicate the atoms constituting the conjugated part and nonconjugated part of donor and acceptor, respectively. The atomic view is produced with use of VMD.56

Figure 4. Plot of the calculated values of the correction factor (1 + δj) in the surface Sj at the distance R from the center of quinone ring.

Sj’s of infinite planes with separation of 0.1 Å which are perpendicular to the line connecting the centers of donor and acceptor molecules. In the case when the infinite plane is cut by the donor molecule, the part of the cut plane is replaced with the surface produced by the outer spheres with van der Waals’ radius of atoms in the donor. The similar thing is made also when the infinite plane is cut by the acceptor molecule. By such procedure, the equivalence of TDA with the sum of pJab in eq 9 is guaranteed for all the surfaces Sj’s. We calculated (1 + δj) as a function of the distance R of the infinite surface. The results are shown in Figure 4. The origin of R is chosen at the center of the quinone ring. The center of the conjugation ring of BPhe is at R ) 13.1 Å. We observe that (1 + δj) is considerably smaller than 1 in all of the values of R: As R increases from 0 Å to 5 Å, (1 + δj) increases stepwise from 0.28 to 0.45. As R increases from 5 Å to 8.3 Å, (1 + δj)


Figure 5. Calculated map of the mean-square interatomic tunneling currents corrected by the quantum interference, 〈Jãb2〉. All of the values of 〈Jãb2〉 are plotted by a pipe model with red color. The others are the same as Figure 3.

decreases stepwise from 0.45 to 0.13. These results indicate that much destructive quantum interference works in the regions near the donor and acceptor molecules. The stepwise change of (1 + δj) happens when the dominant tunneling bond or space is changed to the others as R changes. We calculate 〈Jãb2〉, mean-square interatomic tunneling currents corrected by the quantum-interference, using eq 18 and Figure 4. In this calculation, the correction factor (1 + δj) for 〈Jab2〉 was averaged for the surfaces Sj’s through which the interatomic tunneling current Jab flows. Its result is mapped in Figure 5 by the pipe representation. The width of the pipe is proportional to the relative magnitude of the 〈Jãb2〉. However, the pipe width of Figure 5 is expanded 2.5 times that of Figure 3 for clear visualization of the relative amplitude of 〈Jãb2〉. We find that the two series of the broad pipes exist in the same way as the 〈Jãb2〉 map in Figure 3. However, there are some significant differences between them. The width of the pipe at the connection between C12 in the ring III of BPhe to the hydrogen atom in the Met residue which was abnormally large in Figure 3 is much reduced in Figure 5. The width of the two pipes which connect the two series of broad pipes are smaller in Figure 5. The width of the pipe connecting the oxygen atom bonded to C131 to the hydrogen atom of the Trp residue is also reduced considerably in Figure 5. As a result, Figure 5 represents two rather smooth streams of the tunneling currents. We also calculated the 〈Jãb2〉 map for Et ) -10.0 eV, and we obtained a quite similar map to Figure 5. For the details, see Figure S1, Supporting Information. 4.4. Average Tunneling Routes in the Elastic Tunneling Mechanism. In order to make more detailed analysis of the correlation among the interatomic tunneling currents, we term the significant interatomic tunneling currents which were represented by broad pipes in Figure 5 as shown by the red arrows in Figure 6. In this map, we used a snapshot of the tunneling currents drawn by the arrow whose width is proportional to the magnitude of Jab. Together with the significant interatomic tunneling currents, we also plotted the other interatomic tunneling currents represented by the pink arrows for comparison. It is seen that a series of tunneling currents of



Figure 6. Snapshot of the distribution of the interatomic tunneling current Jab described by the arrow with width proportional to its amplitude. Those Jab’s which satisfy |Jab|/|Jab|max > 0.022 are drawn, where |Jab|max is the maximum value of |Jab|’s. The significant interatomic tunneling currents in the 〈Jãb2〉 map in Figure 5 are represented by the red arrows and pointed by the black arrows. Those Jab’s are termed iT (i ) 1,..., 13) for the tunneling currents through the Trp residue and those Jab’s are termed iM (i ) 1,..., 8) for the tunneling currents through the Met residue. We also term the two Jab’s which flow between the Met residue and the Trp residue as 1TM and 2TM. The other interatomic tunneling currents are represented by the pink arrows.

red arrows (mainly 5M f 4M f 3M f 2M f 1M) through the Met residue flow from the QA side to the BPhe side. In contrast to this, a series of tunneling currents of red arrows (mainly 1T f 2T f 3T f 4T f 5T f 6T f 7T and 8T) through the Trp residue flow from the BPhe side to the QA side. Such opposite flows of the two series of the tunneling currents of the red arrows are typical in this ET system. There are also other broad red arrows 2TM and 7M in this snapshot of the tunneling current. It should be noted that the widths of the pink arrows are not so much narrower than those of the red arrows. As a result, the total current distribution in the map of Figure 6 is much more complex than those in the maps in Figure 3 and Figure 5. Now, we focus on the correlation among fluctuations of Jab’s. For this purpose, we calculate the mutual correlation coefficient between the two interatomic tunneling currents as follows:

rab,cd )

〈(Jab - 〈Jab 〉)(Jcd - 〈Jcd 〉)〉

√〈(Jab - 〈Jab〉)2〉√〈(Jcd - 〈Jcd〉)2〉

(24)

When Jab and Jcd fluctuate synchronously, rab, cd becomes 1 or -1. On the other hand, when Jab and Jcd fluctuate independently, rab, cd becomes 0. Therefore, the value of rab, cd represents how much synchronously the two interatomic tunneling currents fluctuate. In Figure 7, we have drawn the diagram of the absolute value of the fluctuation correlation |rab, cd| among the significant interatomic tunneling currents. Here, we treated the three interatomic tunneling currents 5M, 6M, and 7M altogether in the statistical average since positions of the three hydrogen atoms are replaced many times because of the rotation of the methyl group in the course of the simulation (see section 5). The magnitude of |rab, cd| is represented by the darkness of the red: The larger is the value of |rab, cd|, the darker is the red. We

Figure 7. Diagram of the absolute value of the mutual correlation function |rab, cd| between the important interatomic tunneling currents. The amplitude of |rab, cd| is represented by the darkness of red color.

clearly see that the group of the interatomic tunneling currents 1T, 2T, 3T, 4T, 5T, 6T, 7T, 8T, 9T, 10T, 11T, 12T, and 13T are strongly correlated to each other and are not so much correlated to the other interatomic tunneling currents except that many of them are substantially correlated with 1 TM. Similarly, the group of the interatomic tunneling currents 1M, 2M, 3M, 4M, (5+6+7)M, and 2TM are strongly correlated to each other and are not so much correlated to the other interatomic tunneling currents except that 2M and 1M are appreciably correlated with 8M. The interatomic tunneling current 8M is strongly correlated to the interatomic tunneling current 1TM. On the basis of these results, we may say that the former group of the mean-square interatomic tunneling currents forms the average tunneling route with high correlation, called Trp route. Similarly, we may say that the latter group of the mean-square interatomic tunneling currents forms the other average tunneling route with high correlation, called Met route. An appreciable correlation exists between the tunneling currents of two routes via the interatomic tunneling currents 8M and 1 TM. The numerical values of rab, cd are listed in Table S1, Table S2, and Table S3, Supporting Information. 4.5. Average Tunneling Routes in the Inelastic Tunneling Mechanism. We calculated the variance of the interatomic tunneling currents σab2. The map of σab2 is shown by the purple pipe in Figure S2. The general feature of the distribution of the σab2 with large amplitude is similar to the distribution of 〈Jab2〉 in Figure 3 even through there are some differences: The purple pipe of 1M (in the terminology of Figure 6) in the σab2 map is much narrower than the red pipe of 1M in the 〈Jab2〉 map. Next, we calculated the correction factor (1 + ηj) of eq 20, and the result is shown in Figure S3. The qualitative feature of (1 + ηj) as a function of R is similar to the correction factor (1 + δj) in Figure 4. But quantitatively, the value of (1 + ηj) is 1.5∼2.0 times larger than that of (1 + δj) in the whole range of R. Last, we calculated σab2 using eq 20. We plotted the map of σab2 in Figure 8 by the purple pipe. The pipe width in Figure 8 are expanded 1.6 times than that in Figure 5 for the purpose of clear visualization. We find that the distribution of the pipes with large amplitudes are very similar to that of Figure 5 except some points as follows: The pipes of 5M, 6M, 7M, and 1TM


Figure 8. Calculated map of the variance of the interatomic tunneling currents corrected by the quantum interference, σ˜ ab2. All of the values of σ˜ ab2 are plotted by a pipe model with purple color.

in Figure 8 are relatively broader than those of Figure 5. On the basis of these results, we can say that the similar average tunneling routes, Met route and Trp route, are formed for the ET by the inelastic tunneling mechanism as for the ET by the elastic tunneling mechanism. 5. Discussion In this study, we first calculated the mean-square interatomic tunneling current map, 〈Jab2〉 map. It may be natural to consider that the 〈Jab2〉 map can be a candidate for representing the average tunneling route in the elastic tunneling mechanism. However, the sum of p2〈Jab2〉’s passing through any surface Sj is not equal to 〈TDA2〉 since the quantum interference among Jab’s are not taken into account. Furthermore, an inconvenient thing happened; 〈Jab2〉 of the current 1M in the terminology of Figure 6 became abnormally large (see Figure 3). Analyzing the degree of the quantum interference δj among the interatomic tunneling currents (see Figure 4), we found that the quantum interference in the surface passing the current 1M is extremely destructive (δj ) -0.87). Under such situation, we preferred to make a correction of the quantum-interference effect to 〈Jab2〉. Namely, we expressed 〈TDA2〉/p2 as the sum of 〈Jãb2〉’s which are obtained by multiplying the quantum-interference correction factor (1 + δj) to 〈Jab2〉’s for the surface Sj. As a result, we obtained the 〈Jãb2〉 map in Figure 5 in which the abnormally large contribution from the current 1M is reduced. Strictly speaking, the correction factor (1 + δj) will change to some extent if we choose a surface quite different from the one we have chosen. As to this, however, we think that the surface chosen perpendicularly to the line connecting the centers of donor and acceptor will be most rational. The obtained map of 〈Jãb2〉 has a feature that the amplitude of p2〈Jãb2〉 represents the relative contribution to 〈TDA2〉. Namely, its map represents how much the tunneling current at each bond and each space contributes to the electronic factor 〈TDA2〉. From Figure 5, we find that the contributions from the interatomic tunneling currents through-spaces, 1M, 5M, 6M, 7M, 1T, 7T, and 8T, are not necessarily small or even larger than the contributions from the interatomic currents through bonds. This fact does not indicate that the through-space tunneling is more important than the through-bond tunneling in general, but it reflects a property that the total tunneling current through each surface is conserved (see eq 9) once the tunneling route is determined.

Nishioka and Kakitani Figure 5 represents that the 〈Jãb2〉’s with large amplitudes are limited to small numbers of the interatomic tunneling currents which are termed in Figure 6. By using these small numbers of the interatomic tunneling currents, the two tunneling routes, Trp route and Met route, are formed. Each route is composed of a few major pathways in the Pathways model15–17 which are fused at considerable part to each other. The major pathway is formed by the smallest numbers of bonds and through-spaces. In this sense, the present result gives a physical foundation for the Pathways model. However, the present result embraces much more profound aspect for the tunneling route as follows: (1) The interatomic tunneling currents involved in each route are more coherently coupled to each, even after taking the average of the fluctuations as shown in Figure 7. (2) The numbers of the interatomic tunneling currents involved in the route are quite limited. The other interatomic tunneling currents than the ones termed in Figure 7 do not appreciably contribute to 〈TDA2〉/p2 (see Figure 5). The reason of this feature is as follows: Originally, the magnitude of Jab of other than that of the route was relatively smaller than that of the route. By taking the square of Jab, this difference is enhanced. (3) The mean-square interatomic tunneling currents with large amplitudes in the Trp residue are largely located at one side in the two rings (see Figures 3 and 5). This asymmetrical distribution of the amplitudes is a reflection of the coherent π-electronic structure in the Trp residue. Since the obtained routes are composed of the shortest pathways connecting donor and acceptor and are almost straight, we can appreciate the reason why only the three amino acids were enough in evaluating TDA among the 828 amino acids in the RC.27,29 Furthermore, the almost straight forms of the routes suggest that the donor-acceptor distance is the most important quantity to regulate the electronic factor of the ET rate among different kinds of proteins. This result gives a support for the square-barrier model represented by the “Dutton law”.2 (We called the empirical relation of the exponential dependence of kDA on the distance the “Dutton law”.) We also see that the hydrogen atoms significantly contribute to the route only when the hydrogen atom is absolutely necessary for forming the short route as seen in the interatomic tunneling currents 1M, 5M, 6M, and 7M. We can calculate the total tunneling current flowing into quinone through the Trp residue including all of the small tunneling currents, which will be called JTrp. Similarly, we can calculate the total tunneling current flowing into quinone through the Met residue including all the small tunneling currents, which will be called JMet. Each total current contributes to TDA in the form pJTrp and pJMet. The sum rule

TDA ) p(JTrp + JMet)

(25)

is satisfied. Using the 1 000 001 data points, we produced histogram for the distributions of the values of TDA, pJTrp, and pJMet in the range from -50 × 10-5 eV to 50 × 10-5 eV. The result is plotted in Figure 9A; red box for TDA, green box for pJTrp, and blue box for pJMet. The width of the box is 0.5 × 10-5 eV. The total number of boxes is 200. The population is the number of the data points involved in the box. We found that the distribution of pJTrp is mostly located at the positive region of the current with large skewness having a long tail to the positive side. In contrast to this, the distribution of pJMet is mostly located at the negative region of the current with large skewness having a long tail to the negative side. The width of the distribution of the pJMet is smaller than that of pJTrp. The distribution of TDA is centered at around zero with broad width.


J. Phys. Chem. B, Vol. 112, No. 32, 2008 9955 In Table 1, we listed the calculated values of mean, meansquare, variances, standard deviations, and coherence parameters for the tunneling energy Et ) -9.5 eV. We see that the value of p〈JTrp〉 is positively large while the value of p〈JMet〉 is negatively large. As a result of the different sign of pJTrp and pJMet, the 〈TDA〉 is close to zero even if it is slightly positive. We also see that p2〈JTrp2〉 is about twice larger than p2〈JMet2〉 and 〈TDA2〉. Hence, it should be emphasized that 〈TDA2〉 is much smaller that the sum of p2〈JTrp2〉 and p2〈JMet2〉 even if the fluctuation of the tunneling current through the Trp route is almost independent of that through the Met route. This is because the destructive quantum interference between JTrp and JMet contributes to lower the value of 〈TDA2〉. The variance σ2 and the standard deviation σ are larger in the order of TDA, pJTrp, and pJMet. The coherence parameter C is the quantity first defined by Balabin and Onuchic for TDA as follows:27

C)

〈TDA 〉2

(28)

〈TDA2 〉

In analogy with this, we define the coherence parameter for pJTrp and pJMet as follows

CTrp ) Figure 9. (A) Histogram of the distributions of pJTrp, pJMet, TDA. The bin is 0.5 × 10-5 eV. (B) Histogram obtained by the convolution of the histograms of pJTrp and pJMet (black box) in comparison with the histogram of TDA (red box).

Its shape is almost symmetrical. We examine whether the fluctuation of pJMet is independent of that of pJTrp. For this purpose, we denote the histograms of pJTrp, pJMet, TDA as f(Xi), g(Yi), and h(Zi), where Xi, Yi, and Zi are the values of pJTrp, pJMet, and TDA at the ith box, respectively. We choose arbitrarily Xi, Yj, and Zk which satisfy the following relation

Zk ) Xi + Yj

(26)

We take the convolution of f(Xi) and g(Yj) as follows 200

F(Zk) )

∑

1 f(X )g(Zk - Xi) N i)1 i

(27)

If the histogram F(Zk) coincides with h(Zk), we can say that the two kinds of fluctuations of pJTrp and pJMet are independent to each other. We calculated F(Zk) using the histograms f(Xi) and g(Yj) in Figure 9A. The calculated histogram of F(Zk) is drawn by the black box in Figure 9B. Comparing it with h(Zk) of the red box, we find that F(Zk) is quite similar to h(Zk) except for a small difference at the peak region of the histogram: The histogram of h(Zk) is a little higher and shifted to the positive side at the peak region than that of F(Zk). This difference will be due to a small correlation of the interatomic tunneling currents between the Trp route and the Met route via the interatomic tunneling currents 8M and 1TM as we have seen in subsection 4.4. Therefore, we can say that the fluctuations of the total tunneling currents through the Trp residue and Met residue are almost independent to each other. The tunneling currents through the Trp route due to the limited numbers of the interatomic tunneling currents occupy the most part of JTrp. Similarly, the tunneling currents through the Met route occupy most part of JMet. Based on these facts, we can say that the fluctuations of the tunneling currents through the Trp route are almost independent of the tunneling currents through the Met route.

〈JTrp 〉2 〈JTrp2 〉

CMet )

〈JMet 〉2 〈JMet2 〉

(29)

The calculated values of the coherence parameters of CTrp and CMet is rather close to 1, while that of C is much smaller than 1. This result is easily understood from Figure 9A. The distributions of pJTrp and pJMet are located mostly at the positive and negative sides, respectively, and so their coherence parameter becomes close to 1. On the other hand, the center of the distribution of TDA is located at almost zero with broad width, and so its coherence parameter becomes much smaller than 1. In Table 2, we listed the calculated values of the same quantity as Table 1 for the tunneling energy Et ) -10.0 eV. We find that most parts of the results of Table 2 are similar to those of Table 1 except for one point. The value of the coherence parameter C for TDA is 0.0102 for Et ) -10.0 eV while 0.0870 for Et ) -9.5 eV. Namely, the coherence parameter C for Et ) -10.0 eV becomes about nine time smaller than that for Et ) -9.5 eV. This result is obtained because pJTrp is positively decreased a little and pJMet is negatively increased a little and then 〈TDA〉 is decreased about three time as the tunneling energy Et is decreased from -9.5 to -10.0 eV. This result indicates that the variation of Et does not affect appreciably most of the tunneling quantities, but it affects the value of C very much when the destructive interference is significant (C , 1). In other words, unless the value of Et is correctly evaluated, the calculated value of C is considered as an estimation of the order of magnitude so far as C is less than about 0.1. From Table 1 and Table 2, we find that the value of the variance σDA2 is similar to the value of 〈TDA2〉. This fact indicates that the value of the electronic factor for the ET by the inelastic tunneling mechanism is similar to that for the ET by the elastic tunneling mechanism in this ET system. We have also seen that the σ˜ ab2 map is similar to the 〈Jãb2〉 map. Then, a specific feature of the inelastic tunneling mechanism will not appear in the electronic factor, but it appears in the nuclear factor.37,38 The Met route is formed by the tunneling currents starting from the methyl group attached to C5 of quinone to the carbon atom C12 in the ring III of bacteriopheophytin. The molecular orbital (MO) coefficients of the three hydrogen atoms of the methyl group is largely regulated by the hyperconjugation effect with the π-conjugated ring of quinone. When the methyl group



TABLE 1: Calculated Values of the Various Tunneling Quantities for the Tunneling Energy Et ) -9.5 eV

TDA pJTrp pJMet

mean (× 10-4 eV)

mean-square (× 10-8 eV 2)

variance (× 10-8 eV 2)

standard deviation (× 10-4 eV)

coherence parameter

0.315 1.193 -0.878

1.139 2.201 1.157

1.040 0.778 0.386

1.020 0.882 0.622

0.0870 0.646 0.666

TABLE 2: Calculated Values of the Various Tunneling Quantities for the Tunneling Energy Et ) -10.0 eV

TDA pJTrp pJMet

mean (× 10-4 eV)

mean-square (× 10-8 eV 2)

variance (× 10-8 eV 2)

standard deviation (× 10-4 eV)

coherence parameter

0.104 1.055 -0.951

1.063 1.883 1.371

1.052 0.769 0.467

1.026 0.877 0.683

0.0102 0.592 0.660

is rotated during the MD simulation, the manner of the hyperconjugation is altered and the values of the MO coefficients are changed even in those signs, as shown in eq 23. Figure 10 shows how the rotational angle θ of the methyl group of quinone is changed with time. We confined the rotational angle θ to the range from -180° to 180°. We see that the large change of θ occurs in 1-10 ps. We also observe that θ resides at mostly 0°, 120°, and -120°. This is because the large, negative partial charge of the oxygen atom of the quinone attracts the positive partial charges of hydrogen atoms of the methyl group. It is found that the fluctuating rotation of the methyl group takes place many times within the sampling time of 1 ns. If we do not take into account the variable hyperconjugation effect due to the methyl rotation properly, we are led to quite different results. Indeed, when we calculated the pJMet due to the Met route by fixing the MO coefficients of the three hydrogen atoms of the methyl group, we obtained a result that the distribution of the values of pJMet is symmetrical with the center at zero value of pJMet. This result is quite different from the present result of Figure 9A in which the center of the distribution is largely shifted to the negative value of pJMet. Then, if the variable hyperconjugation effect is not taken into account, the distribution of TDA is shifted much to the positive value, and the coherence parameter C becomes rather close to 1. This result is quite different from the value of Table 1. Therefore, it is instructive to properly take into account the variable hyperconjugation effect when the methyl group plays a significant role in the electron tunneling. In the early stage of this study, we have conducted the MD simulations without restriction for the phytil chain of quinone.

Figure 10. Time-dependence of the rotational angle θ of the methyl group attached to C5 of the quinone ring. The origin of θ is chosen at the position when one of the C-H bonds is parallel to the quinone ring plane. The angle is confined to the region from -180° to 180°.

In the rather long run of simulations, more than 500 ps, we sometimes observed that the quinone plane abruptly slides to the side of BPhe. By carefully investigating the X-ray crystallographic structure of 1AIJ,41 we found that most of the phytil chain is exposed outside the surface of the protein. One LDAO molecule was found to be attached to the phytil chain. Together with the ditch of protein, the LDAO molecule covered the region of C7-C20 of the phytil chain to protect it from exposure to the outside of the protein. Recent X-ray crystallographic study clarified that the phytil chain is attached by two LDAO molecules51 and they protected much more region of the phytil chain from exposure to the outside of the protein. Then, we expect that many LDAO molecules will be attached to the phytil chain, and they completely protect the phytil chain from exposure to the outside of the protein. Covered by the LDAO molecules in the crystal or in solution, the motion of the phytil chain will be much restricted. In the present simulation study, we restricted the position of carbon atoms in the region C21 to C31 of the phytil chain, and we found that the sliding motion is ceased. The average position of the quinone ring was in good agreement with that obtained by the X-ray crystallographic study.41 The above experience of the simulation will give us a hint for the mechanism of how the quinone molecule is inserted into the binding site of the RC at the time of biosynthesis of the complete RC. Namely, the quinone ring will splash through the lipid molecules surrounding the protein of RC until it reaches the binding site. After the quinone ring arrives at the binding site, lipid molecules will totally cover the phytil group outside the protein surface, bringing about the stability for the motion of the phytil chain and quinone ring. It is believed that many kinds of fluctuations of protein conformation ranging from femtoseconds to seconds exist.52–54 In our present study, only the average was taken for the fluctuation of 1 ns time scale. Therefore, our average should not be called the thermal average. The conformation fluctuation of protein of the long time scale is concerned with the transition among protein conformations trapped at deep potential minima. If the ET rate is more rapid than the transition rate of the protein conformation trapped at the deep potential minima, the ET rate should be evaluated separately corresponding to each protein conformation at the deep potential minimum. This ET rate may be called the “primitive” ET rate. The thermally averaged ET rate would be obtained by averaging the “primitive” ET rates with weight of occupation probability at each deep potential minima. In the ET from BPhe- to QA, the ET rate was observed as 2 × 10-10 s.55 Even if it may be the thermally averaged ET rate, each “primitive” ET rate is very rapid. The time range of 1 ns for the average used in our study is larger than the above ET time and will be much shorter than the time of protein conformation change corresponding to the transitions among

Average Electron Tunneling Route the deep potential minima. Therefore, the averaging procedure in our study will be pertinent to the “primitive” ET in the protein conformation trapped at the global minimum, since our simulations started from the optimized conformation of the protein. It will be valuable to examine the physical reasoning of the coherence parameter C, based on the abundant data of the interatomic tunneling current obtained in the present study. As written in eq 28, the parameter C measures how much 〈TDA2〉 is larger than 〈TDA〉2 due to the fluctuation of TDA. This parameter has been used to judge whether the quantum interference among the tunneling currents through the pathways is constructive or destructive depending on the value of C:27 When C is close to 1, mostly constructive interference works. When C is much smaller than 1, mostly destructive interference works. The condition for C to be much smaller than 1 is as follows: The distribution of TDA is nearly centered at the zero value of TDA, and the distribution is broader. Then, it will be necessary to examine by what mechanism such a distribution of TDA is correlated with the quantum interference among the tunneling currents. In Figure 9, we have shown that the tunneling current JTrp through the Trp residue and the tunneling current JMet through the Met residue have broad distributions of their values even if they are expressed as the sum of the tunneling currents of constructive interference (CTrp ) 0.592 and CMet ) 0.660). This fact indicates that the broad distribution of the tunneling currents is not due to the large fluctuation of the quantum interference among the interatomic tunneling currents but is mostly due to the fluctuation of the amplitudes of the interatomic tunneling currents (let us be reminded of the data in Figure 7 that the fluctuations of the interatomic tunneling currents are more correlated to each other and proceed almost synchronously in each route). Although we do not know the origin of the fluctuation of the amplitudes of the interatomic tunneling currents at the present time, this kind of fluctuation always exists independently of the quantum-interference effect. As we see in Figure 9, the distributions of the tunneling currents in the two routes are located at the positive and negative sides of the tunneling current, respectively. This fact will indicate that the distribution of the tunneling current in each route which is made of the tunneling currents of constructive interference cannot be centered near the zero value. In general, the distribution of TDA can be centered near the zero value only when more than two routes coexist whose distributions of the tunneling currents are centered at the opposite sides of the zero value of the tunneling current. In other words, when the distribution of TDA is mostly centered at the zero value, TDA is made of more than two routes whose distributions of the tunneling currents are at the opposite sides of the zero value of the tunneling current. This fact actually indicates that the destructive interference works in the case that the distribution of TDA is centered near zero value and the coherence parameter C becomes much less than 1. Therefore, the present microscopic calculations by means of the interatomic tunneling currents provide a physical reasoning why the coherence parameter can be used as a judgment of the quantum interference among the tunneling currents. 6. Conclusion In this paper, we developed a theory for obtaining the average electron tunneling route in which the quantum-interference effect among the interatomic tunneling currents are taken into account. The novelty of the present theory is that the average tunneling route is directly related to the electronic factor 〈TDA2〉 in the elastic tunneling mechanism and σDA2 in the inelastic tunneling mechanism. As a result, the average tunneling route was

J. Phys. Chem. B, Vol. 112, No. 32, 2008 9957 expressed with use of quadratic quantities of the interatomic tunneling currents instead of linear quantities of the interatomic tunneling currents which were done before. The quadratic expression of the tunneling route using the interatomic tunneling currents exaggerates the significant interatomic tunneling currents in the route. Applying this theory to the ET reaction of BPhe- f QA in BPRC of Rhodobacter sphaeroides, we found that much destructive interference remains among the interatomic tunneling currents even after taking the average. The average tunneling route in the elastic tunneling mechanism was visualized by drawing the 〈Jãb2〉 map with use of the pipe model. Examining the 〈Jãb2〉 map, we found that the two groups of small numbers of interatomic tunneling currents which were mutually coupled with high correlation in each group have large amplitude of 〈Jãb2〉 and form the two routes, Trp route and Met route. We also found that the fluctuation of the total tunneling current through the Trp route is almost independent of that through the Met route. The average tunneling route is useful to examine the availability of the proposed models so far. For example, we notice that each of the two routes visualized by the 〈Jãb2〉 map is composed of a few major pathways in the Pathways model which are fused at considerable part to each other. This result gives a physical foundation for the reason why the Pathways model worked so well. We also found that the obtained routes are almost straight. This property will become a support for the simple square-barrier model or the “Dutton law” for the ET’s in various proteins. We also produced the σ˜ ab2 map for the electronic factor in the inelastic tunneling mechanism and obtained the similar average tunneling routes as those in the elastic tunneling mechanism. Acknowledgment. This work was supported by the Grantin-Aid on Scientific Research (C) to T.K. from the Ministry of Education, Culture, Sports, Science and Technology of Japan. Supporting Information Available: Figure of the 〈Jãb2〉 map for Et ) -10.0 eV. Map of the variance of the interatomic tunneling currents σab2. Figure of the correction factor (1 + ηj) as a function of the distance of the surface Sj from the center of donor R. Tabels of the correlation coefficients rab, cd of the significant interatomic tunneling currents in the Trp route and Met route. This material is available free of charge via the Internet at http://pubs.acs.org. Note Added after ASAP Publication. This article was published on the web on July 16, 2008. Due to a production error, some of the author’s corrections were not included. The correct version was published on July 22, 2008. References and Notes (1) McConnell, H. M. J. Chem. Phys. 1961, 35, 508. (2) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 355, 796. (3) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265. (4) Beratan, D. N.; Onuchic, J. N.; Hopfield, J. J. J. Chem. Phys. 1987, 86, 4488. (5) Newton, M. D. Chem. ReV. 1991, 91, 767. (6) Skourtis, S. S.; Beratan, D. N. AdV. Chem. Phys. 1999, 106, 377. (7) Stuchebrukhov, A. A. J. Chem. Phys. 1996, 105, 10819. (8) Stuchebrukhov, A. A. J. Chem. Phys. 1997, 107, 6495. (9) Gehlen, J. N.; Daizadeh, I.; Stuchebrukhov, A. A.; Marcus, R. A. Inorg. Chim. Acta 1996, 243, 271. (10) Okada, A.; Kikitani, T.; Inoue, J. J. Phys. Chem. 1995, 99, 2946. (11) Stuchebrukhov, A. A. AdV. Chem. Phys. 2001, 118, 1. (12) Stuchebrukhov, A. A. Theor. Chem. Acc. 2003, 110, 291. (13) Bixon, M; Jortner, J. AdV. Chem. Phys. 1999, 106, 35. (14) Winkler, J. R. Curr. Opin. Chem. Biol. 2000, 4, 192.

9958 J. Phys. Chem. B, Vol. 112, No. 32, 2008 (15) Beratan, D. N.; Betts, J. N.; Onuchic, J. N. Science 1991, 252, 1285. (16) Beratan, D. N.; Betts, J. N.; Onuchic, J. N. J. Phys. Chem. 1992, 96, 2852. (17) Regan, J. J.; Risser, S. M.; Beratan, D. N.; Onuchic, J. N. J. Phys. Chem. 1993, 97, 13083. (18) 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. (19) Regan, J. J.; Onuchic, J. N. AdV. Chem. Phys. 1999, 107, 497. (20) Kawatsu, T.; Kakitani, T.; Yamato, T. Inorg. Chim. Acta 2000, 300-302, 862. (21) Kawatsu, T.; Kakitani, T.; Yamato, T. J. Phys. Chem. B 2001, 105, 4424. (22) Skourtis, S. S; Beratan, D. N. J. Phys. Chem. B 1997, 101, 1215. (23) Page, C. C.; Moser, C. C.; Chen, X.; Dutton, P. L. Nature 1999, 402, 47. (24) Jones, M. L.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. A 2002, 106, 2002. (25) Daizadeh, I.; Medvedev, E. S.; Stuchebrukhov, A. A. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 3703. (26) Wolfgang, J.; Risser, S. M.; Priyadarshy, S.; Beratan, D. N. J. Phys. Chem. B 1997, 101, 2986. (27) Balabin, I. A.; Onuchic, J. N. Science 2000, 290, 114. (28) Kawatsu, T.; Kakitani, T.; Yamato, T. J. Phys. Chem. B 2002, 106, 11356. (29) Nishioka, H.; Kimura, A.; Yamato, T.; Kawatsu, T.; Kakitani, T. J. Phys. Chem. B 2005, 109, 1078. (30) Prytkova, T. R.; Kurnikov, I. V.; Beratan, D. N. J. Phys. Chem. B 2005, 109, 1618. (31) Prytkova, T. R.; Kurnikov, I. V.; Beratan, D. N. Science 2007, 315, 622. (32) Medvedev, E. S.; Stuchebrukhov, A. A. J. Chem. Phys. 1997, 107, 3821. (33) Tang, J. J. Chem. Phys. 1993, 98, 6263. (34) Bixon, M.; Jortner, J. Russian Journal of Electronchemistry. 2003, 39, 3. (35) Troisi, A.; Nitzan, A.; Ratner, M. A. J. Chem. Phys. 2003, 119, 5782. (36) Skourtis, S. S.; Balabin, I. A.; Kawatsu, T.; Beratan, D. N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 3552. (37) Nishioka, H.; Kimura, A.; Yamato, T.; Kawatsu, T.; Kakitani, T. J. Phys. Chem. B 2005, 109, 15621. (38) Nishioka, H.; Yamato, T.; Kakitani, T. Mol. Simul. 2006, 32, 727. (39) Katz, D. J.; Stuchebrukhov, A. A. J. Chem. Phys. 1998, 109, 4960.

Nishioka and Kakitani (40) Kawatsu, T.; Kakitani, T.; Yamato, T. J. Phys. Chem. B 2002, 106, 5068. (41) Stowell, M. H. B.; McPhillips, T. M.; Rees, D. C.; Soltis, S. M.; Abresch, E.; Feher, G. Science 1997, 276, 812. (42) Morikami, K.; Nakai, T.; Kidera, A.; Saito, M.; Nakamura, H. Comput. Chem. 1992, 16, 243 PRESTO(A vectorized molecular mechanics program for biopolymers) (43) 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. (44) Saito, M. Mol. Simul. 1992, 8, 321. (45) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (46) Howell, J.; Rossi, A.; Wallace, D.; Haraki, K.; Hoffmann, R. QCPE 1977, 11, 344 FORTICON8 (Extended Hu¨ckel method program) (47) Kincaid, D. R.; Respess, J. R.; Young, D. M.; Grimes, R. G. ITPACK 2C; University of Texas: Austin, TX, 1999 (a fortran package for solving large sparse linear system by adaptive accelerated iterative methods). (48) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 209. (49) Frisch, J.; Trucks, G. W.; Schlegel, H. B.; Scuseria G. E.; Robb, M. A.; Cheeseman,J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Jr. Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian98, revisionA.9; Gaussian, Inc.; Pittsburgh, PA, 1998. (50) Anderson, A. B.; Grantscharova, E. J. Phys. Chem. 1995, 99, 9149. (51) Koepke, J.; Krammer, E.-M.; Klingen, A. R.; Sebban, P.; Ullmann, G. M.; Fritzsch, G. J. Mol. Biol. 2007, 371, 396. (52) Austin, R. H.; Beeson, K. W.; Eisenstein, L.; Frauenfelder, H.; Gunsalus, I. C. Biochemistry 1975, 14, 5355. (53) Frauenfelder, H.; Parak, F.; Young, R. D. Annu. ReV. Biophys. Biophys. Chem. 1988, 17, 451. (54) Yang, H.; Luo, G.; Karnchanaphanurach, P.; Louie, T-M.; Rech, I.; Cova, S.; Xun, L.; Xie, X. S. Science 2003, 302, 262. (55) Gunner, M. R.; Dutton, L. J. Am. Chem. Soc. 1989, 111, 3400. (56) Humphrey, W. F.; Dalke, A.; Schulten, K. J. J. Mol. Graphics 1996, 14, 33.

JP710689S

Average Electron Tunneling Route of the Electron Transfer in Protein

Recommend Documents