Electron-transfer chain in respiratory complex I

Complex I is a part of the respiration energy chain converting the redox energy into the cross-membrane proton gradient. The electron-transfer chain of iron-sulfur cofactors within the water-soluble peripheral part of the complex is responsible for the delivery of electrons to the proton pumping subunit. The protein is porous to water penetration and the hydration level of the cofactors changes when the electron is transferred along the chain. High reaction barriers and trapping of the electrons at the iron-sulfur cofactors are prevented by the combination of intense electrostatic noise produced by the protein-water interface with the high density of quantum states in the iron-sulfur clusters caused by spin interactions between paramagnetic iron atoms. The combination of these factors substantially lowers the activation barrier for electron transfer compared to the prediction of the Marcus theory, bringing the rate to the experimentally established range. The unique role of iron-sulfur clusters as electron-transfer cofactors is in merging protein-water fluctuations with quantum-state multiplicity to allow low activation barriers and robust operation. Water plays a vital role in electron transport energetics by electrowetting the cofactors in the chain upon arrival of the electron. A general property of a protein is to violate the fluctuation-dissipation relation through nonergodic sampling of its landscape. High functional efficiency of redox enzymes is a direct consequence of nonergodicity.


A. System setup
The initial coordinates used for all simulations of the hydrophilic peripheral arm of complex I (NADH:ubiquinone oxidoreductase) were adopted from the Protein Data Bank entry 3IAM and correspond to bacterium T. thermophilus. 1 The original PDB file contains eight different protein subunits, a 1,4-dihydronicotinamide adenine dinucleotide (NAD) molecule, a flavin mononucleotide (FMN), and nine ironsulfur clusters (ISCs). These cofactors include two Fe 2 S 2 (N1a, N1b) and seven Fe 4 S 4 (N2, N3, N4, N5, N6a, N6b, N7) ISCs. The topological structure was created by separating all subunits, individual molecules, and cofactors and recombining with the addition of hydrogen atoms and application of all force-field parameters using the VMD subprogram psfgen. 2 The force field topology and parameters for the protein subunits were taken directly from the CHARMM36 (with CMAP torsion corrections) force field. 3 CHARMM force field parameters for neutral FMN were taken from studies of the flavin-binding LOV (Light-Oxygen-Voltage) domains in the phototropin receptor protein. 4,5 CHARMM parameters for both NADH and NAD + have previously been published and the parameters corresponding to NAD + were applied here. The parameters for the iron-sulfur clusters (including modified ligated protein amino acids) were taken directly from DFT calculations performed by Chang and Kim. 8 The electron transfer reactions between the ISC cofactors follow the hydride transfer from NADH to FMN. We therefore considered the protonated form of FMN (FMNH) and deprotonated form of NADH (NAD + ). The new charges for the FMNH were calculated with DFT (B3LYP/6-311+G ) implemented in Gaussian09. A single hydrogen atom was added to the N5 nitrogen of the flavin and the ribitol alcohol was removed and replaced by a single hydrogen atom (closed shell). Following the geometry optimization, the partial atomic charges were obtained from a fit to electrostatic potential (ESP) using the "CHelp" option. 9 Once the charges were obtained, the ribitol alcohol was added back to the structure and the residual charge from the extra hydrogen was added back to the corresponding carbon atom to assure unit charge on both the flavin and ribitol substructures.
The ISC cofactors were attached to the protein matrix through applied patches using psfgen for all corresponding ligations, with the patches obtained from the supplemental material provided by Chang and Kim. 8 These patches include both oxidized and reduced charge states for all ISCs plus ligands, along with bond, angle, and dihedral force field parameters.
Finally, after all the charges and topological information (bonds, angles, dihedrals, improper angles, etc) were applied to complex I, TIP3P water was added to the system using the CHARMM force field and the solvate plugin from VMD. Ions were added to the system using the autoionize plugin, with the resulting ionic concentration of ≃50 mM/L. The resulting system consisted of 481641 atoms in an orthorhombic box with x, y, z box lengths of 178.9 × 148.3 × 178.1Å (Fig. S1). The box dimensions were taken after a short 20 ns NPT simulation (details below) and held constant throughout all subsequent NVT simulations.

B. Simulation protocol
All simulations were performed using NAMD 2.9. 10 The initial structure described in the previous section was first simulated (NPT) for 20 ns in order to relax the waters surrounding the complex I with all ISCs in their oxidized states except for the reduced N4. This was followed by 5000 steps of steepest decent minimization. This initial equilibration was followed by 340 ns of production for the charge transfer state N4 − /N5/N6a, with the electron on N4, using the NVT ensemble. The starting configuration for the charge state N4/N5 − /N6a was extracted at the 60th ns of the previous run. The new charges were then applied and a 328 ns NVT simulation was carried out. Finally, the configuration after All production simulations were done using the NVT ensemble with the Langevin temperature control set to 300 K and a damping coefficient of 5 ps −1 . Long-range electrostatic interactions were treated with the particle mesh Ewald technique using a cutoff distance of 12.0Å. A 2.0 fs timestep was used for all simulations, along with a 1 ps saving frequency, except for the initial 60 ns of the charge state N4 − /N5/N6a for which the saving frequency was 2 ps. It was noticed that after approximately 50 ns, some of the ISCs would collapse resulting in unphysical geometries due to a new minimum in the angle and dihedral energies that was not offset by the van der Waals interactions. For this reason, we increased the problematic angle and dihedral force constants (making the structures slightly more rigid) incrementally making it impossible for the ISCs to populate these unwanted and unphysical geometries and assuring the stability of the simulations.
The center-to-center distances between the cofactors obtained from MD simulations are listed in Table S1. The protein matrix does not show major conformations in response to moving electrons. However, the hydration levels of N5 and N6a cofactors change upon the arrival of the electron (see below). In addition, there is a shortening of the distance between the cofactors in the reduced N6a − state.

