A 300-fold conductivity increase in microbial cytochrome nanowires due to temperature-induced restructuring of hydrogen bonding networks.

Although proteins are considered as nonconductors that transfer electrons only up to 1 to 2 nanometers via tunneling, Geobacter sulfurreducens transports respiratory electrons over micrometers, to insoluble acceptors or syntrophic partner cells, via nanowires composed of polymerized cytochrome OmcS. However, the mechanism enabling this long-range conduction is unclear. Here, we demonstrate that individual nanowires exhibit theoretically predicted hopping conductance, at rate (>1010 s-1) comparable to synthetic molecular wires, with negligible carrier loss over micrometers. Unexpectedly, nanowires show a 300-fold increase in their intrinsic conductance upon cooling, which vanishes upon deuteration. Computations show that cooling causes a massive rearrangement of hydrogen bonding networks in nanowires. Cooling makes hemes more planar, as revealed by Raman spectroscopy and simulations, and lowers their reduction potential. We find that the protein surrounding the hemes acts as a temperature-sensitive switch that controls charge transport by sensing environmental perturbations. Rational engineering of heme environments could enable systematic tuning of extracellular respiration.

INTRODUCTION

Electron flow between cofactors in proteins is central to many life processes such as respiration, photosynthesis, and nitrogen fixation (). The protein architecture between cofactors has long been considered a passive tunneling bridge resulting in transport that decays exponentially with distance, thus limiting transport to 1 to 2 nm (). Some proteins can conduct over distances up to 20 nm (). Long-distance (>20 nm) charge transport is usually achieved via a redox potential gradient generated by distinct redox enzymes, which provides an energetic driving force for the transport of charge toward the terminal acceptor (). However, the common soil bacterium Geobacter sulfurreducens transports electrons extracellularly over micrometer distances via “microbial nanowires.” These micrometer-long filaments are homopolymers of cytochromes OmcS (Fig. 1) (, ) or OmcZ (), which contrasts with multiple distinct cytochromes used by other bacteria for extracellular electron transfer up to nanometers (). This extracellular electron transfer is important in global biogeochemical cycling of carbon, nutrients, and metals as well as in regulating the bacterial release of methane to the atmosphere ().

Fig. 1.

Closely stacked hemes in OmcS nanowires provide a continuous path for extracellular electron transport over micrometers.

(A) Transmission electron microscopy image showing a G. sulfurreducens producing OmcS nanowires. Scale bar, 200 nm. Image levels were adjusted selectively in different areas to show the most information. (B) Each pair of parallel-stacked hemes is perpendicular to the next pair in a T-shaped geometry. (C) Upon polymerization, the hemes align to stack along the entire length of the OmcS filament with a 20-nm helical pitch. Each OmcS monomer and corresponding hemes are highlighted in different colors.

The atomic structure of OmcS nanowires revealed seamless stacking of heme cofactors, which provide a continuous path for electron flow (Fig. 1, B and C) (, ). Prior theoretical studies of multiheme cytochromes have assumed that charges hop between heme cofactors (), with the heme-to-heme hopping rates determined by the reduction potentials of the individual hemes and the electronic coupling between each pair. However, each protomer in OmcS must have a similar midpoint redox potential. Therefore, while there may exist a redox potential gradient within a protomer, there remain questions regarding how G. sulfurreducens generates high enough current through such a structure to meet its respiratory needs.

Prior computational studies on OmcS nanowires suggested multiple mechanisms such as hopping (), quantum transport (), and coherence-assisted transport () to explain the experimentally measured conductivity (). However, experimental evidence supporting a particular mechanism is lacking.

In this work, we measure the intrinsic (contact-free) conductivity of individual OmcS nanowires as a function of nanowire length, voltage, temperature, and pH, which reveal an unexpectedly high hopping rate (>1010 s−1), comparable to synthetic molecular wires, with negligible carrier loss over micrometers. We also find a 300-fold increase in nanowire conductance upon cooling, which vanishes upon deuteration. Computations show that cooling causes a massive rearrangement of hydrogen (H) bonding networks, making hemes more planar, thus lowering their reduction potential and increasing conductivity.

RESULTS

Correlating nanowire structure with function using conductivity measurements under conditions used to solve the atomic structure

To correlate nanowire structure with conductivity, we performed measurements of individual nanowires in the same alkaline (pH 10.5) buffer used to solve the atomic structure, rather than pH 7 used previously for conductivity measurements (). Using atomic force microscopy (AFM), we identified OmcS nanowires bridging gold electrodes via the distinct height and helical pitch of OmcS nanowires (Fig. 2). The nanowires showed high conductance over 300 nm at physiologically relevant potentials (Fig. 2, D and B) (). This finding contrasts with most proteins that can transfer electrons only up to 1 to 2 nm, due to the exponential decay of current with distance, characteristic of the tunneling mechanism (). Conductance in noncytochrome proteins has been limited to 20 nm (). In contrast, conductance of OmcS nanowires is more comparable to synthetic molecular systems (–), for which the transport mechanism has been observed to transition from tunneling to hopping for molecular lengths greater than 4 nm ().

Fig. 2.

Correlating nanowire structure with conductivity by measurement under conditions used to solve the atomic structure.

(A) Atomic force microscopy (AFM) image of a nanowire bridging the gap between nanoelectrodes. Scale bar, 100 nm. (B) Nanowire height and (C) pitch measured at locations shown in (A). (D) Current-voltage (I-V) profile for a single nanowire hydrated in the buffer, compared to the buffer alone, measured at ~290 K, and at pH 10.5 used to solve the OmcS structure.

Intrinsic conductivity measurements reveal hopping transport with negligible carrier loss

To test whether OmcS transports electrons via hopping, we measured dc (direct current) conductance as a function of nanowire length using conducting probe AFM (CP-AFM) at physiological pH (). We deposited nanowires across a gold/SiO2 interface (Fig. 3A) and electrically contacted the nanowires at varying distances away from the gold electrode with a conductive AFM tip to measure the electrical resistance (Fig. 3B). We found that the nanowire resistance increased linearly with length (Fig. 3C). Such a linear dependence of resistance on length is consistent with a multistep hopping mechanism (). These length-dependent measurements also helped to determine the contact resistance (RC ~ 1.2 gigaohms), yielding the intrinsic conductivity of OmcS nanowires (~24 mS/cm) where S stands for Siemens (Fig. 3C), comparable to prior studies (, , ).

Fig. 3.

Identification of hopping conductivity in nanowires.

