Protein–water electrostatics and principles of bioenergetics

Abstract

Despite its diversity, life universally relies on a simple basic mechanism of energy transfer in its energy chains—hopping electron transport between centers of electron localization on hydrated proteins and redox cofactors. Since many such hops connect the point of energy input with a catalytic site where energy is stored in chemical bonds, the question of energy losses in (nearly activationless) electron hops, i.e., energetic efficiency, becomes central for the understanding of the energetics of life. We show here that standard considerations based on rules of Gibbs thermodynamics are not sufficient, and the dynamics of the protein and the protein–water interface need to be involved. The rate of electronic transitions is primarily sensitive to the electrostatic potential at the center of electron localization. Numerical simulations show that the statistics of the electrostatic potential produced by hydration water are strongly non-Gaussian, with the breadth of the electrostatic noise far exceeding the expectations of the linear response. This phenomenon, which dramatically alters the energetic balance of a charge-transfer chain, is attributed to the formation of ferroelectric domains in the protein's hydration shell. These dynamically emerging and dissipating domains make the shell enveloping the protein highly polar, as gauged by the variance of the shell dipole which correlates with the variance of the protein dipole. The Stokes-shift dynamics of redox-active proteins are dominated by a slow component with the relaxation time of 100–500 ps. This slow relaxation mode is frozen on the time-scale of fast reactions, such as bacterial charge separation, resulting in a dramatically reduced reorganization free energy of fast electronic transitions. The electron transfer activation barrier becomes a function of the corresponding rate, self-consistently calculated from a non-ergodic version of the transition-state theory. The peculiar structure of the protein–water interface thus provides natural systems with two “non's”—non-Gaussian statistics and non-ergodic kinetics—to tune the efficiency of the redox energy transfer. Both act to reduce the amount of free energy released as heat in electronic transitions. These mechanisms are shown to increase the energetic efficiency of protein electron transfer by up to an order of magnitude compared to the “standard picture” based on canonical free energies and the linear response approximation. In other words, the protein–water tandem allows both the formation of a ferroelectric mesophase in the hydration shell and an efficient control of the energetics by manipulating the relaxation times.

---

David N. LeBard

David N. LeBard received his PhD in the theoretical chemistry group of D. V. Matyushov at Arizona State University (2008). He worked as a postdoctoral research fellow in the Center for Molecular Modeling under M. L. Klein at the University of Pennsylvania (2009). Currently he is still conducting postdoctoral research under M. L. Klein, now at Temple University, in the Institute for Computational Molecular Science. His research interests include protein electron transfer and solvation, large scale molecular modeling with emerging technologies, and ligand-gated ion channels.

Dmitry V. Matyushov

Dmitry V. Matyushov received his MS in Chemical Physics from the Moscow Institute for Physics and Technology and his PhD in Theoretical Physics from Kiev's State University in Ukraine. He was a Lise Meitner Postdoctoral Fellow of the Austrian Science Foundation and did postdoctoral work at Colorado State University and University of Utah. He currently holds a joint appointment in the Chemistry and Physics departments at Arizona State University. Research interests involve electron transport in molecular systems, theory of polar liquid and solvation, and properties of hydrated proteins.

---

1. Introduction

Life is based on condensed molecular systems. The main distinction between this type of materials and conventional solids is the presence of energetic and structural heterogeneity. Both arise from a higher energy concentrated in the network of chemical bonds within the molecules compared to much softer intermolecular potentials. The nonuniform distribution of structural stiffness has the effect of localizing thermal excitations on the molecular scale. This length scale is much shorter than the typical meso-to-macro length scale of excitations in solids (phonons, plasmons, etc.). Signals and information in molecular materials do not propagate over large distances, and the mechanisms of molecular transformations become localized and microscopic. This simple length scale limitation imposes quite significant restrictions on biological design principles typically involving cascades of localized interactions and motions.

While the localized nature of biomolecular interactions is well rooted in the prevailing thinking of biochemistry, the difficulty it immediately creates is how to reconcile the locality of chemistry with typical mesoscopic distances of molecular biology.¹ Given the locality of each step, these length scales require concerted sequences of often large numbers of microscopic events. The movement of redox and photonic energy in respiratory chains and photosynthesis follows the same central mechanism.² Each individual transition step of an electron-transport chain is limited to the distance of electron tunneling not exceeding ≃15 Å, but the total path between the points of energy entry and energy consumption can be much longer, through a large number of electron hops adding up to an overall directional electron current.³ The focus of this report is on this problem, with an eye towards both the fundamental mechanisms of energy transfer in biology's energy chains and the generic principles of their efficiency. Such a theoretical perspective is required not only in order to understand biological electron transport, but also to set the stage for the design of artificial photosynthesis.

Processes involved in biological energy transfer can generally be placed under the roof of a fundamental problem of energy transfer in systems with a broad heterogeneity of rigidity and dominant relaxation modes. The understanding of elastic and electrostatic properties of heterogeneous nono-scale materials is among the fundamental challenges that studies of hydrated proteins pose to us. Whether the development of artificial photosynthesis will entirely focus on homogeneous solid materials or will involve more heterogeneous systems with molecular/nano-crystalline centers dissolved in glassy matrices⁴ remains to be seen. Irrespective of that challenge, there is still a need in theoretical foundations of bioenergetics, which require studies of electronic energy transport between molecular cofactors dissolved in soft molecular solvents.

The “standard picture” of electron transfer between localized states communicates the idea of tunneling when the energy resonance is achieved between dynamically readjusting donor and acceptor electronic states.⁵ When the electron arrives to a cofactor site, the electronic level sinks due to the deformation of the surrounding medium (solvation). The magnitude of this stabilization is quantified by the reorganization (free) energy λ. The subsequent activationless tunneling to the next localization site, the acceptor level, then requires the acceptor to be lower in free energy by λ (Fig. 1). This amounts to imposing the linear response restriction (Marcus theory) that the reaction Gibbs energy ΔG is equal to −λ for an activationless electron hop.⁵

Fig. 1 A schematic representation of the energetics of a redox energy transfer chain. The energy input is provided by either the photochemical excitation or in the form of redox energy. The electron transfer chain connects the point of energy input with a catalytic site where the remaining redox energy is stored in the form of chemical bonds. The free energy λ ≃ 0.7 eV^3,6 is released to the phonon bath at each activationless transition. With the typical energy input of 1.0–1.5 eV, this scheme, if correct, represents an energetically wasteful operation design leaving no redox energy to the catalytic site.

The requirement of unidirectional electron flow in biology's energy chains is typically realized in terms of low barriers for forward reactions.² The electron transport can then be viewed as a downhill motion of the electron, in which it loses the free energy λ in each activationless hop (Fig. 1). The overall energetic efficiency of the electron-transfer chain will then critically depend on the magnitude of λ since the free energy −nλ will be lost to the thermal phonon bath in n activationless steps. This conventional view of a chain of energy dissipation typically results in the free energy loss significantly exceeding what is attained by biological systems. The resolution of this “energy puzzle of life” is a part of our discussion.

While the conceptual framework of Marcus theory can generally be applied to hopping conductivity in condensed media, two views of the role of proteins in promoting electron transport in biology have been put forward. In the first (physicist's) view, the protein matrix is considered a rigid static environment that is hard to polarize. This electrostaticaly inactive medium provides a low-dielectric “grease island” for the electron localization at a sequence of cofactors (hemes, chlorophylls, etc.). In the second (biologist's) view, the protein is viewed as a product of natural selection, highly optimized for its function.

The rigid protein matrix of the physicist’s representation of biological electron transport incorporates vibrationally stiff cofactors, additionally restricted in their motions by the protein. Both components of the medium are then little polarized (deformed) when the electron arrives at the center of localization. High rigidity of the protein/cofactor nuclear bath implies a low reorganization energy λ and hence low energy losses in a sequence of electron hops. This general view of electronic transitions in structurally rigid and electrostatically non-polar protein media is supplemented with the view of non-specific electronic tunneling mostly determined by an average density of chemical groups in the medium.³

The second, philosophically distinct, approach to biology's electron transport anticipates a picture in which both the electrostatic environment and the tunneling probability are optimized for best performance by natural selection.⁶ The extreme end of this line of thought—which is opposite to the picture of a structureless rigid continuum—views proteins as highly perfected machines with optimal energetics and preferential pathways⁷ for electron tunneling. A related notion, coming from viewing biomolecules as the primary target of natural selection, is that water serves merely as “spectator” in the theater of biological action providing solvation and diffusional delivery of the reactants to the points of localized enzymatic events. It is however becoming increasingly clear that the active volume of biological machines is not limited to their structural units and instead extends into their hydration shells,^8,9 an almost trivial statement given that most enzymatic activity takes place at the protein–water interface. A more nontrivial question is how far this action can propagate into the solvent and whether some kind of long-distance (cooperative) action can be assigned to such propagation. The concept of locality of biomolecular interactions is thus ripe for re-examination,^8–10 with the goal of establishing the actual length-scales of their action, particularly when applied to the crowded environment of the cell.

