Protein charge transfer far from equilibrium: a theoretical perspective

Abstract

Potential differences for protein-assisted electron transfer across lipid bilayers or in bio-nano setups can amount to several 100 mV; they lie far outside the range of linear response theory. We describe these situations by Pauli-master equations that are based on Marcus theory of charge transfer between self-trapped electrons and that obey Kirchhoff's current law. In addition, we take on-site blockade effects and a full non-linear response of the local potentials into account. We present analytical and numerical current–potential curves and electron populations for multi-site model systems and biological electron transfer chains. Based on these, we provide empirical rules for electron populations and chemical potentials along the chain. The Pauli-master mean-field results are validated by kinetic Monte Carlo simulations. We briefly discuss the biochemical and evolutionary aspects of our findings.

---

1. Introduction

Electron transfer reactions can be found throughout life on earth. Photosynthesis, cellular respiration, sensing and DNA repair are processes that are vital to many, and for some of these processes to all organisms.^1–3 Biological electron transfer is based upon centers of excess electron localization, such as hemes or iron–sulphur clusters. Frequently, the proteins hosting the cofactors are attached to lipid bilayers or even bridge entire biological membranes. For example, the respiratory complex I is partly embedded in the inner mitochondrial or bacterial membrane^4–6 The photoreaction reaction center of Rps. viridis is a membrane protein that contains a cytochrome subunit with four hemes, exposed to an aqueous environment.⁷ Electron transfer involving this cytochrome in recharging the chlorophyll special pair has been addressed by Dohse et al.⁸ and by Medvedev et al.⁹

Within a cycle of work, both of these enzymes transfer two electrons. The biochemically controlled transfer of electrons is not limited to pairs: the NrfHA nitrite reductase catalyses the reduction of nitrite to an ammonium cation, a reaction that involves the transfer of six electrons. In its minimum functional unit, NrfHA_n (n = 1 or 2) contains nine to 14 heme groups, and the enzyme is attached to the inner membrane of bacteria by a hydrophobic anchor.¹⁰

The structural biology of membrane proteins has seen considerable progress in recent years, in particular due to the availability of high-resolution cryo transmission electron microscopy.^11,12 In addition, the combination of redox titration experiments and in situ spectroscopy has added to insight into the energetics of electron transfer reactions, often down to pinning individual transfer steps and the redox states of single electron acceptors. Usually, these studies are performed at thermodynamic equilibrium in aqueous solution, while the biologically operative protein complexes are associated with membranes and experience a redox potential gradient.¹³ For example, electron transfer in the respiratory complex I is driven by the standard redox potential difference of at least 520 mV between NAD⁺/NADH and 1,4-naphthoquinone-based mobile redox carriers.^14,15 For NrfH, we find a standard potential difference of 650 mV between the nitrite and the ammonium ion.¹⁶ For light-induced processes, potential differences span a considerable fraction of the spectrum of visible light.

The situations under review in this paper are depicted in Fig. 1 as cartoon models. Fig. 1a shows a redox-active protein at equilibrium, the redox states of the excess electron centers of localization can be populated by tuning an external chemical potential, μ_red. In Fig. 1c, the protein is embedded in a membrane and subject to two chemical potentials, to which we will refer to as μ_L and μ_R. Such membranes, e.g. lipid bilayers, separate cell or bacterial compartments in which different physical and chemical conditions can be realized. In the case studied here, this compartmentalization is essential for realizing different chemical potentials on either side of the bilayer, their difference drives the electron transport process. These chemical potentials depend on the standard redox potential and on the concentration of the redox-active species on either side of the membrane. The major aim of our work is the computation of the charge current, the charge carrier populations and the chemical potentials of the electron acceptors as a function of μ_L, μ_R and the energies characteristic of biological charge transfer. The situation of a protein bridging two compartments is related to that of a biopolymer spanning metallic contacts or leads, as shown in Fig. 1b. The major physical difference to Fig. 1c is the presence of a continuum of states above and below μ_L and μ_R rather than a single state representing the chemical potential.

Fig. 1 Cartoon representation of protein charge transfer scenarios discussed in this work. Redox-active protein in equilibrium with its environment, resembling a redox titration (a), non-equilibrium transfer through a protein bridging two metallic leads (b) or a lipid bilayer (c). Protein atoms as grey spheres, chromophore atoms in color.

The setup shown in Fig. 1b is reminiscent of that found in the field of molecular electronics, i.e. small- to medium-sized molecules spanning a small gap or bridging a scanning tunneling microscope tip and a surface.¹⁷ In nanostructures such as molecular junctions conductivity is often accompanied by strong coupling of the transport electrons to the mechanical degrees of freedom such as local vibrational modes or phonons.^18,19 For higher bias voltages, i.e., far from equilibrium, this coupling gives rise to non-conservative forces and unusual dissipative phenomena, which have important consequences. For example, theoretical studies have predicted that in situations far from equilibrium current-induced forces may cause so-called runaway vibrational modes. Furthermore, negative friction may exist, both resulting in mechanical instability of the structures. Theoretical studies of these phenomena have, for example, been based on quantum mechanical rate theories,²⁰ classical concepts such as Langevin equations^21,22 or the hierarchical equations of motion approach.^19,23 The latter is a density matrix approach that provides a fully quantum mechanical framework to study nonequilibrium transport in open quantum systems.

There are, however, major differences between molecular junctions and biopolymers, which render the application of the schemes listed above difficult from a numerical perspective. Proteins and DNA exist in aqueous environment, and the negative DNA glycophosphate backbone charges have to be compensated. As a consequence, these systems are dominated by strong electron–phonon coupling, emerging from a continuous vibrational density of states that ranges down to the Terahertz range of the electromagnetic spectrum. As detailed below, this coupling is expressed in a single characteristic quantity, the reorganization energy λ. For a single site that bridges two metallic contacts, models of this type have been studied by Sowa et al.^24,25 In our paper, we discuss both electron transport chains coupled to single entry and exit levels and chains coupled to metallic leads, as they are intimately related mathematically and computationally.