C. Quantum calculations
The coordinates of a particular iron-sulfur cluster were extracted from the PDB file (PDB:3IAM 1 ) and used as the starting point for the quantum mechanical (QM) calculations. All QM calculations were done with Gaussian 09. 9 The original cysteine (and histidine) residues were removed from the entire ligated irons-sulfur clusters, keeping only the organic sulfur (or nitrogen) atoms and ligated carbon atoms. Hydrogen atoms were then added to fill the valencies. The hydrogen atoms were geometrically optimized keeping all heavy atoms frozen.
To model the antiferromagnetic coupling between the iron atoms within the ISCs, we used the broken symmetry (BS) 11 method by first employing the fragment sub-directive under the guess directive within Gaussian. This directive allows one to specify the α and β spins separately for unpaired electrons in each fragment. All Within the protein matrix, the ligands R/R' are cysteine residues for the iron-sulfur clusters N4 and N6a, while for N5, the ligand R' is a histidine. For the quantum calculations performed here, R/R' were SCH3 in place of the cysteine residues for N4 and N6a and R' was NH3 in place of the histidine residue in N5. Due to the slight asymmetry adopted in the geometry from the original PDB file, quantum mechanical calculations were performed using initial guess Hartree-Fock wavefunctions from all allowed permutations of the charge/spin configurations.
ground and excited state calculations were performed using the BS ZINDO technique. 12,13 In this methodology, the initial guess wave function is found by matching the basis set used in ZINDO with the guess and fragment keywords in Gaussian and then proceeding to the SCF and the excited state calculations.
Since the redox properties of the iron-sulfur clusters are dependent on both the Fe and S atoms, we kept the Fe and inorganic S atoms as separate fragments in the initially created guess Hartree-Fock wavefunction. Therefore, a total of twelve separate fragments were used in the QM calculations; the four iron atoms, the four inorganic sulfur atoms, and the four remaining ligands (SCH 3 /NH 3 ). For all QM calculations, the inorganic sulfur atoms were taken in the state with a formal charge of −2 and with a spin multiplicity of 1, while the charge and spin multiplicity of SCH 3 (NH 3 ) ligands were taken to be −1 and 1 (0 and 1), respectively.
The transition metal iron atom is typically found to occur in one of three possible oxidations states (+2, +3, or +4) with various possible spin states (Fig. S2). The iron atoms within the tetranuclear iron sulfur clusters in complex I are found to exist under physiological conditions in two redox states with formal iron charges (+3,+3,+2,+2) and (+3,+2,+2,+2). 14 The oxidized (Ox) state consists of two ferric (+3) and two ferrous (+2) atoms anti-S3 ferromagnetically coupled and forming a complex with total spin S = 0. The reduced state has a single ferric iron atom and three ferrous atoms with an effective spin S = 1/2. The S = 7/2 reduced state has a comparable low lying ground state energy. 15 This S = 7/2 reduced state configuration was also included in the calculation even though it has been shown that the S = 1/2 state is slightly lower in energy. 15 It was also shown that the four iron atoms within tetranuclear clusters were indistinguishable via the Mössbauer spectroscopy. 16 The polarizabilities for the reduced and oxidized states of the iron-clusters can be found by using the excited state transition energies and the transition dipole moments all of which are provided within the ZINDO framework of Gaussian'09 and the BS-ZINDO methodology used here. 12,13 Explicitly, the polarizability tensor for a given structure can be calculated using where ∆E 0j = E j − E 0 is the difference in energy between excited state j and the ground state and β, γ represent Cartesian coordinates. The scalar polarizibilities reported in Fig. S3 are traces of the corresponding tensors α 0 = 1 3 Tr [α 0 ]. The polarizabilities shown in Fig. S3 correspond to SCFs produced for the N5 cofactor for different initial configurations. The red points within the figure are polarizabilities calculated for the dominant spin state along the classical trajectory, which is determined from to the lowest energy eigenvalue of the complete Hamiltonian.
For comparison, we calculated the polarizabilities for the cofactors using the double and triple zeta basis sets (cc-pDVZ and cc-pTVZ) which include polarization functions and the full configuration interaction singles (CIS) method included in Gaussian along with the broken symmetry framework described in the main text under the Methods section (designated here as BS-CIS). One can see in Fig. S4 that the dominant spin state contribution (shown by the red open points) is very slowly converging. Because of the computational overhead (large amounts of memory and very long computational times) needed for the full CIS calculations, slowly converging polarizabilities, and the quality of the results provided by BS-ZINDO for a fraction of the computational cost, the BS-ZINDO calculation was used for all main results in this study. The full energy gap calculation was also carried out using this methodology including 50 excited states in contrast to the included 100 excited states for the BS-ZINDO method. The energy gap statistics (Stokes shift and variational reorganization energies) for the BS-CIS including 50 excited states is given in Table S3 and are consistent with the BS-ZINDO calculations.