Neither of the above views of the role of proteins in electron transfer are free from contradictions. If proteins merely serve as a rigid non-polar environment, cofactors and active sites of enzymes should be buried well within the protein matrix to avoid solvation by highly polar water. This requirement would create a significant design problem for the inter-protein communications and diffusional delivery of the reactants. On the other hand, the notion of proteins as perfect machines optimized by natural selection does not answer the question of which mechanisms and properties make them so much more efficient than commonly used synthetic redox couples. It is a goal of this communication to summarize some fundamental design principles which provide proteins with advantages as electron-transfer media. We advocate here the view that in optimizing energetic efficiency natural selection has targeted the activation free energy (exponential factor of the rate) instead of tunneling pathways (rate preexponent).

The picture of the energetics of redox energy transfer in biology presented here departs in many ways from the traditional paradigm. We suggest that many functions of the redox energy transfer in biology do not rely on the rigid, non-polar interior of the proteins and instead make full use of the flexible protein–water interface. One might anticipate that the softness of the interface should make electron transport extremely inefficient due to strong solvation and electron localization at electron-transfer cofactors (large λ). This is in fact avoided due to a particular structure of water at proteins' surfaces allowing the combination of sluggish dynamics with strongly non-Gaussian statistics of the electrostatic potential. Both these factors severely reduce the amount of free energy released to heat in an activationless electronic transition and allow energetic efficiency sufficient for multiple electron hops within the electron transport chain.

What rises to prominence from this analysis is the interfacial water forming a mesoscale ensemble around a typical protein. In terms of electron transport properties, proteins do not seem to be unique machines of natural selection. Instead, they make use of the collective orientational motions of waters at the mesoscopic interface producing non-Gaussian electrostatic fluctuations at the active sites of redox reactions and cofactors of energy chains. This also implies that, rather than being buried inside the protein matrix, the cofactors and active sites are allowed to come forward to the protein's interface to take full advantage of the large breadth of electrostatic noise reducing activation barriers of electronic transitions.

2. Standard picture of molecular electron transport

As mentioned above, macroscopic molecular systems are energetically and structurally heterogeneous. The main consequence of that for the transport of energy is a significant extent of localization of both the electronic states and the nuclear excitations. The electronic states of most redox agents are well localized within the corresponding molecular moieties, and electronic transitions occur as a result of resonance hops between the states of localization. The probability of each hop is proportional to squared matrix element of the electronic overlap V and the Franck–Condon factor 〈δ(ΔE(Q))〉. The latter describes the probability of bringing the energy gap ΔE(Q) between the acceptor and donor electronic states to zero by thermal fluctuations (phonons and overdamped modes) of nuclear coordinates Q. The energy gap X = ΔE(Q) then naturally becomes the reaction coordinate for radiationless transitions,^11–14 with the transition point at X = 0.

The energetics of electronic transitions has been historically described following the route suggested by Franck–Condon factors studied by optical spectroscopy of molecules. The activation barrier is obtained from two free energy surfaces¹⁵F_i(X) crossing at the point of zero energy gap X = 0 when tunneling becomes allowed. Each free energy surface corresponds to a particular redox state enumerated here as i = 1 and 2 and is defined by the following equation

exp[−βF_i(X)] ∝ 〈δ(ΔE(Q) − X)〉_i, (1)

where β = 1/(k_BT) is the inverse temperature. Here, the angular brackets denote an ensemble average over the statistical distribution of the system characterized by the Hamiltonian H_i(Q). Since ΔE(Q) = H₂(Q) − H₁(Q), energy conservation requires a linear relation between the two free energy surfaces, F₂(X) = F₁(X) + X; the curvatures of F_i(X) have to coincide at each value of the reaction coordinate X.

The stochastic energy gap coordinate ΔE(Q) depends on many nuclear coordinates, and so it seems reasonable, from the central limit theorem, to assume that X obeys Gaussian statistics^16,17 with the average X_0i = 〈ΔE〉_i and variance σ²_i = 〈(δΔE)²〉_i in each electronic state. The free energy surfaces are then parabolas as functions of the reaction coordinate X

The linear relation between F_i(X) and the requirement that they intersect at the transition point X = 0, F₁(0) = F₂(0), put a number of restrictions on the parameters of the model, σ₁ = σ₂ = σ, ΔX = X₀₁ − X₀₂ = βσ² and ΔF = F₀₂ − F₀₁ = X̅, where the mean of the average energy gaps is X̅ = (X₀₁ + X₀₂)/2.

The distribution of the energy gaps in electron-transfer systems can be probed by optical spectroscopy that reports the optical Franck–Condon factor 〈δ(ΔE − hν)〉_i as a function of the photon frequency ν. The average energy gap X_0i = 〈ΔE〉_i then corresponds to the maximum of an optical transition band and ΔX = X₀₂ − X₀₁ is the Stokes shift. Half of the Stokes shift can be used to define the reorganization energy λ^St = ΔX/2, where the superscript “St” emphasizes that the first moments of the stochastic variable ΔE(Q) are used here. Alternatively, one can use the variance of the energy gap to determine the reorganization energy λ^var = βσ²/2. According to the linear relation between F_i(X),

λ^St = λ^var (3)

when the fluctuations of the energy gap ΔE(Q) are Gaussian.

The free energy of activation for the forward electronic transition can be found as ΔF^† = X²₀₁/4λ^var. When X₀₁ = 0, the forward transition is activationless, and the reaction occurs with the negative reaction free energy

The reorganization energy λ, which for Gaussian fluctuations of the energy gap is equal to both λ^St and λ^var, thus has an important physical meaning of the energy dissipated into the heat at each activationless electronic transition (Fig. 1). The minimization of λ must therefore be an essential component of photosynthetic designs aiming at both the unidirectional electron flow and the energetic efficiency.

3. Intermolecular interactions and Gaussian statistics

The notion of the Gaussian fluctuations of the energy gap between the donor and acceptor electronic states does not specify the origin of molecular interactions or nuclear modes contributing to the thermal noise. The interaction potentials for molecular systems are more or less established and can be broadly divided into Coulomb and non-polar interactions. Coulomb interactions are caused by inhomogeneous distribution of molecular charge often represented by either partial atomic charges or by multipolar expansions. Non-polar interactions are related to the ability of the electronic density to alter in an external electric field and are typically given in the form of dispersion and induction forces. The overall dependence of the donor–acceptor energy gap on nuclear coordinates Q can therefore be given by the following sum of the interaction terms

ΔE(Q) = ΔE₀ + ΔE^disp(Q) + ΔE^ind(Q) + ΔE^C(Q). (5)

Here, ΔE₀ is the vacuum energy gap between the cofactors and “disp”, “ind”, and “C” stand, correspondingly, for dispersion, induction, and Coulomb interactions with the protein and the surrounding liquid.

The distance dependencies and effective ranges of action of these interactions are quite different: r⁻⁶ for dispersion, r⁻⁴ for induction, and r⁻¹ for the Coulomb interaction energies. There are therefore two questions relevant to this significant distinction in the interaction range: (i) whether each interaction component obeys the Gaussian statistics and (ii) whether the relation between the first and second cumulants of the interaction energy, following from the fluctuation-dissipation theorem^18,19 (see below), applies to all types of interaction energies.²⁰

The Gaussian statistics of Coulomb interactions is often described by the linear response approximation. The basic idea is sketched in Fig. 2. If a small probe charge q is placed in a polar liquid, its electric field perturbs the orientations of the dipoles in the liquid to the extent proportional to the solute's field. If the charge is then removed and the dipoles of the liquid are allowed to fluctuate, their fluctuations will follow the same statistical rules as those imposed by an external perturbation. The average potential 〈ϕ〉 of the liquid at the position of the probe charge can then be connected to the fluctuations of the liquid in the absence of the external charge and one gets the classical version of the fluctuation-dissipation theorem,^18,19 −〈ϕ〉 = βq〈(δϕ)²〉 = βq〈(δϕ)²〉₀. The second bracket here, 〈…〉₀, denotes the average taken in the absence of the probe charge. The equality between two variances implies that the spectrum of the solvent fluctuations is not perturbed by the solute electric field in the linear-response approximation. The first equation establishes the equality between λ^St and λ^var [eqn (3)] and the second equation leads to the equality between the variances of the energy gap in electronic transitions, σ₁ = σ₂. The linear response is therefore equivalent to the Gaussian statistics.

Fig. 2 Cartoon representation of the linear response of a dipolar particle to the external field E of a negatively charged solute. The average dipole 〈m〉 is oriented along the field and its magnitude is proportional to the field. The dipolar response function χ is the proportionality coefficient. Thermal fluctuations of the same dipole in the absence of the external field (solute without charge) carry information about the response function χ (fluctuation-dissipation theorem).

The electronic states at the cofactors interact with a heterogeneous protein–water medium characterized by vastly different viscoelastic and dielectric properties. The two components, water (s) and protein (p), contribute additively to the first moment of ΔE(Q) in eqn (5). One then finds

λ^St = λ^St_s + λ^St_p. (6)

The variance reorganization energy contains, apart from the direct contributions from the water (λ^var_s) and protein (λ^var_p) components, the cross term λ^var_sp from cross-correlations between the water and protein interaction energies