As a biopolymer that exhibits long-range electron transfer, DNA has been in the spotlight for many years. Two different types of experiments have been in the focus of research. Light-induced DNA fragmentation has been used in the work of Barton, Giese, Schuster, Wagenknecht, to name a few.^26–29 Here, a transition metal complex intercalated into DNA induces charge separation upon irradiation, and an excess hole can travel along the DNA strand until it is captured by a (GC)₃ moiety, with consecutive oxidative fission in this place. DNA amplification and consecutive chromatography has revealed a notably large length scale of charge hopping up to ∼200 Å. In a second type of experiments, short DNA strands bridge a metal surface and a scanning electron microscope tip, permitting the observation of current–voltage curves. Upon the application of a voltage drop of 2–4 V along the strands, conductivities amount to some 10 nA.³⁰ Both processes have been modeled theoretically and computationally.^31–33 While for the latter, the current–voltage curves have been reproduced quantitatively, approximations have been made, such as postulating a linear voltage drop along the chain, or restrictions upon the total number of charge carriers that populate the chain. With application to a different kind of biopolymers, electron transfer proteins, we present a scheme that overcomes these constraints and approximations.

The remaining part of this work is organized as follows. In the next section, we present the model, which is based upon hopping conduction in a polarizing environment. We also detail the computational approach that is used to solve the resulting Pauli-master equations at constant flow, and a kinetic Monte Carlo approach that introduces a stochastic element to the description of the transport process. In the third section, analytical results are given for a specific set of parameters. In addition, computational results are presented in the third section for setups that model proteins spanning membranes and nanoscopic setups. In the final section, the results are discussed in a general context, and conclusions are derived.

2. Model and methods

2.1. Model background and Pauli-master equations

A description of nonequilibrium transport between centers of excess electron trapping has to make at least some reference to the energetics of Marcus' theory of electron transfer,^34,35 to which it can be reduced at equilibrium. The three relevant parameters are the electronic coupling between the centers, t, the driving force, ΔG, and the reorganization energy, λ. For two centers at equilibrium, the corresponding nonadiabatic transfer rate is given by:In the following, we will abbreviate the preexponential factor as k⁰. For charge transfer in proteins, reorganization energies typically amount to ∼1 eV,^36–38 an energy scale that has to be compared to the thermal energy of 25 meV at room temperature. As a consequence, electron transfer is thermally activated and can be interpreted as small polaron hopping from the perspective of solid state physics.^39,40 For clarity only, the reduction of the activation barrier by the electronic coupling, t, is not considered here. In the range of parameters used here, this reduction is not relevant. For biological systems like the NrfH protein studied here, we typically have reorganization energies of 0.5 to 1 eV and couplings less than 50 meV.⁴²

Away from equilibrium, we have to take into account the drop in electronic chemical potential (or the voltage) along the transport chain, which modifies the effective local driving force and leads to

for sites located within the chain, i.e. those not connected to leads. We have now indexed the sites participating in the transfer step and note the symmetries k⁰_ij = k_ji⁰, ΔG_ij = −ΔG_ji, λ_ij = λ_ji and Δμ_ij = −Δμ_ji.

As the centers of localization accommodate excess charge on a molecular scale, they exhibit a strong on-site electron–electron repulsion, usually limiting the number of excess electrons or holes to a single one. From a mean-field perspective, this local blockade effect can be taken into account using the local excess charge carrier population p_i. From this point of view, charge transfer from a site i to a site j occurs with a probability proportional to p_i(1 − p_j). In the end, we will drop the mean-field character of the populations p and simulate charge transfer between states characterized by binary integer charges using a kinetic Monte Carlo algorithm.

All the elements described above can be arranged into a system of Pauli-master equations:

We investigate a linear arrangement of sites for convenience, as depicted in Fig. 2, and we note in the remaining part of the paper whenever this affects the physical properties of the system. We consider a chain of redox-active sites at constant flow, tantamount to ṗ_i = 0 for all sites i. For a quantity that is conserved, e.g. the charge, the net flow into each site from the left (positive contributions to the rhs of eqn (3)) and to the right (negative contributions to the rhs of eqn (3)) are equal. As a consequence, we have

k_i−1,ip_i−1(1 − p_i) − k_i,i−1p_i(1 − p_i−1) = k_i,i+1p_i(1 − p_i+1) − k_i+1,ip_i+1(1 − p_i). (4)

Using the decomposition indicated in eqn (2), we writewhich we abbreviate as A_i − B_i = C_i − D_i in the sequence of the terms in eqn (5).

Fig. 2 Conductivity model used in this work. A chain of n sites 1,…,i, j,…,n is connected to a left (L) and right (R) lead, which exhibit a constant potential. Hopping rates k_ij within the system, and k_L, k_R, to and from the leads. The leads may either be metallic junctions (Fig. 1b) or consist of single levels (Fig. 1c).

For a chain coupled to a left and a right single-level redox potential μ_L and μ_R (the situation depicted in Fig. 1c), the following equations hold for the entry and the exit site. For the first site, we have

In a similar way, we can write for the final siteWe now discuss the presence of metallic leads, where the density of electronic states is continuous (Fig. 1b). In electrochemistry, this model is known as the Marcus-DOS or the Marcus–Hush–Chidsey theory.⁴¹ We make use of the approximation of a constant density of states and a Fermi distribution, F_FD(ε), that is represented by a Heavyside step function of the electronic energy, implying that F_FD(ε) = Θ(μ_L − ε) at the left and F_FD(ε) = Θ(μ_R − ε) at the right electrode.

For the left lead, we now have to modify the f for the entry process

andfor the exit process. Similar relations hold for the right electrode. We only apply the zero-temperature approximation to the Fermi distribution of the lead electrons, and we believe that the error introduced in this way is much smaller than that of assuming a constant lead electronic density of states.

For small populations, small reorganization energies and small chemical potential differences (the latter two compared to k_BT), we recover Ohmic conductivity and Kirchhoff's current law by transforming the chemical potentials μ into voltages U,

and the quantities k⁰_ij become conductances. The accuracy of a linear approximation will be discussed below in the context of eqn (14), which has the same series expansion.