(A) Measurement strategy for the length-dependence of the single nanowire resistance, at 300 K and at pH 7, using CP-AFM. The current is measured between the gold electrode and AFM tip as the tip contacts the nanowire at varying distances away from the electrode. Scale bar, 100 nm. (B) Representative I-V curves along a nanowire at various distances. Inset: I-V collected on the gold electrode and the insulating SiO2 substrate (black). (C) Resistance versus nanowire length. Red dashed line is a linear fit (R2 = 0.74). (D) CP-AFM data plotted as flux (at 0.1 V) versus the inverse of the number of hopping steps. Red dashed line is a linear fit (R2 = 0.80). (E) The product of the charge flux and the number of hopping steps plotted against the number of hopping steps. The slope of the linear fit indicates the rate of carrier loss. All error bars represent SEM. n = 9.

To further understand the hopping conductivity, we analyzed our data using an expression for hopping flux (J) of the form J = (khop/N)e−α, where khop is the hopping rate, α is related to carrier decay, and N is the number of heme-to-heme hopping steps (). Using the bias of 0.1 V and the rise per protomer of 4.67 nm (axial distance between protomers), we converted the measured resistance versus distance into flux (J) versus 1/N (Fig. 3D) [J = 0.1 V/(R ∙ ℯ); N = (d/4.67 nm) ∙ (6 − 1) hemes]. Our measurements showed that α is close to zero, suggesting negligible loss of charge carriers over micrometer distances (Fig. 3E). Furthermore, the measured hopping rate in nanowires (3.2 × 1010 s−1) is comparable to that in synthetic porphyrin wires (). This high hopping rate and negligible carrier loss in nanowires explains unique ability of G. sulfurreducens to transport electrons efficiently over micrometers.

A 300-fold increase in nanowire conductivity upon cooling

Further insight into the conduction mechanism was obtained by measuring the intrinsic conductivity of individual nanowires as a function of temperature at pH 10.5. To eliminate artifacts due to contacts or electrode polarization, we used a four-point probe (4-probe) configuration (). Upon cooling from room temperature (300 K), the conductivity of the nanowires initially increased exponentially by approximately 300-fold, followed by a transition to exponentially decreasing conductivity (Fig. 4A). Notably, this behavior was reversible and did not show hysteresis. This finding is unexpected because instead of increasing, hopping conductivity decreases exponentially upon cooling, known as Arrhenius behavior. Thus, individual OmcS nanowires showed a temperature dependence of conductance similar to living biofilms of a G. sulfurreducens strain CL1 () engineered to overexpress OmcS ().

Fig. 4.

OmcS demonstrates a non-Arrhenius temperature dependence.

(A) Temperature dependence of the conductivity of OmcS nanowires determined experimentally, by 4-probe measurements on nanoelectrodes, at pH 10.5 used to solve the OmcS structure. The red dashed lines are exponential fits to the data. Error bars represent SD. The data are representative of four independent measurements. (B) Temperature dependence of the conductance of a film of deuterated OmcS nanowires. The measurements were performed in a 2-probe configuration on nanoelectrodes with 300-nm spacing. (C) Temperature dependence of the conductance of a film of air-oxidized and chemically reduced OmcS nanowires. The film was measured on a gold interdigitated array (IDA) with 5-μm spacing. (D) UV-visible spectra of OmcS nanowires at varying temperature. The spectrum of chemically reduced hemes is displayed in gray for reference. ITO, indium tin oxide. (E) Temperature dependence of the UV-visible spectrum of OmcS nanowires under bias on a transparent FTO IDA with 5-μm spacing. a.u., absorbance units.

Measurements of nanowire films in a 2-probe configuration also showed the non-Arrhenius temperature dependence, albeit with a deviation from the trend between 280 and 300 K, likely due to contact resistance (Fig. 4B). Notably, the non-Arrhenius trend vanished upon deuterating the nanowires (Fig. 4B). The conductivity was weakly dependent on temperature and showed a sharp transition around 220 K, likely due to the protein glass transition (). Therefore, this low-temperature transition is distinct from the conductivity increase near physiological temperatures. The deviation from the trend at higher temperatures can likely be explained by the difference in zero-point energy for H and deuterium (D) bonds. Deuterium substitution is known to change the strength of H-bonds, with weaker H-bonds transitioning to stronger D-bonds (–). Stronger H-bonds experience a smaller increase in strength upon deuterium substitution than weaker H-bonds (). This change in bond strength results in a smoother free energy landscape (). The vanishing of the increase in conductivity upon cooling after deuterium exchange suggests an important role for H-bonding in nanowire conductivity, as we reported for amyloid proteins ().

Our findings contrast with previous temperature-dependent conductivity measurements of proteins, where deuteration shows a stronger temperature dependence due to dominant electron transport through H-bonds and the increased rigidity of the mostly β sheet protein, which increases the reorganization energy (). OmcS contains a very small number of β strands (<6%), suggesting a previously unknown role of H-bonding in modulating protein conductivity.

Our findings also contrast previous reports of non-Arrhenius temperature dependence in films of cytochrome c3 and hydrogenase (), where the conductivity increase upon cooling was due to an increased tendency of hydrogenase to reduce cytochrome c3 at lower temperatures. In contrast, the ultaviolet (UV)–visible spectrum of air-oxidized OmcS nanowires did not change upon cooling (Fig. 4D), suggesting an absence of substantial change in the redox state of OmcS. Although reduced films showed higher conductance (Figs. 4C and 5E), both reduced films and air-oxidized films showed a non-Arrhenius trend (Fig. 4C). Even under current carrying conditions, only a small, reduced population was evident (Fig. 4E). Therefore, the cooling-induced increase in conductivity is intrinsic to OmcS and not due to a change in the redox state of the hemes.

Fig. 5.

Computations and experiments reveal facile hole transport in OmcS nanowires.

(A) Free energy landscape versus standard hydrogen electrode (SHE) for electron and hole transfer in OmcS obtained by QM/MM electronic structure calculations. Error bars are SEM. (B) Reorganization energies for electron and hole transfer models of OmcS. Solid markers were computed using QM/MM-derived energies. Hollow markers were computed from molecular mechanics–derived energies. Error bars are SEM. (C) Electronic coupling between heme pairs as a function of minimum edge-to-edge distance. The dashed lines are exponential fits to the data. The shaded areas differentiate the electronic coupling values corresponding to the parallel stacked heme pairs and the T-shaped heme pairs. Error bars are SEM. (D) Computed charge mobility and conductivity. Solid bars were computed using QM/MM-derived reorganization energies, and striped bars were computed using MM-derived reorganization energies. (E) Conductance comparison of an air-oxidized versus chemically reduced film of OmcS nanowires reveals higher hole conductivity. Error bars are SEM.

Computational modeling captures experimentally measured conductivity

To obtain a molecular level description of the OmcS charge transport mechanism, we complemented our experimental results with computational modeling. Using the near-atomic resolution cryo–electron microscopy structure of OmcS nanowires () as a starting structure, we used the CHARMM36 () force field to build 12-heme (dimer) models of OmcS under fully oxidized (electron transfer regime) and fully reduced (hole transfer regime) conditions. Molecular dynamics (MD) simulations of these models at 310 K provided sampling of thermally accessible configurations of OmcS. We used hybrid quantum/classical [quantum mechanics (QM)/molecular mechanics (MM)] calculations on an ensemble of these configurations to compute heme ionization energies (details in Materials and Methods). From these energies, we estimated the reduction potential of each heme.