λ^var = λ^var_s + λ^var_p + λ^var_sp. (7)

The corresponding linear response relations then read^21,22λ^St_s,p = λ^var_s,p + (1/2)λ^var_sp.

The rules of the linear response approximation, connecting the first and second cumulants of the Coulomb interaction potential, do not apply to short-ranged dispersion and induction interactions.²³ The first cumulant of ΔE^disp/ind, the average solute–solvent interaction energy, can, to a first approximation, be calculated by averaging it over the liquid structure unperturbed by the interaction.²⁴ The second cumulants are typically small and can often be neglected, resulting in 〈ΔE^disp/ind〉₁ ≃ 〈ΔE^disp/ind〉₂. Even though the contribution of the nonpolar interactions to the average energy gap can be significant,²⁵ their effect is mostly reduced to a uniform shift of both X₀₁ and X₀₂ without affecting either λ^St or λ^var. A small, but not negligible, contribution of the induction interaction to λ^var can, however, be identified for charge separation in the bacterial reaction center.^25,26

The Coulomb interactions thus mostly determine the values of the reorganization energies and, consequently, the free energy released as heat by each electron hop [eqn (4)]. The rules of linear response apply to the statistics of the electrostatic potential produced by the solvent within the solute if the perturbation of the solvent structure by the solute is small and the spectrum of the solvent fluctuations is not altered. This might not always be the case, in particular for hydration. When the solute establishes strong hydrogen bonds with the water solvent, or its size is sufficiently large to perturb the network of water's hydrogen bonds,²⁷ one can expect deviations from the Gaussian statistics of the electrostatic interactions.

Large solutes are indeed known to significantly affect the interfacial water structure.^9,27 Unsatisfied hydrogen bonds at the water-solute interface cause an energetic push for in-plane orientations of the water dipoles,²⁸ altering the average interfacial polarization. It was indeed found that λ^var deviates upward from λ^St for small solutes establishing strong hydrogen bonds with water.²⁹ Qualitatively similar, but quantitatively more dramatic upward deviation was found in a number of recent numerical simulations of hydrated proteins.^{21,22,26,30–32} These studies have reported gigantic reorganization energies from the energy gap variance, much exceeding the values from the Stokes shift

λ^var ≫ λ^St. (8)

This inequality breaks one of the basic relations of the linear response approximation and therefore carries quite significant consequences for the entire energetic balance of electron transport chains, as we discuss below. In a nutshell, this inequality implies that −ΔG = λ ≪ λ^St, λ^var in eqn (4) and less free energy is lost to heat.

4. Reorganization energies and Stokes shift dynamics

The electrostatic properties of the protein–water interface discussed here have been extracted from extensive molecular dynamics (MD) simulations of several hydrated proteins. The main body of the data came from analyzing the dynamics and thermodynamics of a half redox reaction for depositing an electron to the active site of metalloprotein plastocyanin from spinach.^21,22,32 In addition, MD simulations of the hydrated bacterial reaction center of Rb. Sphaeroides bacterium were used to analyze the sequence of photoinduced primary charge-separation transitions.^26,32 For all redox-active proteins, anomalies in the statistics of the electrostatic potential were caused by water. We therefore first discuss the water component of the electrostatic response at varying temperature summarized in Fig. 3.

Fig. 3 Water component S_s(t) of the Stokes shift correlation function at different temperatures indicated in the plot (a) and the reorganization energies λ^var_s (circles) and λ^St_s (diamonds) vs. temperature (b). The inset in (a) gives the relaxation times of the exponential tail of S_s(t); the inset in (b) shows the Binder parameter,³³δ_B = 1 − 〈(δΔE^C_s)⁴〉/[3〈(δΔE^C_s)²〉²], calculated from the water component ΔE^C_s of the Coulomb energy gap. T* and T_tr indicate, correspondingly, the onset of collective fluctuations and the structural transition of the hydration shell. The dashed horizontal line indicates the typical range of temperatures, T_D ≃ 200–240 K, assigned to a dynamical transition in proteins.^34,35 The results were produced from MD simulations of the metalloprotein plastocyanin in its oxidized redox state.