Detailed balance is satisfied under those conditions where Boltzmann statistics can be applied, i.e. in the absence of external driving forces and for a number of particles that is much smaller than the number of states, i.e. p_i ≪ 1.

2.2. Numerical approach

To solve the Pauli-master equations at constant flow, we have to fulfill A_i − B_i − C_i + D_i = 0 for each site, including the corresponding equations for the entry and the exit site. For non-reversible models, as appropriate for high differences between μ_L and μ_R, A_i − C_i ≃ 0 holds.

A _i , B_i, C_i and D_i depend on the n populations, which can be assembled into a vector, p⃑, and on the n local chemical potentials, [small mu, Greek, vector]. The resulting system of n equations is nonlinear, and it is subject to the constraints 0 ≤ p_i ≤ 1 ∀i. For vanishing equilibrium driving forces, ΔG_ij = 0, and a linear topology, μ_L ≥ μ_i ≥ μ_R∀i holds. In this situation, the chemical potentials also have to obey μ_i+1 ≤ μ_i. For a given [small mu, Greek, vector], we have

A_i − B_i − C_i + D_i = 0 = g_i(p⃑), (11)

which can be solved iteratively by Newton's method with a reasonable initial guess p⃑₀,

p⃑_N+1 = p⃑_N − J⁻¹g(p⃑_N). (12)

Here, N is the index of the iteration step, J is the Jacobian matrix, which contains the partial derivatives of the rows g_i with respect to p_j, J_ij = ∂g_i/∂p_j, which can be evaluated analytically or numerically. The eqn (12) can be written as a linear system of equations,

J(p⃑_N+1 − p⃑_N) = −g(p⃑_N), (13)

that can be solved at each iteration step. As the populations enter the eqn (5)–(7) linearly and quadratically, only a small number of iterations (typically less than ten) are required to achieve convergence to ten digits. To address the problem of computing the chemical potentials, we require that the total current through the system becomes a maximum. We note that this requirement has to hold strictly in the Ohmic regime,⁴³ and otherwise can be interpreted as an upper limit to the current. The current is computed as A_i − B_i times the elementary charge for an arbitrary site i. As we consider both the rate equations and look for a maximum of the conductivity, the problem is not underdetermined any more, and we find unique solutions. We have probed the uniqueness by (i) always recovering the analytical solutions numerically with randomized initial conditions and (ii) using different initial conditions for the cases that can only be studied numerically.

To solve for p⃑ and [small mu, Greek, vector] simultaneously, we follow the strategy of Koslowski and Wilkening applied in the Ohmic regime with an added local blockade and backward transfer.⁴⁴ Starting from an initial guess for both quantities, the μ_i are subject to a Monte Carlo (MC) optimization. In each MC step, a random site i is selected, and the corresponding chemical potential μ_i is altered by ±Δμ_i, which is drawn from a uniform distribution. The change is accepted if μ_i fulfills the constraints listed above – if applicable – and if the conductivity increases. For the old and the new μ_i, we solve for p⃑_old and p⃑_new by Newton's method. Eqn (13) is solved by the dgesl and dgefa subroutines of the linpack library.⁴⁵ For small system sizes, the potentials may also be scanned by n coarse-grained outer loops in a brute force manner.

2.3. Kinetic Monte Carlo simulations

We apply a kinetic Monte Carlo simulation in order to overcome the mean-field approximation to the site populations, p_i, which can be continuous, 0 ≤ p_i ≤ 1. In the following, we describe the local populations by a binary (or spin-like) variable, s_i, that can only exhibit values of zero or unity, as appropriate to strong electron trapping in the regime of large electron–phonon coupling. Each electronic state can be represented by a sequence of such variables, S_m = {s₁, s₂, s₃,…, s_n}. For example, S₆ = {0, 1, 1, 0} = {s₁, s₂, s₃, s₄} defines a state with four sites populated by two electrons at sites 2 and 3. The value of m is computed by accumulating the powers of two, weighted by s_i, starting from the left: . We take into account all possible one-electron hops between the sites along the chain and the two reservoirs. The sample state S₆ can turn into S₁₄ = {0, 1, 1, 1} by accepting an electron from the right lead, or into S₇ = {1, 1, 1, 0} by the same process involving the left lead. These processes occur with a rate of k_R and k_L, respectively. In addition, S₆ can turn into S₅ = {1, 0, 1, 0} by electron hopping from site 2 to site 1 with a rate k₂₁ or into S₁₀ = {0, 1, 0, 1} by electron hopping from site 3 to site 4 with a rate k₃₄.

We make use of a rejection Monte Carlo algorithm,^46,47 starting with a randomly populated system at t = 0, and we calculate its index m. We choose one of the states [small script l] out of L accessible by m, and compute the rate k_{m[small script l]} by identifying the corresponding electron hopping process and its rate. The move is accepted with a probability k_{m[small script l]}/k_max, where k_max is the maximum of all hopping rates. We compute a second random number 0 < r ≤ 1 and increment the time by Δt = −(Lk_max)⁻¹ ln r. We loop through the process until a time limit is reached or a predefined number of electrons has passed through the system.

3. Results

3.1. Analytical solutions

To discuss the analytical solutions, we consider electron transfer chains connected to single levels, as sketched in Fig. 1c. For uniform, site-independent values of ΔG and λ (including the entry and exit processes) and for k⁰_L = 2k⁰₁₂ = ⋯ = 2k⁰_i,i+1 = ⋯ = 2k⁰_n−1,n = k⁰_R, it is straightforward to write down an analytical solution of the eqn (5)–(7). We find p_i = 1/2 for all sites i and a linear chemical potential drop along the system, μ_i = μ_L − iΔ[small mu, Greek, tilde], with Δ[small mu, Greek, tilde] = Δμ/(n + 1) = (μ_L − μ_R)/(n + 1). We use this solution to illustrate characteristic features of the current-chemical potential curves and their dependence on one of the characteristic energies of Marcus' theory, the reorganization energy λ. The current amounts toThe nature of the solutions becomes particularly clear for values of Δμ that are considerably larger than those usually applied experimentally due to electrochemical reactions. They are illustrated in Fig. 3 for λ = 0.75 eV and two different numbers of excess electron centers. As a function of the external chemical potential drop, we find maxima of the modulus of the current at Δμ ≃ ±(n + 1)λ. The presence of extrema is a consequence of Marcus' theory as the basis of local transport rates, it reflects the crossover from the so-called normal to the inverted regime. This occurs for a vanishing exponent at equilibrium (eqn (1)), or for a vanishing exponent of one of the two summands on the right hand side of eqn (14), first line, while the other is negligible. It is interesting to note that the maximum is shifted towards higher chemical potentials with an increasing number of centers of electron localization, n (vice versa for the minimum).