We found that the six hemes in OmcS have a wide range of reduction potentials (Fig. 5A and tables S1 and S2), similar to QM/MM studies of other multiheme cytochromes (). Our results differ from prior theoretical estimates of differences in heme potentials in OmcS (), which predicted a narrower range of values using continuum electrostatics calculations on an energy minimized conformation. When we replaced the explicit protein/solvent environment with a polarizable dielectric continuum model, the range of computed reduction potentials was diminished and resulted in a distinct free energy landscape (fig. S2). This result demonstrates the importance of the protein electrostatics in tuning the reduction potentials of c-type hemes. This result also suggests that our methodology may be particularly sensitive to polarization of the electron density in response to the protein environment.

We used the same ionization energies from our QM/MM calculations to estimate reorganization energies (Fig. 5B). We also computed reorganization energies from MM-derived energies, using a framework that was designed to address nonergodic protein dynamics or conformational sampling that occurs on time scales slower than the electron transfer reaction (Fig. 5B) (). Given substantial debate within the literature regarding the necessity of this approach (, ), we present our results as a direct comparison to experiments. In general, the MM-derived reorganization energies were smaller and, thus, give rise to greater charge transfer rates (tables S5 and S6) and higher mobility estimates (Fig. 5D). These values are closer to experimental values. Therefore, further comparisons between computed and experimentally determined conductivity will use the MM-computed reorganization energies.

Calculations of the electronic coupling between adjacent heme pairs demonstrated that coupling is much less than the reorganization energy (Fig. 5, B and C, and tables S3 and S4). Therefore, the electron transfer reactions are nonadiabatic and can be described by Marcus theory (). The values that we obtained for the coupling were higher for the parallel-stacked heme pairs than for the T-shaped heme pairs (Fig. 5C). This is expected as the parallel orientation of the porphyrin rings allows closer distances between hemes and increases overlap between the pi-orbitals. An exponential fit for the donor-acceptor coupling (|HAD|) of the form ∣HAD∣ = A exp [ − β(r − r0)/2] (where β is the distance decay constant, r is the minimum edge-to-edge distance between donor and acceptor, r0 is the van der Waals distance of 3.4 Å, and A is a pre-exponential factor) yielded distance decay factors of 0.69 ± 0.3 Å−1 and 0.59 ± 0.3 Å−1 for electron transfer and hole transfer couplings, respectively. These decay constants are smaller than previously reported values for multiheme cytochromes (), due to the high coupling between T-stacked hemes 1 and 2 (9.1 meV), likely due to the proximity of cysteine 12 to heme 2 (fig. S3) ().

Using these computed reduction potentials, reorganization energies, and electronic couplings in Marcus theory, we computed charge transfer rates between each heme using the hopping model (tables S5 and S6) (). From these values, we estimated electron and hole mobilities (μ) using kinetic Monte Carlo simulations (Fig. 5D). Our modeling predicted the hole conductivity to be approximately an order of magnitude higher than electron conductivity (Fig. 3D). This result is consistent with prior studies () and with our finding that a chemically reduced film of OmcS nanowires had ~6-fold higher conductance than air-oxidized nanowires (Figs. 4C and 5E). The increased conductivity under reduced conditions may have physiological implications for the anoxic growth environment of G. sulfurreducens, allowing rapid extracellular transport of hole carriers to meet metabolic needs.

Cooling decreases heme reduction potentials, which could explain increased conductivity

Non-Arrhenius hopping conductivity suggests a temperature-dependent activation energy, indicating a dynamic, temperature-sensitive, free energy landscape for OmcS nanowires. To determine the effect of cooling on the reduction potentials of OmcS, we recomputed heme reduction potentials and reorganization energies for electron transfer using an ensemble of configurations from MD simulations at 270 K. We found that cooling decreased reduction potentials for hemes 1, 2, and 3, with heme 3 showing the largest shift in potential (Fig. 6A and table S7). Changes in reorganization energy were relatively small (<0.1 eV) (Fig. 6, A and B, and tables S1 and S7), and the electronic coupling did not change with temperature (Fig. 6C). Thus, the reduction potentials dominate the temperature dependence of the activation energy in OmcS. The computed charge mobility was ~77-fold higher at 270 K (Fig. 6D), in quantitative agreement with the experiments (Fig. 3A), validating our computational methodology.

Fig. 6.

Temperature-induced decrease in reduction potentials enhances nanowire conductivity.

(A) Computed free energy landscape for electron transfer at 310 and 270 K. (B) Computed reorganization energies for OmcS models at 310 and 270 K. (C) Computed electronic coupling between heme pairs at 310 and 270 K. (D) Computed temperature dependence of charge mobility and conductivity in OmcS.

Cooling massively restructures H-bonding network in OmcS nanowires

We next performed a comparative study of our two models at 310 and 270 K, with the goal of identifying the features of OmcS that allow its conductivity to be tuned over such a large range. The loss of the non-Arrhenius temperature dependence in response to deuterium exchange (Fig. 4B) indicated the role of H-bonding network in nanowire conductivity. Temperature is reported to influence the strength of H-bonds (). We evaluated whether the temperature alters H-bonds within OmcS, which could modify protein conformation and the electronic structure of hemes (–).

On average, our OmcS dimer model has 6.2% more H-bonds at 270 K relative to 310 K (Fig. 7A). However, this value alone fails to capture the change in composition of the H-bond network (i.e., the distribution of donor-acceptor pairs). To quantify the difference in composition, we computed a H-bonding frequency matrix, whose elements A(i, j) described the frequency at which a H-bond was realized between residues i and j. From this matrix, we use the norm of the matrix to define the characteristic H-bond frequency (CHF). When applied to the H-bond frequency matrix of a particular simulation, this metric provides information regarding the number of H-bonds present at any point in time. Therefore, the CHF for the simulation at 270 K is larger than the CHF for the simulation at 310 K (Fig. 7B). We also applied this analysis to the difference matrix, whose elements are composed of the differences of the matrix elements computed at 270 K and those computed at 310 K. Each difference matrix element provides a direct measure of how prevalent a particular H-bonding interaction is at each temperature. The difference CHF therefore provides a quantification of the difference in H-bonding structure at the two temperatures.

Fig. 7.

Temperature induces restructuring of hydrogen bonds.

(A) Average number of H-bonds in a selection of 200 snapshots from MD simulations of an OmcS dimer model. ****P < 0.0001. (B) Characteristic H-bonding frequency (norm of the H-bonding frequency matrix) for the h-bonding frequency matrices computed at 310 and 270 K. (C) H-bonding frequency difference matrices for the hemes and protein within 5 Å of each heme. Green bars represent the difference CHF for the displayed matrices. Representative H-bonds are shown in the vicinity of heme 2 at (D) 270 K and (E) 310 K.