Fig. 3a shows the normalized correlation function S_s(t) of the Coulomb energy gap ΔE^C_s(Q(t)) arising from hydration water (water's component of the Stokes-shift correlation function). There is a clear growth of the slow exponential relaxation with increasing temperature (Fig. 3a). The exponential relaxation time (inset in Fig. 3a) passes through a spike at T_tr ≃ 220 K, reminiscent of critical slowing down at a point of phase transition. Further, the Binder parameter,³³ devised to distinguish between the first- and second-order phase transitions, also shows a downward spike at the same temperature pointing to a weak first-order phase transition²¹ (inset in Fig. 3b).

If a structural transition occurs in the sub-ensemble of the hydration waters, one expects variances of observables coupled to the corresponding order parameter to pass through upward spikes at the transition temperature.³³ This is exactly what is seen for the reorganization energy λ^var_s obtained from the variance of ΔE^C_s(Q(t)): it passes through a spike at T_tr, reaching the value of ≃8 eV at the maximum. This set of consistent observations points to some structural transition happening at T_tr and affecting the electrostatics of the hydration shells closest to the protein's active site. The exact nature of this transition is not clear although the non-Gaussian electrostatics of the protein–water interface at high-temperatures is clearly linked to polarized (ferroelectric) domains in the protein's hydration shell.²²

The transition temperature T_tr ≃ 220 K falls in the range of temperatures commonly associated with the dynamical transition in proteins, T_D ≃ 200–240 K.^34–36 This latter cross-over marks the change in the temperature slope of several observables, most notably atomic mean-square displacements^34,37,38 and dielectric absorption.³⁹ Several explanations of the observations were put forward. Doster and co-workers^35,37 proposed that a kinetic arrest of the primary α-relaxation of the medium, i.e. a true dynamic freezing of molecular translations described by mode-coupling theories,⁴⁰ is behind the cross-over. An alternative, more trivial, explanation is in terms of a local β-relaxation mode of the hydration shell entering the instrumental observation window at T_D.^36,38,41,42

The second, instrumental resolution, scenario proposed to explain laboratory dynamical transition^36,38,41,42 is clearly consistent with our data. The growth of the reorganization energy λ^var_s starts when a secondary relaxation mode of the hydration layer enters the observation window set up by the length of the simulation trajectory.³² This cross-over, however, occurs at a lower temperature of T* ≃ 200 K, before the structural transition at T_tr (Fig. 3b). The temperature T* is slightly above the glass transition temperature T_g ≃ 180 ± 15 K assigned to partially hydrated proteins,⁴² but should not be confused with it since we are dealing with an observation window ∼ 10⁻⁸ s, much shorter than 10² s used to define the glass transition, and, in addition, λ^var_s couples to a secondary relaxation process.

The appearance of a collective solvent mode in the observation window restores ergodicity of the corresponding fluctuations lost at lower temperatures where the linear response approximation, λ^var_s ≃ λ^St_s, is reasonably accurate. It means that the high-frequency wing of the Stokes-shift fluctuations, still unfrozen at low temperatures, actually follows the Gaussian statistics. On the contrary, the slow collective fluctuations steadily increase the breadth of the electrostatic noise above T*, resulting in a dramatic rise of λ^var_s(T).

The overall temperature dependence of the reorganization energy (except for the narrow spike at T_tr) follows from the rate of restoring ergodicity due to shortening, with increasing temperature, of the principal relaxation time of the mode responsible for the growth of λ^var_s(T). Slow fluctuations are not observable if their relaxation times are outside the observation window τ_obs. This “instrumental” cutoff can be incorporated as a step-wise frequency filter in the frequency spectrum of the water Stokes shift correlation function S_s(ω).^43,44 The breadth of the fluctuations at a given observation time is then given by integrating S_s(ω) over the frequencies above the lowest cutoff k = τ⁻¹_obs

λ^var_s(k) = 2λ^var_s∫^∞_kS_s(ω)dω. (9)

The non-ergodic reorganization energy λ^var_s(k) defined by this equation reduces to its equilibrium value λ^var_s in the limit τ_obs → ∞ of the canonical ensemble.

Eqn (9) is a general result applicable to electronic transitions in media with distributed relaxation times.²⁰ One still needs to stress that the effectiveness of altering the energetics (e.g., the activation barrier) by cutting the low-frequency part of the fluctuation spectrum critically depends on the portion of the overall free energy residing in that part of the spectrum. It is this latter property, that is a large weight of slow-relaxation motions in the overall Stokes shift, that makes proteins effective users of dynamics for the goal of free energy tuning.

The relative weights of the fast and slow relaxation components are remarkably different between small hydrated dyes and hydrated proteins (Fig. 4). While the fast ballistic component of S_s(t) takes ≃80% for the former,⁴⁵ it is only ≃20–30% for the latter. The relaxation time of the exponential tail is also very slow, ≃ 150 ps for primary bacterial charge separation and ≃650 ps for the reduction of plastocyanin. It is this slow relaxation tail that freezes in at T* making the ballistic, Gaussian modes the only contributor to the reaction activation at low temperatures. The more elastically rigid protein has a larger ballistic relaxation component in the protein correlation function S_p(t) (Fig. 4). It is therefore the water part of the system response where most anomalies are found.

Fig. 4 Normalized Stokes-shift correlation function of the water (S) and protein (P) components of the energy gap ΔE(Q) [eqn (5)]. Shown are the correlation functions of the oxidized state of plastocyanin (PC)²² and of the photoexcited state of the Rb. Sphaeroides bacterial reaction center (RC) dissolved in water using a detergent micelle.^25,26

5. Protein–water interface

Electrostatics of most polar liquids are reasonably reproduced by the statistics and dynamics of their dipolar polarization.⁴⁶ Correspondingly, the problem of the electrostatics of the protein–water interface can be addressed by looking at the structure and dynamics of the dipolar polarization field in the interfacial region. How large this interfacial region can be is the first question that needs to be addressed.

The experience learned from hydration of small molecules indicates that the water structure remains mostly unperturbed beyond the first two hydration layers.⁴⁷ The situation might be quite different for proteins. Neutron scattering,⁴⁸ depolarized light scattering,⁴⁹ and recent dielectric data in the terahertz frequency domain^50,51 have all suggested that water structure around proteins is perturbed to a distance of approximately ≃10–20 Å into the bulk.

Our simulation results generally support this conclusion, emphasizing however that this assessment applies only to second cumulants of the observables studied here. In other words, the variance of a property coupled to the interfacial polarization is sensitive to the collective modes of the interface, while the average of the same property mostly misses them. This important distinction needs to be recognized when comparing the results of different experimental techniques. We also emphasize that this distinction does not arise within the realm of linear response where the time autocorrelation function and the relaxation of the average carry the same dynamical information¹⁹ (cf. to the discussion of Fig. 2).

In compliance with this general rule, the Stokes-shift reorganization energy λ^St_s, which reflects the average electrostatic field of the medium, follows the standard framework of conventional solvation theories. This can be demonstrated by looking at how the value of λ^St_s is accumulated within a growing water shell of thickness r around the protein. The molecules of water are included into the r-shell if the distance from the oxygen to the nearest protein atom is smaller then r (Fig. 5). The function λ^St_s(r) saturates to λ^St_s calculated from the entire simulation box within the first two solvation layers around the protein (Fig. 6). On the contrary, the reorganization energy λ^var_s(r) develops a long range of convergence at T > T*, extending nearly five solvation layers into the bulk (Fig. 6). It drops to the level of λ^St_s(r) only at low temperatures T < T* where it also loses its long-ranged character. The development of a wide breadth of electrostatic fluctuations is therefore linked to a long length-scale of the corresponding nuclear mode.

Fig. 5 Cartoon of the hydration layer of thickness r around a protein's surface taken from a snapshot of the oxidized plastocyanin simulation in the NVT ensemble. The protein is rendered in cartoon representation overlayed with a transparent van der Waals surface. Water molecules, rendered as van der Waals spheres, are assigned to the r-shell if their oxygen atoms are within the distance r from the nearest protein atom. In this figure, a thickness of 6 Å was used, while foreground shell waters within 6 Å were removed to provide a more clear image of the protein.

Fig. 6 Reorganization energies λ^var_s(r) and λ^St_s(r) calculated for Red and Ox states of wild type plastocyanin by including waters from a shell of thickness r from the protein surface (Fig. 5). The calculated values are divided by their corresponding values at r = 20 Å. The vertical dashed line indicates the thickness of two solvation shells around the protein within which λ^St_s(r) nearly reaches its bulk value. The data are obtained from NPT simulations with N_s = 5886 waters in the simulation box.²²

A large value of λ^var is caused by extensive fluctuations of the dipole moment of the hydration layers closest to the protein's active site. Fig. 7 shows the distribution of the first-shell dipole moment of two redox-active proteins, plastocyanin and bacterial reaction center, and two redox-inactive proteins, lysozyme and ubiquitin, which we have studied as benchmarks.²² The distribution of the first-shell dipole moment around redox-active proteins is not only much broader, but is also anisotropic, significantly deviating from the isotropic Maxwellian distribution. This is a result of breaking the hydration layers into ferroelectric domains, probably also responsible for the spike in λ^var at T_tr (Fig. 3).

Fig. 7 Distribution of the dipole moments M^I_s of the first-shell waters and the dipole moments M_p (inset) of the proteins listed in Table 1.

The large-scale dipolar fluctuations of the hydration water seem to correlate with fluctuations of the protein dipole (Fig. 7 and Table 1). Elastic deformations of the hydration shell might therefore be linked to the ability of the protein's low-frequency vibrations to modulate the protein dipole moment. This information is encoded into the combination of the protein fold and the distribution of surface charge caused by ionizable surface residues. One might therefore speculate that redox-active proteins fold to create the conditions for large-scale thermal fluctuations of their charge distributions. This property is not needed for redox-inactive proteins and the corresponding dipole moment variance is much lower, as is the case of lysozyme and ubiquitin (Table 1).

Table 1 First and second moments of the dipole moment magnitudes of proteins and their first hydration shells

Protein^a	〈Mp〉	〈(δMp)^2〉	〈M^Is〉	〈(δM^Is)^2〉	κ ^Ib	R eff
a Average dipole moments in D and variances in kD^2. b Compressibility of the first hydration shell calculated as κ = 〈(δN^I)^2〉/〈N^I〉, where N^I is the number of first-shell waters. c Effective radius of the protein (Å) calculated by rolling the probe sphere of radius 2 Å on the van der Waals surface of the protein. d 25 ns NVE simulations^22 with Ns = 27 918 (ubiquitin) and Ns = 27 673 (lysozyme) waters in the simulation box. e Wild type protein in the oxidized state. 25 ns NPT simulations^22 with Ns = 21 076 waters in the simulation box. f 10 ns NPT simulations of the reaction center of Rb. Sphaeroides in the charge-separated state with the electron transferred from the special pair to the M-side bacteriochlorophyll.^26 The simulation cell contains Ns = 10506 waters.
Ubiquitin^d	256	0.91	48	2.5	0.22	17.2
Lysozyme^d	150	1.44	77	1.4	0.23	20.0
Plastocyanin^e	249	24.0	582	438	0.20	18.2
Reaction center^f	996	56.9	705	609	1.19	50

---

The qualitative distinction between the polarity of the hydration layers around redox-active and redox-inactive proteins is illustrated by the distance profiles of their dielectric constants ε(r) (Fig. 8). This function accounts for the dielectric response of the hydration layer of thickness r and is calculated from the variance of the water dipole moment M_s(r) within the r-shell (Fig. 5). The standard equation valid for the tin-foil implementation of Ewald sums in the simulation protocol is used: ε(r) = 1 + 4π〈[δM_s(r)]²〉/[3V(r)], where V(r) is the volume of the water shell. The absolute value of the dielectric constant depends on the shape of the dielectric sample⁴⁶ requiring corresponding shape corrections. Given the complex shape of r-shells replicating the protein's van der Waals surface, absolute values of ε(r) are not available to us and only the change Δε(r) relative to the value calculated at r = 20 Å is given in Fig. 8. The main result of this calculation is the qualitative difference in the polarity profiles between redox-active and redox-inactive proteins. While polarity is depleted near the surface of lysozyme, it is enhanced for plastocyanin. This enhanced surface polarization also propagates further into the bulk compared to the dielectric depletion layer of two other proteins.

Fig. 8 Change in the dielectric constant of the water shell of thickness r (Fig. 5) relative to the corresponding values at r = 20 Å for plastocyanin (PC), ubiquitin (Ubiq), and lysozyme (Lys).

We would like to emphasize that the anomalies in the structure and dynamics of the hydration layers discussed here apply to dipolar orientations only. We have not observed any unusual patterns in the density fluctuations of the interfacial water. The compressibility of the first hydration shell, somewhat below that of bulk water,²² is indeed consistent for all globular proteins (Table 1), with the exception of a higher compressibility of water around the reaction center. This latter result might be the consequence of a severe disruption of the water structure by the detergent micelle enveloping the membrane-exposed, hydrophobic belt of the reaction center complex.

6. Energetics of bacterial photosynthesis

The reaction center of bacterial photosynthesis is a membrane-bound protein complex which is dissolved in water both experimentally and in simulations by surrounding its membrane-bound hydrophobic belt with a detergent micelle. MD simulations^25,26 of this complex on the time-scale of 10–20 ns have revealed large magnitudes of λ^var far exceeding λ^St [eqn (8)], in qualitative agreement with the results for redox metalloproteins.^30–32 In contrast to half redox reactions studied for these latter systems, electronic transitions in the reaction center occur between centers of electron localization at bacteriochlorophyll and bacteriopheophytin cofactors (Fig. 9). The two electron hops within the primary charge separation event (from P to B_L and from B_L to H_L)⁵² occur on a short time-scale of about 3 ps (Rb. Sphaeroides bacterium). The rate constant of primary charge separation k_CS then sets up the observation frequency k_CS = τ⁻¹_obs that severely cuts into the frequency spectrum of the corresponding Stokes-shift dynamics. The non-ergodic reorganization energy obtained by using the entire (water and protein) Stokes-shift correlation function in eqn (9) is significantly reduced from λ^var = 2.36 eV from a 10 ns simulation trajectory to λ^var(k_CS) = 0.36 eV on the time-scale of the reaction (Fig. 10).

Fig. 9 Schematic arrangement of cofactors in the bacterial reaction center. P is the special pair, and B and H are monomeric bacteriochlorophylls and bacteriopheophytins, respectively. Electron transfer in wild-type reaction centers occurs almost exclusively along the L-branch of cofactors (subscript “L”), while the M-branch (subscript “M”) is mostly inactive.

Fig. 10 Reorganization energy λ^var(k) [eqn (7)] as a function of the reaction rate k = 1/τ_obs at 300 K obtained by using the entire (protein and water) Stokes-shift correlation function in eqn (9) (upper panel). The inset shows the experimental charge-separation rates⁵³vs. temperature. The lower panel shows the experimental rates of charge separation obtained at different temperatures⁵³ (circles) and hydrostatic pressures⁵⁴ (diamonds) against the average energy gap X₀₁(T,P). Both the rate constant k^exp_CS(T,P) and the average energy gap X₀₁(T,P) are normalized to their values at T = 300 K and P = 1 bar. The temperature dependence of X₀₁(T) at 1 bar was obtained from MD simulations. Experimental isothermal compressibility⁵⁴χ_T = 15 Mbar⁻¹ was used to find the pressure dependence of X₀₁(P) at 300 K. The theoretical temperature dependence of the rate (solid line) was produced by using the temperature slope of the induction shift 〈ΔE^ind〉 and the Stokes shift dynamics S(t) from MD simulations. The theoretical pressure dependence of the rate (dashed line) was produced by changing the induction shift with the experimental compressibility and assuming pressure-invariant Stokes shift dynamics. More details of the calculations are given in the ESI.†

The dependence of the Franck–Condon factors on the reaction rate poses a significant conceptual dilemma as the canonical ensemble and the relevant equilibrium thermodynamic functions cannot be used to determine the activation barrier (see also the discussion in section 7). The standard practice of condensed-phase formulations of the transition-state theory based on the canonical free energy barrier needs to be replaced with an alternative theory.

One formulation of such an algorithm is to use a frequency-dependent linear response function with the step-wise frequency cutoff,⁴⁴ as in eqn (9). This approach sets up the activation barrier depending on the reaction rate that needs to be found from solving a self-consistent kinetic equation. For electronic transitions following the recipe of Fermi's golden rule one obtains the electron-transfer rate k_ET by self-consistently solving the equation

k_ET ∝ V² exp[−X₀₁(k_ET)²/(2σ(k_ET)²)]. (10)

This simple prescription for Golden-rule quantum transitions allows generalization to more complex non-exponential kinetics driven by the coupling of the reaction coordinate X to overdamped (diffusional) modes of the thermal bath. The activation process in this case is a diffusional motion along the non-ergodic free energy surface F(k,X) which itself depends on the rate constant through the non-ergodic cutoff of the frequency spectrum. The free energy surface is then globally distorted compared to the canonical free energy F(X) (Fig. 11a).

Fig. 11 Panel (a): Cartoon of the free energy surfaces along the electron-transfer reaction coordinate X obtained in the canonical ensemble, F(X), and in the restricted canonical ensemble,⁴⁴F(k,X), when slow nuclear modes are frozen in. The shaded area at the top of the barrier of F(X) shows the region where dynamical solvent effect,^12,55 amounting to the diffusional drag on passing the activation barrier, is active. The activation barrier of the reaction, ΔF(k), in the non-ergodic formulation depends on the reaction rate k which needs to be calculated by self-consistently solving eqn (10) for the exponential Golden-rule kinetics and diffusional eqn (11) and (12) for the non-exponential kinetics in terms of P(X,t). Panel (b): Self-consistent solutions of eqn (11) and (12) for a number of reaction center mutants studied in ref. 56. The solid lines are the the results of theoretical calculations²⁵ and the points are experimental measurements. The initial slope of the population decay defines the rate constant k_ET used to calculate the free energy surfaces F(k_ET,X) in eqn (12).

The dynamics of the population P(X,t) of the initial reactants' state is given by Fokker–Planck dynamics along F(k, X)

∂P(X,t)/∂t = [L(k_ET,X) − k(X)]P(X,t). (11)

In this equation, L(k_ET,X) is a diffusional operatorwith the diffusion constant D(k_ET) and k(X) representing the Golden-rule rate of multi-phonon transitions induced by localized molecular vibrations at each value of the classical reaction coordinate.⁵⁵Eqn (11) and (12) are solved iteratively by defining the rate k_ET from the initial slope of the reactants' population P(t) = ∫⁰_−∞P(X,t)dX (Fig. 11b). The result is a self-consistently determined and non-exponential population P(t) which agrees quite well with the laboratory measurements on a number of reaction center mutants reported in ref. 56 (Fig. 11b).

One needs to emphasize that the dynamical arrest affecting the free energy profile is clearly distinct from the dynamical solvent, or Kramers friction,⁵⁷ effect entering the preexponent of the reaction rate.^12,55,58 The solvent effect of the Kramers kinetics is a consequence of diffusional friction restricting the number of successful barrier crossings at the barrier top (shown by the shadow region in Fig. 11a). In contrast, the dynamical non-ergodic arrest produces a global alteration of the free energy surface over which diffusion occurs (F(k,X) in Fig. 11a). The dynamical (Kramers) solvent effect is a part of the Fokker–Plank description given by eqn (11) and (12), and our formalism is reduced to the latter⁵⁵ when non-ergodic corrections can be neglected and F(k,X) becomes F(X).

Both the polar (Coulomb) and the non-polar (induction-dispersion) interaction potentials enter the instantaneous energy gap between the donor and acceptor electronic states [eqn (5)]. One then wonders what are the relative contributions of these potentials to different components of the activation free energy. The induction interaction typically makes a significant contribution to the average energy gap X₀₁, and is actually a dominant part of it in the case of bacterial charge separation. Physically, 〈ΔE^ind〉 is the difference of free energies of polarizing the protein–water medium by the electric field of the electron localized at the acceptor and donor electronic states.²³ If average density ρ and polarizability α are assigned to the molecular groups of the medium, 〈ΔE^ind〉 for charge separation can be estimated from a simple equation^25,43

〈ΔE^ind〉 = −e²4πραg, (13)

where e is the elementary charge and g = (2R_D)⁻¹ + (2R_A)⁻¹ − R⁻¹ is a geometrical factor, also appearing in the Marcus equation for the polar reorganization energy.⁵ Here, R_D and R_A are the effective radii of the donor and acceptor, respectively, separated by the distance R. The product ρα in eqn (13) can be expressed in terms of the refractive index of the protein. With the typical value of the protein refractive index n_prot = 1.473 and the geometric parameters R_D = R_A = 5.6 Å, R = 11.3 Å for the bacteriochlorophyll cofactors, eqn (13) yields 〈ΔE^ind〉 = −1.09 eV, in a good agreement with −1.03 eV from MD simulations.²⁶ This latter value far exceeds the contribution to the average energy gap X₀₁ from the Coulomb interactions, 〈ΔE^C〉 = −0.24 eV, making induction interactions the main solvent component in X₀₁.

The Coulomb interaction is the main contributor to the equilibrium reorganization energy λ^var, most of which is frozen on the picosecond time-scale of primary charge separation. The remaining 0.36 eV of non-ergodic reorganization energy from fast ballistic dynamics (upper panel in Fig. 10) splits into 0.12 eV from fluctuations of the induction gap ΔE^ind(Q) and 0.24 eV from fluctuations of the Coulomb gap ΔE^C(Q) [eqn (5)]. Under the conditions of a severe non-ergodic cutoff of the Coulomb reorganization, fluctuations of the induction interaction become a significant part of the reorganization energy.

Several electron-transfer reactions in the photosynthetic sequence are known to accelerate with lowering temperature.^53,59 This unusual temperature dependence (inset in the upper panel and circles in the lower panel in Fig. 10) has attracted much attention in the field. The traditional explanation requiring near-zero activation barrier,⁶⁰ and therefore X₀₁ ≃ 0, is not supported by recent laboratory⁵⁶ and numerical²⁵ evidence showing X₀₁ ≃ 0.17 eV. It turns out that the temperature variation of X₀₁(T) arising from protein's contraction is sufficient to explain the observed kinetic data.²⁵ Contraction increases the density of polarizable groups ρ in eqn (13), making the induction shift more negative on cooling. Since induction interactions dominate X₀₁, this component of electron-transfer energetics becomes critical for explaining the effect of external conditions on the rate constant.

In order to demonstrate the dominant effect of X₀₁(T) on k_CS(T) we have plotted in Fig. 10 the experimental rates from Fleming et al.⁵³ (circles) against X₀₁(T) from MD simulations (see ESI† for more details). The theoretical solid line is calculated by iteratively solving eqn (10) with λ^var(k_CS) obtained from eqn (9). The agreement between the calculated and experimental rates is good in the range of temperatures where quantum effects can be neglected.

Further, if contraction of the protein matrix mostly contributes to the variation of the rate in the temperature scans, an increase of the hydrostatic pressure must produce a comparable acceleration of the primary rate. Exactly this trend was obtained by Timpmann et al.⁵⁴ (diamonds in the lower panel in Fig. 10). When the experimental compressibility of the reaction center is used in eqn (13) to calculate the pressure variation of X₀₁(P), the solution of eqn (10) (dashed line in Fig. 10) again gives a trend consistent with observations. Either cooling or compression in the range of conditions studied so far have reduced the average donor–acceptor energy gap to about half of its physiological value, yielding comparable accelerations of the rate.

The approximate two-fold symmetry of the reaction center (Fig. 9) suggests that two branches of cofactors might be used for the electron transport. This does not happen, and only the right-hand (L) branch is significantly active. The most probable scenario is the asymmetry in the anchoring of the bacteriochlorophyll cofactors by protein's aminoacids, which lifts the energy gap ΔE₀ of charge separation to B_M by ≃ 0.4 eV²⁶ compared to ΔE₀ of charge separation to B_L. The energy shift slows the reaction down to ≃200 ps, sufficient to unleash collective nuclear solvation frozen for the fast transition along the L-branch. This change sets in yet another reaction mechanism absent for other electron-transfer reactions in the reaction center.

The special pair of the reaction center is a sandwich of closely spaced and electronically coupled bacteriochlorophyll cofactors P_M − P_L, where P_M and P_L are, correspondingly, the M and L subunits of the special pair. The electronic state of the photoexcited special pair, P*, is a superposition of the covalent, (P_M − P_L)*, and ionized, P⁻_M − P⁺_L, states. A small energy gap between these two states makes P* highly polarizable, with the largest principal component of the polarizability tensor along the dipole moment of the ionized state. This asymmetry of the molecular polarizability translates itself into an order-of-magnitude asymmetry in the reorganization energy λ^var between the L and M branches (Table 2).

Table 2 Energetic parameters (eV) of primary charge separation (P* →B) and charge recombination (B → P*) along L and M branches of the reaction center of Rb. Sphaeroides

Property	P* → BL	BL → P*	P* → BM	BM → P*
a Solvent component of the average energy gap, Xs = 〈ΔE^C〉 + 〈ΔE^ind〉; ΔE0 = 1.43(L), 1.87(M) eV. b Equilibrium reorganization energy is shown in parentheses.
X s	−1.26	−2.80	−1.19	−2.26
λ ^St	0.77		0.54
λ ^var	0.36(2.36)^b	1.69	4.98	2.92

---

A simple picture of what happens is as follows. The overall solvent reorganization energy is proportional to the average square of the sum of the permanent, m_CT, and induced, p, electron-transfer dipoles,⁶¹λ^var ∝ 〈(m_CT+p)²〉, where m_CT is the dipole moment created by moving the electron from the donor to the acceptor site. The prevalent direction of the fluctuating dipole p relative to m_CT will determine the value of the reorganization energy. Since the average induced dipole points in the direction of moving the electron along the M branch, the corresponding reorganization energy significantly exceeds that along the L branch.

The free energy surfaces of charge separation along the L and M branches are summarized in Fig. 12. Electron transfer to B_M is uphill, with an activationless backward reaction. This state is therefore never realized as a reaction intermediate, slowing down electron transport along the M branch. The physiological importance of this branch, and the sandwich structure of the special pair, are not clear; a strong electronic coupling between two subunits of the special pair is not preserved for Photosystem II.¹ However, this example points to yet another mechanism of regulating the pathway of electronic transitions by tuning the relative orientations of the electron-transfer dipole and anisotropic molecular polarizability. Once again, as in the case of polarization fluctuations at the protein–water interface, the statistics of the energy gap coordinate X become highly non-Gaussian,⁶¹ as is reflected in the non-parabolic free energy surface F₁(X) in Fig. 12.

Fig. 12 Free energy surfaces of the primary charge separation P* → B in the bacterial reaction center along the L and M branches (Fig. 9). The surfaces along the L branch were calculated on the 4 ps observation window, while the entire simulation trajectory of 10 ns was used to calculate the free energy surfaces along the M branch. The asymmetric shape of F₁(X) (marked as P* — B_M) for the M-branch charge separation is caused by a very large and highly anisotropic polarizability of this state.^26,61 The vertical dotted line indicates the transition state X = 0 when the resonance of the donor and acceptor electronic states is reached.

7. Protein and water: two “non's” of protein electron transfer

There are obviously two sides, protein's and water's, to the question of what is unique about hydrated proteins as basic machines of life. The specificity of the protein fold is clearly central to its function as a catalyst. One still wonders what are the general properties of proteins distinguishing them from small molecules. Scattering experiments point to a major distinction in a large density of low-frequency vibrations overlapped with a broad quasielastic background from overdamped conformational fluctuations.⁶² These two bands, generic for globular proteins, are typically well separated for small molecules. Returning to the discussion in the Introduction, low-frequency vibrations with frequencies <25 cm⁻¹, which mostly contribute to the observable mean square displacements of the protein atoms,^63,64 fall between fully localized vibrations of small molecules and fully delocalized phonons of solids. This long-ranged cooperativity of atomic displacements not only makes the conformational response of proteins more delocalized than of small molecules, but also dramatically expands the low-frequency wing of the spectrum of protein relaxation. This low-frequency wing often overlaps with the rate constants of biochemical processes and that situation creates the necessity to describe adequately the first of two “non's” of biological electron transport, the non-ergodic reaction kinetics.

The standard definition of the canonical ensemble requires a certain separation of the time-scaled involved. Following Feynman's definition,⁶⁵ the canonical ensemble is defined by requiring that “all the ‘fast’ things have happened and all the ‘slow’ things have not”. This is illustrated in Fig. 13 by placing the observation frequency ω_obs = τ⁻¹_obs in the gap between the fast and slow lumps in the density of relaxation frequencies g(ω). Hydrated proteins present a challenge to this picture by allowing a significant density of relaxation modes at the frequency ω_obs = k of a biochemical process determined by its reaction rate k. The standard canonical prescription fails in this case.

Fig. 13 The definition of the canonical ensemble requires the observation frequency ω_obs = τ⁻¹_obs to be well separated from both the fast and slow modes in the system.⁶⁵ The “equilibrium” canonical average is then performed over all possible configurations of the fast modes. Protein electron transport and biochemical reactions often have rates k = ω_obs well inside a significant density of relaxation modes. The standard canonical prescription fails at these conditions. The dynamically restricted ensemble⁴⁴ [eqn (14)] provides an alternative approach to define the configurational average over the manifold of fast modes in the system.

Since the canonical ensemble is postulated (or derived from the second law of thermodynamics¹⁸) one faces the dilemma of either to change the postulate⁶⁶ or to introduce some dynamical restrictions on the top of the postulated canonical ensembles.⁵⁵ The first approach employs the idea of a restricted canonical ensemble⁶⁶ that limits the phase space available to the system by some set of criteria. Using the dynamical restriction⁴⁴ discussed here amounts to limiting the phase space available to the system by the frequencies exceeding ω_obs. This requirement then modifies the phase-space element Π_j,ωdQ_ω,j to Π_j,ω>ωobsdQ_ω,j and the dynamically restricted ensemble becomes:⁴⁴

where H is the system Hamiltonian.

This formulation leads to eqn (9) when the reorganization energy is concerned and to eqn (10)–(12), for the reaction rate calculations. The existence of slow modes is required to produce a substantial effect on the reaction rate in this formalism, and it is likely the coupling between the protein and water that is responsible for the slow Stokes-shift dynamics⁶⁷ (Fig. 3 and 4). The second “non” of protein electron transfer, the non-Gaussian electrostatics, is where the focus shifts to the water side of the protein–water tandem.

The interest in water as an active medium of bioenergetics has picked up in recent years⁹ driven by two notions: (i) the statistics of the protein fluctuations are highly influenced by their hydration shells³⁶ and (ii) the perturbation of the water structure induced by the protein penetrates surprisingly deeply into the bulk.^48–50 One of course needs to stress that the length-scale of water perturbation is heterogeneous and depends on the property examined. For instance, the experimental data seem to converge to the notion that water's interfacial density profile is perturbed only to a short range,⁹ but the perturbation of the dipolar orientational structure is far more long-ranged (Fig. 7 and 8). The peculiar role of water here is its ability to significantly preserve hydrogen bond structure in the hydration layers, which probably contributes to the creation of dynamical ferroelectric domains observed in simulations of redox-active proteins. It is currently not clear if other polar liquids might possess this ability to form a distinctly new mesophase in the interfacial region.

It appears that both “non's” peculiar to the protein–water duo, the non-ergodic kinetics and the non-Gaussian electrostatics, are employed by living machines to reduce the wasteful release of redox energy into the heat of the phonon bath. The non-ergodic cutoff of nuclear solvation reduces the reaction Gibbs energy ΔG on ultra-fast reactions times of photosynthetic charge separation. The subsequent steps in the electron-transfer chain occur on the nanosecond time-scale where non-ergodic solvation cutoff becomes inefficient and the non-Gaussian statistics of water's fluctuations turns into the leading mechanism of reducing the Gibbs energy lost in electronic hops.

The Gibbs energy lost to heat in a near-activationless transition is given by eqn (4). The released free energy would be then −λ^St = −1.07 eV for activationless reduction of plastocyanin and −λ^St = −1.48 eV for two successive hops of charge separation in the reaction center if the standard picture of linear-response free energies was applied (Fig. 14). In particular for the latter, that Gibbs energy loss not only much exceeds the experimental value of −0.25 eV,⁶⁸ but also puts the charge-separated state below the ground state of the special pair, which is photoexcited by a 1.37 eV photon. As we have discussed in the Introduction, such large losses of the free energy in electronic hops would make the energetics of electronic transport in molecular systems prohibitively inefficient. The very ability of living organisms to transfer the redox energy requires a mechanism reducing the heat losses.

Fig. 14 Energetic scheme of bacterial charge separation in Rb. Sphaeroides. The special pair is photoexcited by a 1.37 eV photon leading to charge separation to H_L in a sequence of two successive hops,⁵² to B_L and then to H_L (Fig. 9). The experimentally measured⁶⁸ Gibbs energy drop is ΔG = −0.25 eV. This number is exceptionally well reproduced by the combination of MD simulations with the non-Gaussian paradigm of electron transfer.²⁶ In contrast, the traditional theories predict a Gibbs energy drop of ΔG = −λ^St = −1.48 eV, placing the charge-separated state below the ground state of the special pair.

Non-Gaussian electrostatics discovered in our simulations shift the energetics of electronic transitions in the desired direction. The energy loss in eqn (4) is significantly reduced due to the breakdown of the linear response resulting in inequality 8. The energetic advantage of the non-Gaussian statistics can be characterized by the number of electron hops in the non-Gaussian mechanism per one hop in the Gaussian picture. Since ΔG = −λ^St applies to an activationless transition in the Gaussian picture, this non-Gaussian parameter becomes

χ_NG = λ^var/λ^St. (15)

This parameter is equal to ≃ 3 for a number of electronic transition in bacterial reaction center (Table 2) and χ_NG ≃ 8 for plastocyanin. It appears that up to an order of magnitude increase of the charge-transfer efficiency can be achieved in protein electron transfer.

Both mechanisms of enhancing the efficiency of electron transport focus on a deeper understanding of the reorganization energy of electronic transitions, the hallmark of the classical Marcus theory of electron transfer.⁵ The collective, highly correlated nature of the dipolar polarization of the protein's hydration shell is behind both the dynamic and statistical anomalies of the protein–water electrostatics. The long tail of the protein Stokes-shift dynamics^67,69 allows an efficient dynamical control of the thermodynamics through the non-ergodic cutoff of polar solvation. No “generic” value⁷⁰ can therefore be assigned to the reorganization energy of protein electronic transitions without specifying the reaction time-window. On the other hand, the non-Gaussian statistics of the electrostatic fluctuations do not allow to describe the activated kinetics with only one reorganization energy. At least two reorganization parameters, λ^St and λ^var, are required. The number of parameters needed for a full description of the global shape of F_i(X) might grow since non-Gaussian electrostatics also imply non-parabolic free-energy surfaces.²⁰

We note in passing that a number of recent studies of heterogeneous exchange of electron between a metal electrode and a protein immobilized on a self-assembled monolayer^71–73 have resulted in very low reorganization energies λ ≃ 0.1–0.4 eV. These values are hard to reconcile with the results of simulations and measurements of homogeneous rates.⁶ However, these data are consistent with the general framework of the non-Gaussian statistics of electronic transitions once one has realized that the effective reorganization energy λ = (λ^St)²/λ^var of eqn (4) is reported by these techniques.

The short-time dynamics and energetics of electron-transfer proteins do not appear to be highly specific and/or optimized. Both show some generic properties of disordered viscous materials, such as an increased density of vibrational states¹⁰ around ≃20–40 cm⁻¹ known as the “boson peak” in structural glasses. This phonon bath is purely vibrational at [small nu, Greek, macron] > 25 cm⁻¹, but becomes increasingly overdamped at lower energies.^74,75 The corresponding ballistic component of solvent reorganization follows the linear response with heat release that cannot be reduced by the non-ergodic cutoff when transitions occur on the picosecond time-scale. The ballistic reorganization energy, ≃0.3–0.4 eV for bacterial charge separation, sets up the lowest free energy penalty for moving electrons in molecular energy chains. Some of this penalty can be reduced by non-Gaussian statistics of electrostatic fluctuations at longer time-scales of transitions. However, phonons affect the non-polar interaction potentials as well. The corresponding reorganization energy, ≃0.12 eV for bacterial charge separation, will always be released to heat in an activationless transition. Bacterial charge separation is therefore pushed into the normal region of electron transfer (X₀₁ ≃ 0.17 eV) to off-set this loss of free energy.

Where specific structural design of a redox protein comes into the picture is in inducing non-Gaussian electrostatic fluctuations. We have found that the appearance of a large dielectric response of the protein–water interface correlates with a large variance of the protein's dipole moment (Table 1). First-shell waters around proteins establish strong hydrogen bonds with the protein surface groups, conceivably even stronger than corresponding hydrogen bonds in the bulk.⁷⁶ Low-frequency vibrations of the protein, modulating the positions of charged surface residues, and with that the protein's dipole moment, will initiate large fluctuations of the hydration layer dipole. A redox protein is then effectively wrapped in the “aluminium foil” of a highly polar hydration layer.

The propensity to form ferroelectric domains might however be an intrinsic property of hydration shells at the surface of a sufficiently large solute.⁷⁷ What is a set of requirements for the formation of a highly polar hydration shell is currently not clear. However, the appearance of this type of response of the protein–water interface resolves several puzzles of bioenergetics. Apart from the overall efficiency of the electron-transfer chains, shallow free energy surfaces of electron transfer satisfying inequality (8) reduce the sensitivity of the activation barrier to the variations of the reaction Gibbs energy ΔG. This feature provides a more robust design of the energy chains since small mutations and variations between different organisms,⁶⁸ as well as environmental conditions, produce smaller effects on the reaction rates.

The last note is on biomimetics.^78,79 What can be learned from studies of natural systems for the design of artificial photosynthesis? The common notion of high optimization of biological energy chains achieved by evolution offers little hope that these designs can be implemented for cheap solar energy harvesting. An unspoken belief behind this notion is that the basic principles of elementary charge transfer are well understood and what remains is to combine several elementary steps in a chain of transitions resulting in a long-lived charge-separated state. This view, implying that natural selection has pushed natural systems into optimizing the reaction Gibbs energies, might be not entirely true.⁸⁰ Natural selection might have operated not only in the space of optimizing the redox potentials, but also in the space of principles of operation. One then wonders if any such principles have been selected which are not a part of the current “standard picture”.

Our studies suggest that such new principles, not seen for small redox molecules, exist for hydrated proteins and have in fact been implemented in natural redox machines. Besides commonly anticipated role of the reaction Gibbs energy, at least three novel types of control of electron transport have been unraveled: (i) adjustments of dynamical time-scales, (ii) non-Gaussian electrostatic fluctuations, and (iii) reorganization energy asymmetry due to highly anisotropic polarizability of the cofactors. It appears that Nature has a broad set of tricks to control the flow of electrons and does not limit its arsenal to two crossing equal-curvature parabolas. Whether these new operational principles are transferable to artificial photosynthetic design remains to be seen.

It also appears that proteins viewed as basic mechanistic elements of biology's energy chains are neither rigid, non-polarizable matrices nor highly tuned molecular machines. They operate on time- and length-scales³ dictated by three major steps of the redox energy transfer: energy input, charge separation, and catalytic reaction (Fig. 1). Like many biological machines operating under conditions of large-scale fluctuations,^81,82 proteins transfer charge by utilizing the broad fluctuation spectrum of the interfacial water with directionality driven by a slight downward gradient of the reaction Gibbs energy in the direction of electron flow.

Acknowledgements

This research was supported by the DOE, Chemical Sciences Division, Office of Basic Energy Sciences (DEFG0207ER15908). We are grateful to Arvi Freiberg for sharing with us the results of ref. 54 and to Dor Ben-Amotz for many useful comments on the manuscript. CPU time was provided by a number of allocations through the TeraGrid Advanced Support Program (TG-MCB080071, TG-MCB080116N, TG-ASC090088).

References

1. D. G. Nicholls and S. J. Ferguson, Bioenergetics 3, Academic Press, London, 2002.
2. A. J. Hoff and J. Deisenhofer, Phys. Rep., 1997, 287, 1.
3. D. Noy, C. C. Moser and P. L. Dutton, Biochim. Biophys. Acta, Bioenerg., 2006, 1757, 90.
4. D. Venkataraman, S. Yurt, B. H. Venkatraman and N. Gavvalapalli, J. Phys. Chem. Lett., 2010, 1, 947.
5. R. A. Marcus and N. Sutin, Biochim. Biophys. Acta, 1985, 811, 265.
6. H. B. Gray and J. R. Winkler, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 3534.
7. D. N. Beratan, J. N. Betts and J. N. Onuchic, J. Phys. Chem., 1992, 96, 2852.
8. Y. Levy and J. Onuchic, Annu. Rev. Biophys. Biomol. Struct., 2006, 35, 389.
9. P. Ball, Chem. Rev., 2008, 108, 74.
10. F. Gabel, D. Bicout, U. Lehnert, M. Tehei, M. Weik and G. Zaccai, Q. Rev. Biophys., 2002, 35, 327.
11. M. Lax, J. Chem. Phys., 1952, 20, 1752.
12. L. D. Zusman, Chem. Phys., 1980, 49, 295.
13. M. Tachiya, Chem. Phys. Lett., 1989, 159, 505.
14. G. King and A. Warshel, J. Chem. Phys., 1990, 93, 8682.
15. A. Warshel and W. W. Parson, Q. Rev. Biophys., 2001, 34, 563.
16. R. A. Marcus, J. Chem. Phys., 1956, 24, 966.
17. M. Marchi, J. N. Gehlen, D. Chandler and M. Newton, J. Am. Chem. Soc., 1993, 115, 4178.
18. L. D. Landau and E. M. Lifshits, Statistical Physics, Pergamon Press, New York, 1980.
19. J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, Amsterdam, 2003.
20. D. V. Matyushov, Acc. Chem. Res., 2007, 40, 294.
21. D. N. LeBard and D. V. Matyushov, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 061901.
22. D. N. LeBard and D. V. Matyushov, J. Phys. Chem. B, 2010, 114, 9246.
23. D. V. Matyushov and R. Schmid, J. Chem. Phys., 1996, 104, 8627.
24. H. C. Andersen, D. Chandler and J. D. Weeks, Adv. Chem. Phys., 1976, 34, 105.
25. D. N. LeBard, V. Kapko and D. V. Matyushov, J. Phys. Chem. B, 2008, 112, 10322.
26. D. N. LeBard and D. V. Matyushov, J. Phys. Chem. B, 2009, 113, 12424.
27. D. Chandler, Nature, 2005, 437, 640.
28. C. Y. Lee, J. A. McCammon and P. J. Rossky, J. Chem. Phys., 1984, 80, 4448.
29. P. K. Ghorai and D. V. Matyushov, J. Phys. Chem. A, 2006, 110, 8857.
30. M.-L. Tan, E. Dolan and T. Ichiye, J. Phys. Chem. B, 2004, 108, 20435.
31. J. Blumberger and M. Sprik, Theor. Chem. Acc., 2006, 115, 113.
32. D. N. LeBard and D. V. Matyushov, J. Phys. Chem. B, 2008, 112, 5218.
33. K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics, Springer-Verlag, Berlin, 1992.
34. F. G. Parak, Rep. Prog. Phys., 2003, 66, 103.
35. W. Doster, Eur. Biophys. J., 2008, 37, 591.
36. H. Frauenfelder, G. Chen, J. Berendzen, P. W. Fenimore, H. Jansson, B. H. McMahon, I. R. Stroe, J. Swenson and R. D. Young, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 5129.
37. W. Doster, S. Cusack and W. Petry, Nature, 1989, 337, 754.
38. P. W. Fenimore, H. Frauenfelder, B. H. McMahon and R. D. Young, Proc. Natl. Acad. Sci. U. S. A., 2004, 101, 14408.
39. Y. He, P. I. Ku, J. R. Knab, J. Y. Chen and A. G. Markelz, Phys. Rev. Lett., 2008, 101, 178103.
40. J. P. Boon and S. Yip, Molecular Hydrodynamics, Dover Publications, Inc., New York, 1991.
41. R. M. Daniel, J. L. Finney and J. C. Smith, Faraday Discuss., 2003, 122, 163.
42. S. Khodadadi, A. Malkovskiy, A. Kisliuk and A. P. Sokolov, Biochim. Biophys. Acta, Proteins Proteomics, 2010, 1804, 15.
43. D. V. Matyushov, J. Chem. Phys., 2005, 122, 084507.
44. D. V. Matyushov, J. Chem. Phys., 2009, 130, 164522.
45. R. Jimenez, G. R. Fleming, P. V. Kumar and M. Maroncelli, Nature, 1994, 369, 471.
46. B. K. P. Scaife, Principles of dielectrics, Clarendon Press, Oxford, 1998.
47. A. W. Omta, M. F. Kropman, S. Woutersen and H. J. Bakker, Science, 2003, 301, 347.
48. R. Giordano, G. Salvato, F. Wanderlingh and U. Wanderlingh, Phys. Rev. A: At., Mol., Opt. Phys., 1990, 41, 689.
49. S. Perticaroli, L. Comez, M. Paolantoni, P. Sassi, L. Lupi, D. Fioretto, A. Paciaroni and A. Morresi, J. Phys. Chem. B, 2010, 114, 8262.
50. S. Ebbinghaus, S. J. Kim, M. Heyden, X. Yu, U. Heugen, M. Gruebele, D. M. Leitner and M. Havenith, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 20749.
51. B. Born, S. J. Kim, S. Ebbinghaus, M. Gruebele and M. Havenith, Faraday Discuss., 2009, 141, 161.
52. W. Zinth and J. Wachtveitl, ChemPhysChem, 2005, 6, 871.
53. G. R. Fleming, J. L. Martin and J. Breton, Nature, 1988, 333, 190.
54. K. Timpmann, A. Ellervee, A. Laisaar, M. R. Jones and A. Freiberg, Ultrafast processes in spectroscopy, 1998, p. 236.
55. H. Sumi and R. A. Marcus, J. Chem. Phys., 1986, 84, 4894.
56. H. Wang, S. Lin, J. P. Allen, J. C. Williams, S. Blankert, C. Laser and N. W. Woodbury, Science, 2007, 316, 747.
57. C. W. Gardiner, Handbook of Stochastic Methods, Springer, Berlin, 1997.
58. P. F. Barbara, T. J. Meyer and M. A. Ratner, J. Phys. Chem., 1996, 100, 13148.
59. C. Kirmaier and D. Holten, in Photosynthetic Bacterial Reaction Center: Structure and Dynamics, ed. J. Breton and A. Vermégio, Plenum, New York, 1988, vol. 149, pp. 219–228.
60. M. Bixon and J. Jortner, Adv. Chem. Phys., 1999, 106, 35.
61. D. V. Matyushov and G. A. Voth, J. Chem. Phys., 2000, 113, 5413.
62. J. C. Smith, Q. Rev. Biophys., 1991, 24, 227.
63. N. Go, T. Noguti and T. Nishikawa, Proc. Natl. Acad. Sci. U. S. A., 1983, 80, 3696.
64. F. G. Parak and K. Achterhold, J. Phys. Chem. Solids, 2005, 66, 2257.
65. R. P. Feynman, Statistical Mechanics, Westview Press, Boulder, CO, 1998.
66. R. G. Palmer, Adv. Phys., 1982, 31, 669.
67. T. Li, A. A. Hassanali and S. J. Singer, J. Phys. Chem. B, 2008, 112, 16121.
68. M. Volk, G. Aumeier, T. Langenbacher, R. Feick, A. Ogrodnik and M.-E. Michel-Beyerle, J. Phys. Chem. B, 1998, 102, 735.
69. S. K. Pal and A. H. Zewail, Chem. Rev., 2004, 104, 2099.
70. C. C. Moser, C. C. Page and P. L. Dutton, Philos. Trans. R. Soc. London, Ser. B, 2006, 361, 1295.
71. J. Hirst and F. A. Armstrong, Anal. Chem., 1998, 70, 5062.
72. A. Kranich, H. Naumann, F. P. Molina-Heredia, H. J. Moore, T. R. Lee, S. Lecomte, M. A. de la Rosa, P. Hildebrandt and D. H. Murgida, Phys. Chem. Chem. Phys., 2009, 11, 7390.
73. D. E. Khoshtariya, T. D. Dolidze, M. Shushanyan, K. L. Davis, D. H. Waldeck and R. van Eldik, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 2757.
74. M. Diehl, W. Doster, W. Petry and H. Schober, Biophys. J., 1997, 73, 2726.
75. M. Marconi, E. Cornicchi, G. Onori and A. Paciaroni, Chem. Phys., 2008, 345, 224.
76. F. Demmel, W. Doster, W. Petry and A. Schulte, Eur. Biophys. J., 1997, 26, 327.
77. A. D. Friesen and D. V. Matyushov, arXiv1004.1728v2.
78. J.-M. Lehn, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 4763.
79. B. Brushan, Philos. Trans. R. Soc. London, Ser. A, 2009, 367, 1445.
80. C. C. Page, C. C. Moser and P. L. Dutton, Curr. Opin. Chem. Biol., 2003, 7, 551.
81. P. Hänggi and F. Marchesoni, Rev. Mod. Phys., 2009, 81, 387.
82. F. C. MacKintosh and C. F. Schmidt, Curr. Opin. Cell Biol., 2010, 22, 29.

---

Footnote

† Electronic supplementary information (ESI) available: Calculations of the temperature and pressure effects on the rate of bacterial charge separation. See DOI: 10.1039/c0cp01004a

---