Fig. 3 Analytical solution (eqn (14)) for the current as a function of the external drop in the chemical potential for a reorganization energy of λ = 0.75 eV and n = 1 (red, solid line) and n = 3 (black, solid line) centers of excess electron localization. Open circles, ○, depict the results of kinetic Monte Carlo simulations for n = 3. The underlying model is shown in Fig. 1c. We also present numerical results for metallic leads and n = 3, as depicted in Fig. 1b. The results are shown as a black chain line.

It is interesting to compare these analytical results obtained for systems coupled to two redox potentials (Fig. 1c) to numerical values computed for those attached to metallic leads (Fig. 1b). For identical values of k⁰ and λ as used above, the current-chemical potential curve is shown in Fig. 3 for a system consisting of three chromophores. As compared to the Marcus result, the system attached to conducting leads is slightly shifted to higher values of Δμ, and it saturates slightly above Δμ = λ(n + 1) to a maximum current. This behaviour is due to the presence of metallic states in a uniform density of states below and above the Fermi level, which are always available as a source or drain of electrons.

Eqn (14) can be expanded into a Taylor series around Δμ = 0. With λ = 1 eV, at room temperature (k_BT = 0.025 eV) and with k₀ = 0.25 (in arbitrary units), we have two leading terms of 20 exp(−10)Δ[small mu, Greek, tilde] and 3400 exp(−10)Δ[small mu, Greek, tilde]³/3. For n = 1, a linearization gives rise to a relative error that is less than 10% for Δμ < 80 meV. For uniform chains, this range becomes more extended with increasing n. This energy scale has to be compared to ±2 V that can be achieved in a bio-nano setup or to the biologically relevant differences in redox potentials of up to 500 meV referenced above.

3.2. Numerical solutions

Nonuniform intrasite transfer rates k⁰_ij can be addressed numerically, and we use these computations in simple systems to test our computational methods and to check the results against physical concepts and our physical intuition. Following these calculations, we compute the charge transfer characteristics for two bacterial membrane-bound proteins, each containing four heme molecules as centers of excess charge localization.

For a two-site system, we have three rates, of which we set two to k⁰ and use the third (k⁰_L, k⁰_R or k⁰₁₂) as a bottleneck by decreasing the rate by a factor of 100. We also study an additional symmetric system with a small but equal entry and exit hopping rate. We apply external chemical potential differences of 0.5 and 1 eV, the former defining a maximum of biological ground state protein electron transfer, the latter lying in the range of nanoscopic setups involving narrow-band metals or doped semiconductors. In Table 1, we list the populations and the chemical potentials calculated for these values of Δμ. All of these calculations refer to the model with single level contacts depicted in Fig. 1c.

Table 1 Populations p_i and chemical potentials μ_i for two-site models showing one or two bottlenecks. For the systems, L and R refer to left and right leads, the numbers refer to the sites. The underlying model is depicted in Fig. 1c and consists of an electron transfer chain coupled to single levels. Hopping rates k⁰ are indicated by a line (−), reduced hopping rates of k⁰ × 10⁻² by a dot (·). Note that μ_R is set to zero

System	μ L (eV)	p 1	p 2	μ 1 (eV)	μ 2 (eV)
L-1·2-R	0.5	0.66	0.34	0.43	0.07
L-1·2-R	1.0	0.77	0.23	0.80	0.20
L·1-2-R	0.5	0.38	0.48	0.16	0.06
L·1-2-R	1.0	0.26	0.46	0.42	0.17
L-1-2·R	0.5	0.52	0.62	0.44	0.34
L-1-2·R	1.0	0.54	0.74	0.83	0.58
L·1-2·R	0.5	0.40	0.60	0.27	0.23
L·1-2·R	1.0	0.36	0.64	0.58	0.42

---

Regardless of the Δμ value applied, we always find the largest drops in the chemical potential between the sites – including the leads – that exhibit the smallest hopping rates. The hopping rates of the first and the last model listed in Table 1 are symmetric w.r.t. to an exchange of the hopping rates, this symmetry is reflected both in the populations and chemical potentials. Populations add up to unity for the symmetric models, and the individual contributions can deviate considerably from the analytical solutions for uniform models with p = 1/2, as computed above. This deviations increase with an increasing Δμ. It is important to note that nonuniform hopping rates lead to changes both in the populations p_i and in the chemical potentials μ_i, as compared to the uniform case and its analytical solution. In retrospect, this behaviour justifies the strategy of treating both quantities as variables, and of designing a scheme that permits their simultaneous computation.

We have also tested our numerical approach for two charge proteins of biological relevance. Our choice of the proteins has been motivated by our previous experience with the systems and the availability of reliable computed energy landscapes for electron transfer. Furthermore, we believe that NrfH is a very suitable candidate for a protein in bio-nano setup. It has a membrane anchor that can be attached to a lipid covering a conducting metal surface on one end, and a free valence on the heme iron on the other end that can be contacted e.g. by a His–Ni moiety.

First, we turn to an electron transfer protein that has been mentioned in the introduction, the cytochrome subunit of the photoreaction center of the purple bacterium Rps. viridis. This protein contains four heme groups, and its geometry is displayed in Fig. 4 using the X-ray structure of Deisenhofer et al.⁷ Following the approach of Gnandt et al.,³⁸ the driving forces and the reorganization energies have been computed by a numerical solution of the Poisson–Boltzmann equation applying the Delphi program package.^48,49 For the atomic charges and radii that enter these computations, we have used the values of the Amber force field.^50,51 Within the protein and the lipid bilayer, a dielectric constant of ε = 3.5 has been used; for water, we take the room temperature value ε = 78.4. Rather than computing the electronic coupling that enters eqn (1), we apply the empirical Dutton–Moser rule^36,52 to estimate the value of k⁰. The rule relates k⁰ to the edge-to-edge distance between two chromophores, r, via