We found that the difference CHF was large, as its value was more than half the CHF of the individual models and larger than the difference between the CHF values computed at 270 and 310 K (Fig. 7B). This is indicative of a massive restructuring of the H-bond network upon cooling. Decomposition of the individual H-bond interactions revealed that only 7% of these interactions were unique to a particular temperature. Rather, 59% of the interactions were more frequent at 270 K (fig. S4). Thus, the H-bonding structure of OmcS is dominated by a different set of donor-acceptor pairs at lower temperatures. This result suggests that there exist distinct states of the H-bonding network, defined by local minima in the H-bonding free energy landscape. This interpretation is consistent with the weak temperature dependence that we observed with deuterated OmcS (Fig. 4B). Deuteration creates a smoother free energy landscape. Therefore, the preference for particular H-bonding states will be less defined, resulting in less contrast between the H-bonding networks at high and low temperatures. The increase in H-bonding frequency upon cooling that we observed in our model is typical of secondary structures such as α helices and β sheets (). However, OmcS has no β sheets and only 13% of helices (). Therefore, the restructuring of the H-bonding network and its role in nanowire conductivity may be unique to OmcS.

Changes to the H-bonding structure local to the hemes will have the largest effect on heme reduction potentials. To focus attention on the H-bonding interactions most likely to alter heme energetics, we computed the frequency difference matrices for residues within 5 Å of each heme (Fig. 7C). Notably, the H-bond network around heme 4, whose predicted reduction potential remained unchanged at 270 K, was most robust under changing temperature. The remainder of the hemes experienced various changes to the H-bond network in their vicinity.

For example, an additional H-bond with the propionate group of heme 2 is realized at 270 K (Fig. 7, D and E). However, the complexity of the H-bond networks surrounding the hemes makes it difficult to assign shifts in reduction potential to individual H-bond interactions. Therefore, we analyzed how H-bonding–induced changes in the protein structure surrounding the hemes affected heme conformation and electrostatics.

Cooling makes hemes more planar, as revealed by Raman spectroscopy and simulations

Heme conformation, particularly heme ruffling (Fig. 8A), has been demonstrated to affect heme reduction potential (, ), with increased ruffling corresponding to lower reduction potentials. Ruffling distortions are induced by the structure of the protein surrounding the hemes (, ) and are therefore readily influenced by the H-bond network. Quantification of the ruffling distortion from our MD simulations revealed that, on average, the ruffling distortion was lower at 270 K than at 310 K (Fig. 8B). This can readily be explained by the thermal depopulation of the heme ruffling mode. However, the cooling-induced shift toward heme planarity was not uniform, as hemes 2, 3, and 6 did not exhibit decreased ruffling. Therefore, in addition to thermal depopulation, changes to the structure of the protein surrounding the hemes appear to have altered heme ruffling. Notably, the computed shifts in ruffling correlated with the computed shifts in reduction potential [coefficient of determination (R2) = 0.65; Fig. 8C].

Fig. 8.

Hemes become more planar upon cooling.

(A) Illustration of the heme ruffling distortion. (B) Computed ruffling component of the out-of-plane distortion of the six hemes in OmcS from our models simulated at 270 and 310 K. Error bars are SEM. Data are from an ensemble of 200 snapshots from our MD trajectory. (C) Computed cooling-induced shifts in the reduction potentials and ruffling distortions plotted against each other. Points are labeled to indicate heme identity. Error bars are SEM. Dashed line is a linear fit to the data (R2 = 0.65). (D) Raman spectra of OmcS nanowires. Data were collected under warm (~296 K) and cold (~280 K) conditions. (E) Magnification of the ν2 peak from the Raman spectra in (D). (F) Magnification of the ν10 peak from the Raman spectra in (D). (G) Lorentzian decomposition of the ν10 peak of the Raman spectrum collected under warm (~296 K) conditions. Blue dashed line is the multipeak fit to the data. (H) Lorentzian decomposition of the ν10 peak of the Raman spectrum collected under cold (~280 K) conditions. Blue dashed line is the multipeak fit to the data. (I) Peak positions corresponding to the Lorentzian fit to the ν10 peak of the Raman spectra in (G) and (H).

We further probed heme ruffling experimentally with Raman spectroscopy. The ν2, ν3, and ν10 modes have been associated with heme out-of-plane distortions () and have been directly correlated with ruffling. We collected the Raman spectra at room temperature (~296 K) and under cooled conditions (T ~ 280 K) (Fig. 8D). Both the ν2 and ν10 peaks showed a shift toward higher frequency (Fig. 8, E and F), consistent with a decrease in ruffling (). The ν3 peak was too noisy under cooled conditions and therefore could not be included in our analysis. Of the ν2 and ν10 peaks, ν10 is the most sensitive to nonplanar distortions of the hemes (). We performed Lorentzian decomposition analysis on the ν10 peak to extract information regarding the ruffling of individual hemes (Fig. 8, G and H). Comparison of the decomposition analysis with the computed ruffling distortions revealed agreement between experiment (Fig. 8B) and simulation (Fig. 8I), further suggesting the role of heme distortions in modulating nanowire conductivity.

Cooling increases the electric field at the heme iron, lowering reduction potential

In addition to change in heme distortions upon cooling, we also examined the effect of temperature on electrostatic potential of hemes and protein-induced electric field at the heme iron centers. We found that the direction of these fields was altered upon changing the temperature from 310 K (red arrow) to 270 K (blue arrow) (Fig. 9A). We also found that cooling increased the magnitude of heme electric field, thereby decreasing its reduction potential, in agreement with prior studies using external electric fields (, ). We observed a negative correlation between the computed shift in the electric field magnitude and the computed shift in reduction potential (R2 = 0.51; Fig. 9B). In contrast, the hemes experienced relatively small (<20 mV) shifts in electrostatic potential (Fig. 9C), which did not correlate with the shifts in reduction potential as expected (Fig. 9D). This finding suggests that the change in H-bonding structure altered the orientation of charged residues without significantly changing their distance from the hemes. Therefore, we conclude that a combination of heme ruffling and protein-induced electric field is responsible for the predicted shifts in heme reduction potentials and thus the non-Arrhenius temperature dependence (Fig. 10).

Fig. 9.

Temperature induced changes in protein structure alter the electric field at the heme Fe atoms.