D. Water Penetration
To understand water's role in solvating the cofactors, we calculated the average number of water molecules surrounding the N6a iron-sulfur cluster along the trajectory ( Fig. S5). One can see that no waters are present within a 3Å shell when the electron is on N4 and begin to appear after the electron arrives to N5. In the final state, with the electron on N6a, the average number of waters increases on the time-scale of ∼ 150 ns from nearly 0 to ∼ 6 waters. This increase in hydration corresponds to an increase in fluctuations of the electrostatic potential, as well as an increase in the fluctuations of the electric field.

E. Correlations between the cofactors
Diagonalizing the Hamiltonian matrix along the trajectory produces the lowest-energy eigensatate, with the energy E a g (D/A), among all spin configurations in either Red or Ox state of each cofactor (a = Red, Ox). These four energies allow us to define the instantaneous energy gap 17 sampled along the trajectory. In Table S2, we have separated X into the donor and acceptor contributions in the variance (δX) 2 . This separation allows us to look at individual contributions to the total variance and at the cross term. Except for the N4 − → N5 transition, the total variance of X is largely dominated by the sum of the variances of individual cofactors.  S2)) produced from the QM/MM formalism. The total is the sum of the donor and acceptor components, . All energies are in eV. The mean of the values in the last column for a given reaction yields the reorganization energy λ in Table I in the main text. The QC subscript in classical calculation using the electrostatic potential (X = ∆φ QC ) indicates that the potential was averaged over the quantum center, the cofactor in this case, as well as along the simulation trajectory.

II. Free energy of electron transfer
Here we provide the derivation of Eqs. (8) and (9) in the main text, connecting the free energy of electron transfer to the free energy of Coulomb interaction between the cofactors in the chain. We start with a general problem of electron transfer between two cofactors carrying the charges Z Ox D (at the donor) and Z Ox A (at the acceptor) in their oxidized (Ox) states. The total charge of the donor is therefore Z Ox D − 1 in the initial state and the total charge of the acceptor is Z Ox A − 1 in the final state (Fig. S6). The free energy of electron transfer can be written in the linear response approximation as 18 Here, ∆E gas is the gas-phase energy gap and χ L is the longitudinal response function of the medium, which generally can be given as a non-local kernel. 19 The electric field of the donor-acceptor complex E i in the initial (i = 1) and final (i = 2) states is integrated over the volume occupied by the thermal bath (protein and water). The asterisks notation involves both the volume integral and the scalar product The fields in the initial and final states can be represented as sums of the corresponding donor and acceptor fields where E Ox,Red are the fields of the donor (D) and acceptor (A) in the corresponding redox states. After some manipulations, we obtain where E e D,A is the electric field of the electron localized either on the donor or the acceptor. Further, ∆G D,A s are the solvation free energy differences between the Red and Ox states. By using the Coulomb law, this equation can be transformed to where R DA is the donor-acceptor distance. One has additionally to keep in mind that the free energy of the electrostatic field of the donor acceptor complex is contained in ∆E gas Combining these terms together one gets where ∆E gas is now the difference of the gas-phase energies of the isolated donor and acceptor fragments, excluding their Coulomb interaction. The standard assignment of the longitudinal susceptibility is 19  1 − ǫ −1 eff , where ǫ eff is the effective dielectric constant of the medium. One obtains the final expression (S10) Expressing the first three terms in this equation in terms of the standard (midpoint) reduction potentials of the donor and acceptor, one arrives at Eq. (7) in the main text. For charge separation reactions, Z Ox D = 1 and Z Ox A = 0 and one arrives from Eq. (S10) at the standard Rehm-Weller result used in the literature. 20 The effective dielectric constant ǫ eff in Eq. (S10) is not known and this relation is not practical for most heterogenous systems. An alternative can be thought in finding the interaction energy between the charges directly from simulations.
The free energy of two charges placed in a polarizable medium can be written as the sum of the direct Coulomb interaction between them and the variance of the medium potentials φ 1 and φ 2 at the location of the charges (linear response approximation 21 ) The self variances in this equation (δφ i ) 2 produce solvation free energies and the cross variance is responsible for the potential of mean force, i.e., the alteration of the free energy of bringing the charges from infinity to the distance R ∆G s (R) = q 1 q 2 R −1 − β δφ 1 δφ 2 . (S12) The results of calculations of cross-variances of the medium electrostatic potentials at the cofactors of the charge-transfer chain are listed in Table S4. These results are combined with distances listed in Table S1 to calculate the profile of reaction free energies along the chain shown in Fig. 5 in the main text.