k⁰(r) = k⁰(r₀)exp(−γ(r − r₀)) (15)

with r₀ = 3.6 Å, γ = 0.6 Å⁻¹ and k⁰(r₀) = 10¹³ s⁻¹. The computed input parameters of our numerical scheme are listed in Tables 2 and 3.

Fig. 4 Protein structures used in the work. Amino acid atoms drawn as grey spheres, non-iron heme atoms as red sticks, iron atoms as orange spheres. (a) Cytochrome subunit of the photoreaction center of Rps. viridis after the X-ray structure of Deisenhofer et al.,⁷ protein databank identifier 1PRC, (b) NrfH subunit of the nitrite reductase of D. vulgaris after the X-ray structure of Rodrigues et al.,¹⁰ protein databank identifier 2J7A. The hemes H1 and 335 are closest to the membrane, respectively.

Table 2 Exponential prefactor k⁰, as computed using the empirical Dutton–Moser rule eqn (15), and the characteristic energies of Marcus theory of charge transfer (thermodynamic driving force ΔG and reorganization energy λ) for the cytochrome subunit of the photoreaction center of Rps. viridis. For the enumeration of the hemes see Fig. 4a

Reaction	k ^0 (s^−1)	ΔG (eV)	λ (eV)
L → 333	8.00 × 10^10	0.00	1.0
333 → 334	7.45 × 10^10	0.05	0.81
334 → 336	1.84 × 10^10	−0.27	0.72
336 → 335	8.57 × 10^10	−0.28	0.82
335 → R	8.00 × 10^10	0.00	1.0

---

Table 3 Exponential prefactor k⁰, as computed using the empirical Dutton–Moser rule eqn (15), and the characteristic energies of Marcus theory of charge transfer (thermodynamic driving force ΔG and reorganization energy λ) for the NrfH subunit of the nitrite reductase of D. vulgaris. For the enumeration of the hemes see Fig. 4b

Reaction	k ^0 (s^−1)	ΔG (eV)	λ (eV)
L → H1	1.00 × 10^13	0.00	1.0
H1 → H2	1.32 × 10^13	−0.61	1.37
H2 → H3	4.60 × 10^11	0.00	1.3
H3 → H4	4.32 × 10^12	0.07	1.64
H4 → R	1.00 × 10^13	0.00	1.0

---

We assume that the parameters relevant for the contacts of proteins to metallic leads are essentially governed by the same relations as those relevant to charge transfer within the protein, i.e. a dependence of the tunnel splitting upon the contact-chromophore distance that follows an exponential Dutton–Moser law, a reorganization energy that is determined by the dielectric response of the aqueous environment of the system, and with driving forces that can be tuned by the material of the contacts and the external voltages.

For the contacts, we use ΔG_L = ΔG_R = 0, λ = 1 eV and a k⁰_L entry and a k⁰_R exit hopping rate close to the maximum of k⁰_ij within the chain.

For the cytochrome subunit, we have a monotonous drop of ΔG along the transport chain, a nearly uniform λ and a hopping rate bottleneck between the hemes 334 and 336 as the parameters entering the computation of the currents, the populations and the local chemical potentials as a function of Δμ. Within the computations, μ_R is fixed and defines the zero of the chemical potential, while μ_L is varied. We note that the driving forces depend on the protonation pattern of the cytochrome,⁵³ and that the charge flow in the biological system is still unknown due to the lack of the identification of the electron donor binding site.

The current-chemical potential curves are displayed in Fig. 5 for single-level (cf.Fig. 1c) and metallic (Fig. 1b) contacts. As above, the curve in the positive Δμ branch is shifted to slightly higher chemical potentials for a computation involving metallic leads. The curves are asymmetric w.r.t. Δμ = 0 and exhibit a current that is close to zero in range from −1.5 eV to 0.2 eV. This behaviour partly reflects the extent of the original free energy surface (cf.Table 2): barriers of 0.05 eV (transfer to the left) or 0.65 eV (transfer to the right) have to be overcome to facilitate electron transfer.

Fig. 5 Currents I as a function of the applied difference in the external chemical potential, Δμ for two four-heme protein subunits. Red, solid line: NrfH, single-level contacts; black, solid line: photoreaction center cytochrome, single-level contacts; (both as in Fig. 1c); red, dashed line: NrfH, metallic leads; black, dashed line: photoreaction center, metallic leads (the latter two as in Fig. 1b). Note that for the latter curves, the current has been multiplied by a factor of 100.

Similar curves have been obtained for electron transfer through the NrfH protein, which also constitutes a charge transfer process that is essentially downhill. They are also shown in Fig. 5. Compared to the cytochrome, the currents are about an order of magnitude larger throughout the entire range of Δμ. Here, larger λ values are overcompensated by a tighter packing of the chromophores, which in turn leads to larger values of k⁰. We note that the Amber charges usually lead to an overestimate of the reorganization energies⁵⁴ and that we have not applied the so-called Pekar factor that reduces the reorganization energies. At Δμ = 1 eV, the currents lie in the range of 100 to 500 pA, which puts them into the range of scanning tunnel microscope or patch clamp experiments. We note that currents of 100 pA can be detected with an accuracy of 1 ppm using commercially available instrumentation.⁵⁵ For biological soft-matter systems, ionic currents around 1 pA can be measured using patch clamp techniques,^56,57 electronic currents recorded by scanning tunnel microscope experiments permit a resolution much better than 1 nA.³⁰

While the currents as a function Δμ are as continuous as the local chemical potentials, the local populations show a different behaviour. For NrfH, they are displayed in Fig. 6. Approaching the range of a particularly small conductivity from the left, p₁ and p₄ show values in the range of 0.5 to 0.85, while p₂ and p₃ gradually converge to unity and hence block the conductivity channel by effectively taking out these sites as electron acceptors, cf.eqn (3)–(6). At zero voltage, p₂ starts to decline, and the channel is unblocked with p₁, p₃ and p₄ not far away from one half. We find a very similar behaviour both for single-site and metallic contacts, and overlapping curves reflecting these situations are not shown in Fig. 6.