(A) Illustrations of the vectors representing the electric field (E-field) at the iron atoms of the hemes in OmcS. Vectors displayed in blue were computed from our simulation at 270 K, and the vectors displayed in red were computed from our simulation at 310 K. Insets plot the magnitude of the electric field for the simulations run at 270 and 310 K. Error bars are SD. (B) Cooling-induced shifts in heme reduction potentials and electric field magnitudes plotted against each other. Points are labeled to indicate heme identity. Error bars are SEM. Dashed line is a linear fit to the data (R2 = −0.51). (C) Electrostatic potential at each heme site, computed from our simulations at 270 and 310 K. Error bars are SD. (D) Cooling-induced shifts in heme reduction potentials and electrostatic potentials plotted against each other. Points are labeled to indicate heme identity. Dashed line is a linear fit to the data (R2 = 0.60). ns, not significant. *P < 0.05; **P < 0.01; ****P < 0.0001.

Fig. 10.

A 300-fold increase in conductivity of OmcS nanowires due to temperature-induced restructuring of H-bonding networks.

Cooling restructures H-bonds, which makes hemes planar and increases the electric field at the heme iron. These effects collectively lower the heme reduction potential to increase nanowire conductivity.

DISCUSSION

Electron transfer via c-type cytochrome plays a central role in many chemical and biological processes. However, cytochromes were long-considered monomeric, limiting transfer to over a few nanometers. Here, we show that common soil bacteria Geobacter can transport electrons over micrometers via polymerized cytochrome nanowires with negligible electron loss. Combining theory with measurements of intrinsic, contact-free, conductivity measurements of individual nanowires as a function of nanowire length, voltage, and temperature, we find that a hopping mechanism explains the unexpected length and temperature dependence of nanowire conductivity. This was previously reported only in synthetic nanowires and nonredox proteins where it was limited to distances less than 20 nm (, ).

Alternative models based on quantum transport have introduced the possibility that such a mechanism has a role in OmcS (, ). To date, all studies of quantum transport in OmcS have not incorporated the effect of the protein environment. Such a heterogeneous environment with thermal fluctuations would greatly affect the stability of delocalized electronic states. Our calculations suggest that the protein has an important role in curating unique environments for each heme. Further studies are required to evaluate quantum transport in OmcS nanowires.

Protein environment controls the reduction potentials of c-type hemes (). We demonstrate that this protein control of reduction potentials allows nanowire conductivity to be tuned over a remarkable range. Thus, the OmcS nanowires are a previously unknown class of metalloproteins with highly tunable electronic properties. Our studies can help guide systematically fine-tuning the conductivity of nanowires by modulating protein-induced heme distortions. Using a suite of complementary experiments to measure the conductivity of individual nanowires as a function of molecular length and temperature, and combining with computational studies using MD and QM/MM electronic structure calculations, here we show that the thermally activated vibrations induce out-of-plane heme distortions that modulate activation energy to enhance conductivity by 300-fold. We thus present a strategy to induce cooperative and large-scale conformational changes that modulate heme distortions, which were previously thought to have only local influence.

Our findings have consequences for the physiology of G. sulfurreducens (). OmcS is essential for the reduction of iron oxides, an important electron acceptor in the native environment of G. sulfurreducens, as well during the early growth stages of electricity-producing biofilms in bioelectrochemical systems (). In addition, artificially expressing cytochrome OmcS in photosynthetic cyanobacteria increased catalytic performance in a diversity of processes such as an increase in photocurrent by 9-fold (), increased nitrogen fixation by 13-fold (), and improved photosynthesis by increasing 60% biomass () compared to the wild-type cyanobacteria. These studies highlight the important role of OmcS in light-driven biocatalysis.

Therefore, improving the conductivity of OmcS could help to control bacterial physiology and ecology as well as improve biocatalysis efficiency. The increased nanowire conductivity upon cooling could also help bacteria to compensate for the loss of metabolic rate at cold environments. Furthermore, protein structure and heme geometries can significantly differ after binding to substrates () such as minerals or anodes of microbial fuel cells, potentially enhancing the coupling between the nanowires and their electron acceptors.

Our discovery also has implications in the rational design of protein nanowires (, –). Our mechanistic studies of OmcS nanowires provide a framework for rational design of nanowire structures with tunable conductivity by modulating the redox potential of hemes and heme distortions by leveraging changes in the H-bonding.

MATERIALS AND METHODS

Bacterial strains and culture conditions

G. sulfurreducens wild-type strain principal components analysis (designated DL-1) (American Type Culture Collection, 51573) was obtained from our laboratory culture collection. The cultures were maintained at 30° or at 25°C under strictly anaerobic conditions in growth medium supplemented with acetate (10 mM) as the electron donor and fumarate (40 mM) as the electron acceptor (). These strains were grown under electron acceptor limiting conditions that induces nanowire expression. Nanowires were purified via centrifugation as described previously ().

Nanoelectrode design and fabrication

Gold electrodes separated by a 300-nm nonconductive gap were designed using electron beam lithography. Silicon wafers with 300 nm of thermally grown oxide were used as substrates. The electrodes were first patterned by electron beam lithography on an electron beam resist layer and developed in a solution of cooled isopropyl alcohol. A 30-nm-thick gold film deposited on top of a 5-nm-thick titanium adhesion layer was thermally evaporated on the lithographically patterned wafer. The electron beam resist was removed using N-methyl(l-2)pyrrolidone, resulting in gold nanoelectrodes separated by a 300-nm gap, followed by sequential rinsing with acetone, methanol, and isopropanol to remove any residue. Before usage, the device was washed with distilled water and then rinsed with isopropanol to remove contaminants on the surface. The device was further plasma-cleaned for 1 min and dried with nitrogen flow, yielding a hydrophilic surface. In addition to electron beam lithography, nanoimprinting was used to fabricate electrodes. The devices were further inspected with optical and scanning electron microscopy to ensure that the gap was well separated, and resistance measurements were used to confirm that the electrodes are well insulated from each other (Rgap > 1 teraohm).

Nanowire sample preparation

Cell-free nanowire preparations were imaged with negative-staining transmission electron microscopy to ensure sample quality. Diluted 5 μl of solutions containing OmcS nanowires were placed on gold nanoelectrodes to achieve individual filaments bridging the electrodes (Fig. 2A). Nanowires on the electrodes were imaged with AFM under air-dried conditions, and height measurements were performed to confirm the presence of individual filaments. A comparison of the nanowire diameter to the atomic structure of OmcS allowed confirmation of nanowire identity. Samples were maintained under hydrated buffer environments (150 mM ethanolamine), and pH of the buffer was equilibrated using HCl. Samples were air-dried before data acquisition.

Nanowire reduction was achieved via the addition of 50-fold molar excess (versus hemes) of sodium dithionite to a nanowire suspension. Reduction was confirmed by UV-visible spectroscopy. To form films, 0.5-μl drops of reduced nanowires were deposited on interdigitated electrodes and were allowed to dry overnight in an anaerobic chamber.

Deuterium exchange

Deuterium exchange was performed by drying 10 μl of 30 μM OmcS nanowires suspended in 20 mM ethanolamine (pH 10.5) under vacuum (at least 3 hours) followed by 5 min under N2 flow. The nanowire film was resuspended in pure deuterium oxide (D2O). This process was repeated twice. For the “native OmcS” control, the same process was repeated using Milli-Q water rather than D2O.