Fig. 6 Excess electron populations for the four heme molecules of the NrfH nitrite reductase protein chain as a function of the external chemical potential, Δμ. Red open symbols and line: H1, single-level contacts, black open symbols: H2, single-level contacts, blue open symbols: H3, single-level contacts, green open symbols: H4, single-level contacts, the underlying model is depicted in Fig. 1c. Full symbols (•) refer to coupling to metallic leads, as shown in Fig. 1b.

3.3. Kinetic Monte Carlo

Kinetic Monte Carlo simulations^46,47 have been performed for chain lengths of n = 2 to n = 4 for the single level coupling case depicted in Fig. 1c. The parameters are identical to those situations where an analytical mean-field solution is available, which gives average local populations p_i = 0.5, a linear chemical potential drop along the chain and currents expressed by eqn (14), all as a function of the externally applied chemical potential, Δμ. For each value of Δμ, we have used the chemical potentials of the analytical solutions, but have made the populations a quantity subject to a random hopping process using the rates expressed in eqn (2). For each chemical potential difference, we have computed 10⁵ Monte Carlo steps, which have shown an acceptance ratio of roughly one half. For two sites, a net of about 4000 electrons passed the system during the simulation, reduced to ∼2000 for n = 3. The emerging currents – number of particles that passed the system divided by the total simulation time – are virtually indistinguishable from the mean-field results, they are shown in Fig. 3.

In this way, we are not only able to show that the mean-field solutions for the populations are identical to their numerical mean-field counterparts, but recover the underlying microscopic picture. The objects of hopping transport are electrons, and sites of excess electron localization exhibit a charge of zero or the elementary charge. As a simple argument from dielectric theory suggests that the reorganization energy is proportional to the square of the charge,^34,35 a hopping electron experiences the full reorganization energy λ rather than λ/4, as a mean-field value of p = 1/2 suggests.

We are not aware of a rigorous deviation of the nonlinear master equations from an extended linear scheme. As we obtain the same numerical results for a linear Fock space representation as for the mean-field nonlinear site representation (cf.Fig. 3), we are confident that the introduction of the nonlinear terms is valid.

4. Conclusions

In this work, we have presented a model of steady-state electron transfer in a strongly polarizable medium far from equilibrium. It extends the equilibrium charge transfer theory of Marcus, and it applies to biological charge transfer in DNA or proteins. Focussing on the latter, we use the characteristic energies of Marcus' theory as input parameters to a system of Pauli-master equations for hopping transport between centers of excess electron localization. The electrons are subject to an external chemical potential difference via metallic leads or single levels. They experience local blockade effects and obey Kirchhoff's current law. In the parameter space addressed in our work, we have always checked that the solution has a steady-state character, both in the site representation (mean-field populations) and in Fock space (kinetic Monte Carlo simulation with given chemical potentials). We can, however, not guarantee that this statement holds generally.

Analytical results have been derived for a special set of parameters, from which we conclude a general trend towards an average site occupation of one half and a linear voltage – or chemical potential – drop along an electron transport chain with a linear topology. Introducing bottlenecks, we have to resort to numerical solutions that optimize the occupations and the local chemical populations simultaneously. Around these bottlenecks, both the chemical potentials and the electron populations show a drop. In retrospect, this justifies the setup of a considerable numerical apparatus that treats both quantities on an equal footing. We have verified and rationalized our methodology by relaxing the mean-field approach to one of these variables, the populations, via kinetic Monte Carlo simulations. We have demonstrated the application of the scheme to four-center proteins of biological relevance. Here, the comparatively simple current-chemical potential characteristics reflect a remarkably complex underlying pattern of chromophore populations. We are confident that the extension of our work to a catalogue of charge transfer proteins will lead to new perspectives and microscopic interpretations of biological electron transport.

The computational approach applied is robust and efficient, for the systems studied here the calculations can be performed within a few minutes on a desktop computer. Nevertheless, we find it desirable to have a more controllable numerical procedure available for the optimization of the local chemical potentials, e.g. by a nonlinear simplex method. An important physical aspect still missing are long-range intersite electron–electron interactions, which can also be obtained using dielectric theory computations.⁵⁸ We have not yet touched more complex topologies involving parallel conduction channel or dead-end storage sites, although they can be addressed by the methods described here. Furthermore, the high potentials applied in nanoscopic setup may have an impact on the protein structure and on the protonation pattern, which may in turn effect the characteristic parameters of Marcus' theory.

We note that coherence effects in biological charge transfer have been discussed recently,^59–61 which lie outside the range of Marcus' theory used here. The equations outlined above can be adapted to corrections of the Marcus scheme, provided they can be formulated in the context of simple rate equations.

Finally, we briefly comment on the biochemical and evolutionary aspects of our findings. Multicenter charge transfer chains can efficiently bridge membranes, or they can store electrons for many-electron redox reactions, as in the NrfHA, cytochrome and respiratory complex I protein complexes referenced above. In addition, the current–voltage characteristics are subject to the parameters of the chain and the number of cofactors. They are the subject of an evolutionary optimization process, which can adapt to requirements such as handling an increased oxygen intake on the mitochondrial level, which can increase by a factor of ten to 65 between states of rest and activity for mammals and birds.⁶²

Data availability statement

Programs for the computations involving metallic leads and single states as entry points are made available via the GitHub server (https://github.com/theochem-freiburg/newtp), including input and output files used in the paper as examples.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We gratefully acknowledge financial support by the Deutsche Forschungsgemeinschaft via the research group FOR 5099 ‘Reducing complexity of nonequilibrium systems’, projects P2 and P3. It is a pleasure to thank Thorsten Hugel, Thorsten Friedrich, Jetmir Haxhija, Tanja Schilling, Gerhard Stock, Steffen Wolf and Heinz-Peter Breuer for fruitful discussions.

References

1. B. E. Ramirez, B. G. Malmström, J. R. Winkler and H. B. Gray, The currents of life: the terminal electron-transfer complex of respiration, Proc. Natl. Acad. Sci. U. S. A., 1995, 92, 11949–11951.
2. H. A. Wagenknecht and H. B. Gray, Charge Transfer in DNA: From Mechanism to Application, Wiley-VCH, Weinheim, 2006.
3. K. Schulten, C. E. Swenberg and A. Weller, A Biomagnetic Sensory Mechanism Based on Magnetic Field Modulated Coherent Electron Spin Motion, Z. fur Phys. Chem., 1978, 111, 1–5.
4. U. Brandt, Energy Converting NADH: Quinone Oxidoreductase (Complex I), Annu. Rev. Biochem., 2006, 75, 69–92.
5. J. Hirst, Mitochondrial Complex I, Annu. Rev. Biochem., 2013, 82, 551–575.
6. T. Friedrich, On the mechanism of respiratory complex I, J. Bioenerg. Biomembr., 2014, 46, 255–268.
7. J. Deisenhofer, O. Epp, I. Sinning and H. Michel, Crystallographic refinement at 2.3 Ã… Resolution and Refined Model of the Photosynthetic Reaction Centre from Rhodopseudomonas viridis, J. Mol. Biol., 1995, 246, 429–457.
8. B. Dohse, P. Mathis, J. Wachtveitl, E. Laussermai and S. Iwata, et al., Electron Transfer from the Tetraheme Cytochrome to the Special Pair in the Rhodopseudomonas viridis Reaction Center: Effect of Mutations of Tyrosine L162, Biochemistry, 1995, 34, 11335–11343.
9. E. S. Medvedev, A. I. Kotelnikov, N. S. Goryachev, B. L. Psikha and J. M. Ortega, et al., Protein dynamics control of electron transfer in reaction centers from Rps. viridis, Mol. Simul., 2006, 32, 735–750.
10. M. L. Rodrigues, T. F. Oliveira, I. A. C. Pereira and M. Archer, X-ray structure of the membrane-bound cytochrome c quinol dehydrogenase NrfH reveals novel haem coordination, EMBO J., 2006, 25, 5951–5960.
11. J. Frank, Generalized single-particle cryo-EM – a historical perspective, Microscopy., 2015, 65, 3–8.
12. E. G. Yoga, H. Angerer, K. Parey and V. Zickermann, Respiratory complex I – Mechanistic insights and advances in structure determination, Biochim. Biophys. Acta, Bioenerg., 2020, 1861, 148153.
13. P. Mitchell, Coupling of Phosphorylation to Electron and Hydrogen Transfer by a Chemi-Osmotic type of Mechanism, Nature, 1961, 191, 144–148.
14. F. S. Saleh, M. R. Rahman, T. Okajima, L. Mao and T. Ohsaka, Determination of formal potential of NADH/NAD⁺ redox couple and catalytic oxidation of NADH using poly(phenosafranin)-modified carbon electrodes, Bioelectrochemistry, 2011, 80, 121–127.
15. M. T. Huynh, C. W. Anson, A. C. Cavell, S. S. Stahl and S. Hammes-Schiffer, Quinone 1e⁻ and 2e⁻/H⁺ Reduction Potentials: Identification and Analysis of Deviations from Systematic Scaling Relationships, J. Am. Chem. Soc., 2016, 138, 15903–15910.
16. Y. Guo, J. R. Stroka, B. Kandemir, C. E. Dickerson and K. L. Bren, Cobalt Metallopeptide Electrocatalyst for the Selective Reduction of Nitrite to Ammonium, J. Am. Chem. Soc., 2018, 140, 16888–16892.
17. J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment, World Scientific, Singapore, 2010.
18. D. Segal and A. Nitzan, Conduction in molecular junctions: inelastic effects, Chem. Phys., 2002, 281, 235–256.
19. M. Thoss and F. Evers, Perspective: Theory of quantum transport in molecular junctions, J. Chem. Phys., 2018, 148, 030901.
20. D. Gelbwaser-Klimovsky, A. Aspuru-Guzik, M. Thoss and U. Peskin, High-Voltage-Assisted Mechanical Stabilization of Single-Molecule Junctions, Nano Lett., 2018, 18, 4727–4733.
21. J. T. Lü, M. Brandbyge and P. Hedegard, Blowing the Fuse: Berrys Phase and Runaway Vibrations in Molecular Conductors, Nano Lett., 2010, 10, 1657–1663.
22. N. Bode, S. V. Kusminskiy, R. Egger and F. von Oppen, Scattering Theory of Current-Induced Forces in Mesoscopic Systems, Phys. Rev. Lett., 2011, 107, 036804.
23. C. Schinabeck, A. Erpenbeck, R. Härtle and M. Thoss, Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems, Phys. Rev. B, 2016, 94, 201407(R).
24. J. K. Sowa, J. A. Mol, G. A. D. Briggs and E. M. Gauger, Beyond Marcus theory and the Landauer-Büttiker approach in molecular junctions: a unified framework, J. Chem. Phys., 2018, 149, 154112.
25. J. K. Sowa, N. Lambert, T. Seideman and E. M. Gauger, Beyond Marcus theory and the Landauer-Büttiker approach in molecular junctions. II. A self-consistent Born approach, J. Chem. Phys., 2020, 152, 064103.
26. C. J. Murphy, M. R. Arkin, Y. Jenkins, N. D. Ghatlia, S. H. Bossmann and N. J. Turro, et al., Long-Range Photoinduced Electron Transfer Through a DNA Helix, Science, 1993, 262, 1025–1029.
27. B. Giese, Long-Distance Electron Transfer Through DNA, Annu. Rev. Biochem., 2002, 71, 51–70.
28. H. A. Wagenknecht, Reductive Electron Transfer and Transport of Excess Electrons in DNA, Angew. Chem., Int. Ed., 2003, 42, 2454–2460.
29. G. B. Schuster, Long-Range Charge Transfer in DNA I. Topics in Current Chemistry, Springer, Berlin/Heidelberg, 2004.
30. D. Porath, A. Bezryadin, S. de Vries and C. Dekker, Direct measurement of electrical transport through DNA molecules, Nature, 2000, 403, 635–638.
31. F. C. Grozema, Y. A. Berlin and L. D. A. Siebbeles, Mechanism of Charge Migration through DNA: Molecular Wire Behavior, Single-Step Tunneling or Hopping?, J. Am. Chem. Soc., 2000, 122, 10903–10909.
32. D. Porath, G. Cuniberti and R. D. Felice, Charge Transport in DNA-Based Devices, Springer, Berlin/Heidelberg, 2004, pp. 183–228.
33. T. Cramer, S. Krapf and T. Koslowski, DNA Charge Transfer in an External Field: An Atomistic Approach, J. Phys. Chem. C, 2007, 111, 8105–8110.
34. R. A. Marcus, On the Theory of Oxidation–Reduction Reactions Involving Electron Transfer. I, J. Chem. Phys., 1956, 24, 966–978.
35. R. A. Marcus and N. Sutin, Electron transfers in chemistry and biology, Biochim. Biophys. Acta, Rev. Bioenerg., 1985, 811, 265–322.
36. R. S. Farid, C. C. Moser and P. L. Dutton, Electron transfer in proteins, Curr. Opin. Struct. Biol., 1993, 3, 225–233.
37. H. B. Gray and J. R. Winkler, Electron transfer in proteins, Annu. Rev. Biochem., 1996, 65, 537–561.
38. D. Gnandt, S. Na and T. Koslowski, Simulating biological charge transfer: Continuum dielectric theory or molecular dynamics?, Biophys. Chem., 2018, 241, 1–7.
39. N. F. Mott, Metal-Insulator Transitions, Taylor & Francis, London, 2nd edn, 1990.
40. H. Böttger and V. V. Bryksin, Hopping conduction in solids, Akademie Verlag, Berlin, 1985.
41. C. E. D. Chidsey, Free energy and temperature dependence of electron transfer at the metal-electrolyte interface, Science, 2002, 251, 919–922.
42. S. Na, A. Bauss, M. Langenmaier and T. Koslowski, Thermodynamic integration network study of electron transfer: from proteins to aggregates, Phys. Chem. Chem. Phys., 2017, 19, 18938–18947.
43. J. Jeans, The mathematical theory of electricity and magnetism, chapter IX, Cambridge University Press, Cambridge, 1908.
44. T. Koslowski and A. Wilkening, A combined Kirchhoff-Master equation approach to electronic transport in one-dimensional molecular conductors, Chem. Phys., 2010, 369, 22–26.
45. J. J. Dongarra, C. B. Moler, J. R. Bunch and G. W. Stewart, LINPACK users' guide, SIAM, 1979.
46. W. M. Young and E. W. Elcock, Monte Carlo studies of vacancy migration in binary ordered alloys: I, Proc. Phys. Soc., 1966, 89, 735.
47. A. B. Bortz, M. H. Kalos and J. L. Lebowitz, A new algorithm for Monte Carlo simulation of Ising spin systems, J. Comput. Phys., 1975, 17, 10–18.
48. L. Li, C. Li, S. Sarkar, J. Zhang, S. Witham and Z. Zhang, et al., DelPhi: a comprehensive suite for DelPhi software and associated resources, BMC Biophys., 2012, 5, 9.
49. C. Li, Z. Jia, A. Chakravorty, S. Pahari, Y. Peng and S. Basu, et al., DelPhi Suite: New Developments and Review of Functionalities, J. Comput. Chem., 2019, 40, 2502–2508.
50. J. W. Ponder and D. A. Case, Force Fields for Protein Simulations, in Protein Simulations of Advances in Protein Chemistry, Academic Press, 2003, vol. 66, pp. 27–85.
51. D. Case, I. Ben-Shalom, S. Brozell, D. Cerutti, T. Cheatham and V. Cruzeiro, et al., AMBER 2018, San Francisco, 2018.
52. C. C. Moser, J. M. Keske, K. Warncke, R. S. Farid and P. L. Dutton, Nature of biological electron transfer, Nature, 1992, 355, 796–802.
53. P. Voigt and E. W. Knapp, Tuning Heme Redox Potentials in the Cytochrome c Subunit of Photosynthetic Reaction Centers, J. Biol. Chem., 2003, 278, 51993–52001.
54. T. Steinbrecher, T. Koslowski and D. A. Case, Direct Simulation of Electron Transfer Reactions in DNA Radical Cations, J. Phys. Chem. B, 2008, 112, 16935–16944.
55. F. Stein, D. Drung, L. Fricke, H. Scherer, F. Hohls and C. Leicht, et al., Validation of a quantized-current source with 0.2 ppm uncertainty, Appl. Phys. Lett., 2015, 107, 103501.
56. O. P. Hamill, E. Neher, B. Sakmann and F. J. Sigworth, Improved patch-clamp techniques for high-resolution current recordings from cells and cell-free membrane patches, Pflügers Arch., 1981, 391, 85–100.
57. C. Schmidt, M. Mayer and H. Vogel, A Chip-Based Biosensor for the Functional Analysis of Single Ion Channels, Angew. Chem., 2000, 112, 3267–3270.
58. D. Gnandt and T. Koslowski, Long-range electron-electron interaction and charge transfer in protein complexes: a numerical approach, Phys. Chem. Chem. Phys., 2019, 21, 18595–18604.
59. Y. Eshel, U. Peskin and N. Amdursky, Coherence-assisted electron diffusion across the multi-heme protein-based bacterial nanowire, Nanotechnology, 2020, 31, 314002.
60. Z. Futera, I. Ide, B. Kayser, K. Garg and X. Jiang, et al., Coherent Electron Transport across a 3 nm Bioelectronic Junction Made of Multi-Heme Proteins, J. Phys. Chem. Lett., 2020, 11, 9766–9774.
61. Y. Agam, R. Nandi, A. Kaushansky, U. Peskin and N. Amdursky, The porphyrin ring rather than the metal ion dictates long-range electron transport across proteins suggesting coherence-assisted mechanism, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 32260–32266.
62. C. M. Bishop, The maximum oxygen consumption and aerobic scope of birds and mammals: getting to the heart of the matter, Proc. Biol. Sci., 1999, 266, 2275–2281.

---