Electrical (dc) conductivity measurements

As described previously (), all dc conductivity measurements were performed in either a two- or four-electrode configuration using semiconductor parameter analyzer (Keithley 4200A-SCS) equipped with preamplifiers, allowing 0.1-fA current resolution and 0.5-μV voltage resolution. For two-electrode conductance measurements, a fixed potential was applied across the electrodes, and the current was measured for a minimum of 120 s in sampling mode until the steady-state current was reached. For four-electrode conductance measurements, a fixed current was applied between the two outer electrodes, and the potential difference between two inner electrodes was measured over a period of minimum 120 s in sampling mode until the steady state was reached. Each current was adjusted to keep voltage difference close to 100 mV to eliminate artifacts due to hot carrier injection. The validity of the 4-probe measurements was checked by reversing the polarity of the input current. Forward and reverse currents yielded similar conductivity values, verifying the ohmic contact of the junction. The linearity of the current-voltage (I-V) characteristics was maintained by applying an appropriate low current, and the slope of I-V curve was used to determine the conductance (G = 1/R). All analysis was performed using Igor Pro software (WaveMetrics Inc.).

CP-AFM measurements

As described previously (), to measure the distance dependence of conductance using CP-AFM, 5 μl of buffer solution containing OmcS nanowires was deposited on gold electrodes patterned with electron beam lithography on a silicon wafer grown with 300-nm silicon dioxide similar to those used for four-electrode measurements. The excess buffer was absorbed with filter paper. The sample was air-dried and was mounted on a metal puck. It was further transferred to a heater-cooler sample stage inside AFM (Cypher ES, Oxford Instrument Co.), which enables temperature control over the sample.

AFM and subsequent CP-AFM experiments were performed using soft cantilevers (ASYELEC-01, Oxford Instrument Co.) with a nominal force constant of 2 N/m and resonance frequencies of 70 kHz. The tip was coated with Pt/Ir, altering the tip radius from 7 to 10 nm to 40 to 60 nm. The free-air amplitude of the tip was calibrated with the Asylum Research software, and the spring constant was captured by thermal vibration method. For CP-AFM experiments, the dual gain ORCA holder was used (Asylum Research) to record both low and high current values. The sample was imaged with a Cypher ES scanner using intermittent tapping [ac (alternating current) air topography] mode. AFM showed that gold electrodes were partially covered with nanowires to facilitate CP-AFM measurements (Fig. 3A).

After identifying nanowires across the electrodes and the substrate, the tip was withdrawn and brought into contact with the nanowire by switching on the contact mode. Points along the length of the nanowire were selected on the basis of the AFM image (Fig. 3A), and the corresponding distance from the gold electrode to the AFM tip was measured with Asylum Research Software that was used to operate AFM. With partially gold-coated substrate as the first electrode and a metal-coated AFM tip used as a second mobile electrode, the uncovered parts of the nanowire were contacted, measuring the I-V characteristics as a function of the distance between the tip and the fixed electrode by applying a bias ramp in the range −1 to +1 V (Fig. 3B), and the tip was withdrawn to the next point. To verify proper electrical contact during the measurements, the tip was frequently brought into contact with either gold or SiO2 substrate (Fig. 3B, inset). The resistance versus length curve R (L) was determined along the nanowire by fitting a straight line to the linear regime (−0.2 to 0.2 V) of the I-V data. A loading force of 10 nN was used to obtain stable current, and the same force was used for all experiments reported here. The I-V characteristics were captured in the software, and the analysis was performed in Igor Pro software (Wavematrics Inc.). For each point, at least eight repeats were recorded to demonstrate statistical significance.

The linear regime of the I-V characteristics was used for the resistance calculation. I-V curves were processed using 75-point Savitzky-Golay smoothing function before fitting using Igor Pro software (Wavematrics Inc.). The resistance (G) was calculated using the equation R = V/I. The resistance from different points was plotted as a function of distance, and the data were fit by linear regression. The y intercept revealed a contact resistance of 1.25 ± 3.6 gigaohms, which was about 10-fold lower than the nanowire resistance.

Temperature variation conductivity experiments

As described previously (), a Physical Property Measurement System (DynaCool, Quantum Design) was used to vary the temperature of a sample stage with an electrode of size 1 cm by 1.25 cm. Each sample stage included premapped pads for three sets of current and voltage channels. Electrical contacts between the sample and the specific channel input pads were made by creating wire bonds between the electrodes and pads using the operator-guided wedge-wedge wire bonder (West Bond), which uses ultrasonic energy to attach aluminum wire to contact pads. The sample-to-sample stage wiring connections were tested using either a resistance meter (Fluke) or a source meter (Keithley 2400) before inserting into the system. The dissipative power was kept under 10−6 W to eliminate self-heating effects. All conductance measurements were performed under ambient atmospheric pressure.

Raman spectroscopy

A Horiba LabRAM HR Evolution confocal Raman microscope was used to collect the Raman spectra of an OmcS film, dried on a glass coverslip. The Raman spectra were acquired with a 532-nm laser focused on the protein film with a 100× objective over a 30-s exposure. Temperature control was achieved with a thermoelectric cooling stage that we built by coupling a pair of thermoelectric Peltier elements to a copper (Cu) bar. Heat dissipation was achieved via heat sinks and fans attached to the Peltier device. The temperature of the sample was measured with a solid-state thermometer placed on the glass coverslip, adjacent to the protein film. To generate the cool sample spectrum, the sample was cooled to approximately 280 K during spectrum acquisition. The warm spectrum was acquired at an ambient temperature of 296 K. On each film, at least eight spectra were collected under both warm and cooled conditions. The spectra for each condition were then aligned with each other, smoothed with a 31-point Savitzky-Golay filter, and averaged. We measured a total of three OmcS films.

The Lorentzian decomposition of the ν10 peak between 1610 and 1670 cm−1 was performed using the Igor Pro Multipeak Fit tool. We used a linear baseline, and we constrained the fit such that all peaks should have roughly the same area. We aimed to fit six peaks corresponding to the six hemes. However, the program was only able to fit five.

UV-visible spectroscopy

We collected the UV-visible spectrum of an OmcS film deposited on an interdigitated fluorine-doped tin oxide (FTO) electrode using an Avantes UV-visible spectrometer. A blank FTO electrode was used to measure the background FTO spectrum. The spectrum of OmcS was measured under changing bias and temperature. Temperature was controlled with the thermoelectric Peltier device described above. Bias was applied using a Gamry potentiostat, whose leads were clipped to wires that were attached to the FTO electrode with silver epoxy. The estimation of charge density from UV-visible spectra is described in the “Charge mobility calculations” section.

Molecular dynamics

We performed MD simulations to obtain equilibrium sampling of heme configurations for each redox state of OmcS. The atomic model used in this study is available from the Protein Data Bank (PDB ID: 6EF8). To investigate the effect of cytochrome polymerization on the charge transfer environment, we built two models for simulation: an OmcS monomer (one chain from 6EF8) and an OmcS dimer (two consecutive chains from 6EF8). As the OmcS filament structure was solved at pH 10.5 (), the hemes were modeled in their deprotonated states, and all ionizable residues were protonated according to their standard pKa values (where Ka is the acid dissociation constant).

Each model was simulated in a hole-hopping regime (reduced hemes background) and an electron-hopping regime (oxidized hemes background). To model the state of the system at each electron (hole) transfer step, we built several models in which each heme was individually reduced (oxidized) while keeping the rest of the hemes oxidized (reduced). To avoid potential effects of heme solvation, we changed the redox state of only the central six hemes in the dimer model. The partial charges of the hemes, cysteine residues, and histidine residues were obtained from Barrozo et al. (). Each of our models were solvated with a TIP3P water model in a rectangular box with approximate dimensions of 90 Å by 90 Å by 140 Å. Sodium or chloride counter ions were added to neutralize each system.

MD simulations were performed using NAMD () with the CHARMM36 () force field parameters. The force field parameters for the hemes were obtained from (). We minimized each system followed by a 2.5-ns relaxation of the solvation water box. The models were then equilibrated to either 270 or 310 K in the constant atom count, volume and temperature (NVT) ensemble for 3.5 ns under harmonic restraints with a force constant of 0.1 kcal/mol to the amino acid side chains and a force constant of 1.0 kcal/mol to the protein backbone and the hemes. Production runs were then performed in an constant atom count, pressure and temperature (NPT) ensemble with the Langevin thermostat for 50 ns, with frames being written to the trajectory every 2.5 ps. The root mean square deviation of each system was observed to converge within the first 10 ns of the simulation. We chose 20 evenly spaced snapshots from our trajectories between 10 and 20 ns (500 ps between frames) for free energy and electronic coupling calculations. We found that 20 snapshots were sufficient for convergence of the free energy of oxidation and the electronic coupling (fig. S5). Calculation of the reorganization energy from the MD simulations required longer simulations to achieve convergence. We analyzed every frame excluding the first 10 ns.

QM/MM electronic structure calculations

To determine the free energies of oxidation and reorganization energies for each of the six hemes in OmcS, we performed electronic structure calculations using a hybrid QM/MM electrostatic embedding scheme in the Q-Chem software package (). The QM region containing a single heme and its covalently bound substituents, including axial histidine residues and cysteine residues (fig. S6), was represented with the density functional theory (DFT), while the remaining atoms were included as point charges. Covalently bound cysteine and histidine residues were truncated between the Cα and Cβ atoms. To avoid artifacts associated with the unsatisfied valency of the frontier C atom, we used the link atom approach (), adding a H capping to the Cβ. To avoid overpolarization at the QM/MM boundary, the point charges of the backbone atoms of the residues included in the QM region were set to zero. Therefore, no redistribution of charges was required to keep the net charge constant.

For the QM region, we used the ωB97X-D functional. A LANL2DZ effective core potential was used for the Fe atom, and cc-pVDZ basis sets were used for all C, N, O, S, and H atoms. The ωB97X-D functional has exact long-range exchange and dispersion correction, which eliminates self-interaction errors that may introduce artificial charge delocalization ().

We used this scheme to compute single point electronic energies for each heme from each of the relevant MD snapshots in both oxidized and reduced states. Computing the energies in both an oxidized and reduced state with a fixed nuclear configuration allowed us to obtain vertical detachment energies from which we computed the free energy of oxidation and reorganization energy within a linear response approximation (LRA), as previously described (). Values of the free energy of oxidation and reorganization energy for electron and hole transfer are reported for each heme in tables S1 and S2, respectively. The values of these parameters computed for electron transfer at 270 K are reported in table S7.

Calculation of reorganization energy from molecular mechanics

The reorganization energy as calculated in the LRA is given bywhere 〈ΔE〉RED/OX is the average of the energy gap between the reduced and oxidized states of the heme when in the RED (reduced) or OX (oxidized) equilibrium configurations. In this formalism, initially introduced by Warshel (), the energy gap serves as the reaction coordinate for the charge transfer reaction. The reorganization energy given by Eq. 1 works well when the system is ergodic (obeys the fluctuation dissipation theorem on the time scale of the reaction). However, protein systems have been shown to break ergodicity, and the simplification of the reorganization energy to Eq. 1 becomes insufficient. Matyushov and co-workers (, ) implemented a correction to account for ergodicity-breaking behavior, which uses reaction reorganization energy λr

Here, λSt is equivalent to the reorganization energy in Eq. 1, and it provides information regarding the energetic separation between the oxidized and reduced states of the heme along the reaction coordinate. λRED/OX is defined by the variance of the reaction coordinatewhich describes the width of the free energy profiles for the reduced and oxidized states.

To compute the reaction reorganization energy, we computed the energy gap (∆ERED/OX) using the NAMD pair interaction command. An estimation of ∆E is given bywhere j is the atom index, ∆q is the difference in charge between the oxidized and reduced states of atom j, and ϕEW is the Ewald lattice sum electrostatic potential of the thermal bath (surrounding protein and water). This is implemented by computing a pair interaction calculation on each frame of our MD trajectory, which specifies the heme of interest to have partial charges equal to the difference in charge between the oxidized and reduced states. The electrostatic energy output in the NAMD log file is the energy gap ∆E.

We extended the MD simulations until λr converged. We considered our calculations to be converged when the value of λr fluctuated by less than thermal energy in the last 10 ns of the simulation. The convergence of λr is plotted in fig. S5.

Provided that there has been recent debate regarding the proper way to compute reorganization energies, we have included results from both the LRA methodology and the MM methodology with the reaction reorganization energy correction for nonergodic dynamics. Matyushov and co-workers have presented a ratio of the reaction reorganization energy and Stokes shift reorganization energy for use in determining whether a reaction exhibits nonergodicity: κ = λ/λSt. Typically, a value greater than 1 is indicative of nonergodicity. Apart from heme 3 for hole transfer, the charge transfer half reactions for all hemes have a value of κ greater than 1 (tables S10 to S12). However, all values fall in the range [0.8, 1.34]. These values are lower than the value reported for cytochrome c (κ = 2.26) ().

Rather than picking one methodology over the other to be the correct approach, we have elected to present our results as an objective comparison to experiment. The conductivity that we computed using the MM-derived reorganization energies was higher than the conductivity computed with LRA reorganization energies and thus was in better agreement with experiment. Hence, the comparison of computed conductivities at 310 and 270 K is performed with just the MM-derived reorganization energies.

Electronic coupling calculations by fragment orbital DFT

The electronic couplings for electron transfer and hole transfer were calculated using the ADF () package, which implements the charge transfer integral (CTI) method (). The calculations were performed at a PBE0/DZ level of theory. This combination of functional and basis set consistently produced stable wave functions whose frontier orbitals contained the d-orbitals of the heme iron.

Using the fragment orbital approach, we split each pair of hemes into two fragments. Each fragment was equal in size to the QM region that we defined for the QM/MM electronic structure calculations (fig. S6). As in our electronic structure calculations, the coupling calculations were performed in the presence of external point charges corresponding to the surrounding protein and solvent. Prior work has suggested that the protein environment and solvent have a negligible effect on electronic coupling values (); however, we found that their presence greatly influenced the energetic ordering of the molecular orbitals and the stability of the occupied orbitals. Initially, we ran a single-point calculation on each fragment, followed by the CTI calculation to get the coupling. In brief, the CTI method defines the coupling as follows ()

Here, JAD is the off-diagonal element of the Fock matrix corresponding to the donor and acceptor orbitals (CTI); SAD is the overlap integral of the donor and acceptor orbitals; and eD and eA are the energies of the sites occupied by the charge on the donor and acceptor orbitals, respectively. We used this equation to compute the coupling between all pairs of frontier orbitals [highest occupied molecular orbital (HOMO) − 2 to HOMO] containing substantial d or d character. These are denoted as the dπ orbitals. The total coupling between a pair of hemes was calculated as a combination of the individual couplings between dπ-containing orbitals as described previously ()where denotes the electronic coupling between the ith and jth dπ-containing molecular orbitals of the donor (D) and acceptor (A) fragments. ccorr is a correction factor that corrects for the lack of polarization in fragment orbital DFT (FODFT) calculations. Prior work has calculated ccorr by taking the ratio of the coupling values obtained using constrained DFT (CDFT) and FODFT. Previous studies have used values of 1.75 () and 1.38 (). We used CDFT as implemented in Q-Chem to compute coupling for a subset of our MD snapshots, and we obtained a value of 1.70 for ccorr. This value was used to scale our FODFT coupling values. Note that we computed coupling using CDFT for only a subset of MD snapshots as many of the snapshots did not converge. FODFT allowed for better sampling and thus proved to be the most feasible technique. For calculations of electron transfer couplings, snapshots were selected from the fully oxidized MD trajectory, and for hole transfer couplings, snapshots were obtained from the fully reduced MD trajectory.

Calculation of electron transfer rates

Electron transfer and hole transfer rates were calculated in the nonadiabatic limit using the theory developed by Marcus ()where ΔG0 is the free energy for charge transfer from heme i to heme j, λ is the reorganization energy for charge transfer from heme i to heme j, and HAD is the electronic coupling between hemes i and j. The free energies and reorganization energies used in this calculation are computed from the free energies of oxidation and reorganization energies previously computed via the LRA as follows (for self-exchange reaction)

These values and the electronic coupling values are reported in tables S3, S4, and S8. The computed charge transfer rates are reported in tables S5, S6, and S9.

Charge mobility calculations

We implemented a single-particle kinetic Monte Carlo simulation, adapted from the scheme used by Byun et al. () to compute steady-state flux in the decaheme cytochrome MtrF. We initialized a reduced state on heme 1 of a six-site model with periodic boundary conditions. The system was allowed to evolve for 100,000 iterations with the electron changing sites according to the conditionwhere m is an index of the list rates corresponding to the possible events, ξ1 is a uniformly distributed number on the range (0, 1), and ktotal is defined as the sum of rates for all of the possible events at the given point in time. As there is only one electron in the system, ktotal is the sum of k and k, the forward rate and reverse rate for electron transfer to the hemes adjacent to site i. In addition, m is always equal to 2. The effective time for each iteration was determined to bewhere ξ2 is another uniformly distributed random number on the range (0, 1).

The mobility was determined by recording the mean squared displacement of the electron as a function of the simulation time. The distance per step was determined by the average iron-to-iron distances obtained from our MD simulation. The mean squared displacement was fit to the expressionto extract the diffusion coefficient (D). This value was divided by thermal energy–giving mobility.

Normal coordinate structural decomposition

In our implementation of the normal coordinate structural decomposition (NSD) analysis, we used an unsubstituted planar Cu(II) porphine with D4h symmetry as the reference structure. In the first implementation of NSD, a Cu(II) porphine was used as the Cu(II) metal minimizes internal strains, which lead to in-plane and out-of-plane distortions of the optimized structure (). We optimized the geometry of our reference structure and computed its normal modes with DFT using Gaussian 16. We used a B3LYP functional, a 6-31G(d) basis set for the H, C, and N atoms, and an LANL2DZ basis set with an effective core potential for the Cu atom.

We used a custom Tcl script, executed with VMD (), to compute a mass-weighted difference vector (xdiff = xobs − xref), which describes the displacement of each Cartesian coordinate of each atom in the heme macrocycle, including the metal, from the optimized reference structure. We then used the same script to compute the absolute value of the scalar product of the difference vector with the mass weighted unit vectors (q) describing the normal modes of the reference structure (xdiff • q). We restricted our analysis to the normal modes belonging to the following classes of symmetry: B2u, B1u, A2u, Eg, and A1u. These symmetry groups correspond to the saddling, ruffling, doming, waving, and propellering distortions, respectively. To benchmark the ability of our program to capture these distortions, we compared the results of NSD analysis of the heme in horse heart cytochrome c (PDB ID: 1HRC) using our program to the results given by the publicly available NSD program described by Graves et al. (). Both programs gave similar results (fig. S7), demonstrating that our program can capture the relative contributions of each of the porphyrin distortions. To determine the total out-of-plane displacement, we computed the root mean square of a modified mass-weighted difference vector composed only of the displacements that lie along the unit normal vector () of the porphyrin plane. This was achieved by multiplying the x, y, and z components of with the corresponding entries in xdiff for each atom.

Calculation of electrostatic potential and electric field

We calculated the protein-induced electrostatic potential and electric field using the Delphi program. The point charges of the heme of interest were set to zero, and a site calculation was used to compute the potential and electric field components at each atom site. The reported electrostatic potential was computed as the average potential over all heme atoms. The reported electric field magnitude was taken from the site of the heme iron atom.

Calculation of H-bonding structure

The H-bonding structure of OmcS was calculated with the VMD () HBonds plugin, applied to the fully oxidized OmcS MD simulation. The instantaneous H-bonding count was extracted from the default output file, while the CHF was calculated from the detailed output file. We used a custom Python script to generate a H-bonding matrix whose entries A(i, j) described the frequency at which a H-bond was observed between residues i and j in our MD simulations. The CHF was calculated as the Frobenius norm of the matrix. The difference CHF was calculated as the Frobenius norm of the difference between the matrices corresponding to 310 and 270 K.