Horizontal proton transfer across the antiporter-like subunits in mitochondrial respiratory complex I

Abstract

Respiratory complex I is a redox-driven proton pump contributing to about 40% of total proton motive force required for mitochondrial ATP generation. Recent high-resolution cryo-EM structural data revealed the positions of several water molecules in the membrane domain of the large enzyme complex. However, it remains unclear how protons flow in the membrane-bound antiporter-like subunits of complex I. Here, we performed multiscale computer simulations on high-resolution structural data to model explicit proton transfer processes in the ND2 subunit of complex I. Our results show protons can travel the entire width of antiporter-like subunits, including at the subunit–subunit interface, parallel to the membrane. We identify a previously unrecognized role of conserved tyrosine residues in catalyzing horizontal proton transfer, and that long-range electrostatic effects assist in reducing energetic barriers of proton transfer dynamics. Results from our simulations warrant a revision in several prevailing proton pumping models of respiratory complex I.

---

Introduction

Mitochondrial respiratory complex I (NADH:ubiquinone oxidoreductase) is a large (∼1 MDa) enzyme in the electron transport chain.^1,2 It couples the processes of NADH oxidation and quinone reduction to the translocation of four protons across the inner mitochondrial membrane. The generated protonmotive force is utilized by the respiratory chain for ATP synthesis. Despite a plethora of high-resolution structures of complex I,^3–10 its proton pumping mechanism still remains a mystery.

Mitochondrial complex I comprises 14 core subunits, organized into a characteristic L-shape: a peripheral arm where redox reactions take place, and a membrane part working as a proton pump (Fig. 1A). The core subunits are shared between bacterial and mitochondrial enzymes. However, the latter have several supernumerary (accessory) subunits, the function of which largely remains unknown. In the core membrane domain of complex I, the three antiporter-like subunits (ND2, ND4 and ND5, Fig. 1A) are responsible for transferring protons from the N side to the P side of the membrane.^1,2 Interestingly, these three antiporter-like subunits are interlinked to four other membrane-bound core subunits (ND6, ND4L, ND3 and ND1) by a central hydrophilic axis (Fig. 1A). This ‘river’ of charged and polar residues together with water molecules provide connectivity between the quinone binding tunnel and the terminal end of membrane arm, spanning ca. 200 Å (Fig. 1A). The putative proton transfer pathways in antiporter-like subunits have been investigated by structural biology techniques, site-directed mutagenesis and computational approaches.^3–16 It has been proposed that each of the three antiporter-like subunits pick protons from the N side and release to the P side as a response to redox activity in the quinone binding tunnel. Proton transfer from the N side is assumed to first occur to a conserved Lys residue in transmembrane helix (TMH) 8 (Lys241 in ND2), located in the central part of the subunit (all amino acid numbering corresponds to complex I from Yarrowia lipolytica, unless otherwise stated). Subsequent proton transfer from this to another conserved lysine in ND2 (Lys383 or a glutamate in ND4 but with function shared¹³), and putatively to the P side is driven by the dynamics and/or protonation reactions of conserved Lys/Glu pair from TMH 7/5.^1,15,16 However, how this is achieved at a molecular level remains elusive and debated.^6,15

Fig. 1 Overall architecture of respiratory complex I and its membrane bound antiporter-like subunit (A) structure of mitochondrial complex I from Y. lipolytica (PDB 7O71). The core subunits of the membrane arm (ND1, ND2, ND3, ND4, ND5, ND4L, ND6) are shown in the ribbons representation, whereas rest of the protein as grey surface. Structurally-resolved water molecules (red spheres) and conserved residues (cyan) in the central hydrophilic axis, which connects the Q tunnel (brown surface) with terminal end of membrane domain, are shown. Protein-bound electron transfer components FeS clusters (yellow-orange spheres) and FMN molecule (licorice) are also displayed. Proton and electron transfer are marked in blue and red arrows, respectively. (B) Colored dotted boxes highlight the QM regions (as part of QM/MM free energy simulations) studied in the work. Purple (L region), green (R region), light blue (ND2–ND4L interface) and light green (ND4–ND2 interface). The “big QM region” comprises the protein residues and water molecules from both R and L regions together. All amino acid residues that are part of each of the QM regions are displayed in licorice, and structural waters as red spheres. Water molecules that are not part of QM, but MM, are shown in dotted representation (red).

Recently, proton pumping models in which antiporter-like subunits act as self-containing proton pumps have been challenged by the high-resolution structural data on complex I.⁶ In contrast to a clear hydration observed towards the P side exit in terminal ND5 subunit,^10,14 analogous path in ND2/4 subunits is blocked by bulky hydrophobic residues.¹⁰ This has led to the proposal that only the terminal antiporter-like subunit ND5 is responsible for the release of protons to the P side,^6,10 whereas protons are fed into the central hydrophilic axis from the N side via pathways present in all three antiporter-like subunits.¹⁰

A recent high-resolution (2.4 Å) cryo-EM structure of Y. lipolytica complex I revealed extensive hydration in the middle of the ND2 subunit (Fig. 1B). Refinement with classical molecular dynamics (MD) simulations showed a fully hydrogen bonded connectivity across the antiporter-like subunit.¹⁰ This led to the suggestion that protons can move on the hydrogen-bonded proton wire, and can travel long distances in the membrane domain of complex I. Here, based on classical and quantum mechanical simulations, we show protons can indeed migrate through the hydrated path of antiporter-like subunits and crucially also at the subunit–subunit interface. To track the energetic feasibility of proton diffusion across the entire width of antiporter-like subunits, we performed quantum mechanical/molecular mechanical (QM/MM) free energy calculations, which reveal low energy barriers of proton transfer dynamics in the membrane domain. Our results suggest long range proton transport occurs in the membrane-bound subunits of complex I and justifies the revision of models of proton pumping in the membrane domain of respiratory complex I.

Results

Proton transfer pathways in the ND2 subunit studied by QM/MM MD simulations

We applied QM/MM free energy simulations on the high-resolution complex I structure from Y. lipolytica.¹⁰ We restricted the size of our QM/MM model system to 6 core subunits: ND1, ND2, ND3, ND4, ND4L, and ND6. The hydrated central axis of ND2 subunit was treated with full QM description and the rest of the system was described by MM framework (see Methods). To keep the computationally demanding free energy simulations tractable, we split the ND2 QM region into two parts (see Fig. 1B; S1A and B, ESI†). The first one (R) comprises the region between the highly-conserved Glu131 of conserved Lys/Glu pair from TMH7/5 and the central Lys241 (proximal to Lys/Glu pair, see Fig. 1). The second region (L region) spans the section between Lys241 and terminal Lys383, which is distal to the Lys/Glu pair and closer to ND2/ND4 interface (Fig. 1). These four highly conserved titratable residues (Lys211 and Glu131 of Lys/Glu pair, central Lys241 and terminal Lys383) are all known to be important for the proton pumping activity of complex I.^17,18

Based on pK_a calculations¹⁹ and classical MD simulations, we first modelled these residues in their standard (charged) protonation states and performed unbiased QM/MM MD (molecular dynamics) simulations. We did not observe any spontaneous proton conductance in either of the two ND2 regions (R or L), or in a QM/MM MD simulation of a much larger QM region (Fig. 1B, see methods, Tables S1–S4†). We thus created proton vacancies to drive the forward proton transfer via a Grotthuss “hole” mechanism.²⁰ A proton vacancy at the end of the “L” region was created by neutralizing terminal Lys383. pK_a calculations performed on snapshots obtained from classical simulations demonstrate low proton affinity of Lys383 (Fig. S2†), thereby justifying its modeling as a proton hole. A similar low proton affinity of a buried lysine residue has also been observed in another protein system.²¹ A ∼5.5 ps QM/MM MD trajectory did not result in any spontaneous proton transfer to neutral Lys383 from central Lys241 (protonated) suggesting the presence of energy barriers that could not be overcome in unbiased simulations. However, an excess proton (modelled as H₃O⁺ ion) at several structurally-conserved water positions (Fig. S1A†), resulted in rapid proton transfer to terminal Lys383, within 50–200 fs of the unbiased QM/MM MD simulation runtime. Remarkably, this involved the deprotonation (and re-protonation) of highly conserved tyrosine (Tyr413 and Tyr316) and Ser288 residues in the ND2 subunit.

With the indication that deprotonation of these residues may be responsible for the barrier in proton transfer dynamics, we initiated umbrella sampling (US) simulations to induce proton transfer to neutral Lys383 from the protonated Lys241. The potential of mean force (PMF) profile derived from this setup shows an energy barrier (∼6.5 kcal mol⁻¹, Fig. 2, note the statistical error bars displayed in the figure) with a favorable proton position on Lys383. The obtained activation energy barrier is similar to the one obtained previously based on calculations on a lower resolution structure of mammalian complex I.¹⁶ The additional analysis of the hydrogen bonding patterns along the pathway suggests the preservation of the distinct hydrogen bond network in initial, final, and intermediate stages of proton transfer (ESI, Fig. S3 and S4†). Interestingly, in place of an asparagine (Asn292) in Y. lipolytica complex I, mammalian complex I consists of a titratable residue His186 (Mus musculus complex I numbering, ESI, Fig. S5†). Given that the energetics obtained in two studies are comparable, we suggest a titratable residue at this location is not necessary, instead, proton transfer would follow the path discussed here involving protonation/deprotonation of highly conserved Tyr413 (Fig. 2 and S1A, ESI†). Indeed, proton transfer energetics obtained on the path from Tyr413 to Lys383 revealed a similar barrier of 6–7 kcal mol⁻¹ (ESI, Fig. S6†), where spontaneous re-protonation of Tyr413 occurred from protonated Lys241 with Tyr316 participating in the transfer (ESI, Fig. S1A and Movie S1†). We further analyzed the simulated trajectory and noted that the two energy maxima displayed in Fig. 2 correspond to the deprotonation of tyrosines, which are transiently stabilized by positively charged lysines (transient dipoles, see Movie S2†). Moreover, the proton transfer via tyrosines and serine occurs in a semi-concerted fashion, as also observed in other enzymes.^22,23 These data not only allow us to identify the cause of energetic barrier in this region, but also that the proposed proton transfer route is a likely one.

Fig. 2 Proton transfer from central to terminal lysine (L region). Proton transfer occurs from Lys241 (protonated) → H₂O → Tyr316 (neutral) → H₂O → Ser288 (neutral) → Tyr413 (neutral) → H₂O → H₂O → Lys383 (neutral) as depicted in the top panel. The bottom panel shows PMF (potential of mean force) profile of proton transfer derived from our QM/MM umbrella sampling simulations. Purple shaded area around the bold purple line depicts the bootstrapping error range (see also methods). Occupational histograms for the reaction coordinates are shown in ESI, Fig. S17A.† The inset depicts the core membrane subunits of respiratory complex I: ND4 (dark blue), ND2 (pink), and ND4L (green). The purple circle denotes the QM region where the proton transfer is investigated.

Overall, our QM/MM MD simulations on high-resolution structure of Y. lipolytica complex I suggest that proton transfer from central Lys241 to terminal Lys383 occurs with reasonable energetics and involves protonation/deprotonation reactions of conserved tyrosine residues (Tyr316 and Tyr413), which are ideal candidates for site-directed mutagenesis in bacterial complexes.

Horizontal proton transfer from conserved Lys/Glu pair to the conserved central lysine

The R region (Fig. 1B) has been suggested to catalyze proton transfer,¹⁰ but no charge transfer dynamics have been studied. Therefore, it remains unknown if proton migration on this path is energetically feasible or not. Explicit proton transfer in this section of the ND2 subunit is critical for those models of proton pumping which propose proton migration along the entire length of membrane arm of complex I.^6,10 Therefore, we launched equilibrium unbiased QM/MM MD simulations of the R region with Lys241 modeled neutral (Tables S1 and S3†). Similar to the L region, we saw no spontaneous proton transfer from conserved Lys/Glu pair to the created hole position at Lys241, highlighting the presence of high energy barriers.

Based on the critical role of tyrosine residues found above in the L region, we studied the proton transfer in R region in short-to-long fragments using QM/MM free energy simulations (Table S3† and setups 1–3,9). Proton transfer on the pathway connecting conserved Tyr225 to central Lys241 occurred with an activation energy barrier of ∼4.0 kcal mol⁻¹ with equally favorable proton positions on both residues (Fig. 3A). However, despite protonation of Lys241 from Tyr225, subsequent propagation of the “hole” in the direction of Lys211 did not occur, which had now become “stuck” on Tyr225, suggesting the existence of activation energy barrier in the region from Tyr225 towards Lys/Glu pair (Fig. 1B and S1B, ESI†).

Fig. 3 Free energy profiles of proton transfer from conserved Lys/Glu pair to central lysine (R region). (A) Tyr225 (neutral) → H₂O → Lys241 (neutral) pathway. (B) Lys211 (protonated) → H₂O → H₂O →Tyr157 (deprotonated) pathway. (C) Proton transfer from Tyr157 (neutral) through H₂O to Tyr225 (deprotonated). The left panel shows energetics when Glu131 of conserved Lys/Glu pair is deprotonated, and when protonated (right panel). The bold green lines are PMF profiles with light green shaded area highlighting the bootstrapping error range (see also methods). Occupational histograms for the reaction coordinates are given in ESI, Fig. S17B–E.† Inset shows ND4, ND2, and ND4L subunits in blue, pink, and green surface representation, respectively. The QM region where the proton transfer is studied is marked with the dark green circle.

To investigate what will drive proton transfer from Lys211 of Lys/Glu pair, we initiated additional free energy simulations on the structurally-conserved Tyr157-H₂O-Tyr225 pathway (Fig. 1B and S1B, ESI†), where Tyr225 was modelled anionic, and Tyr157 neutral (Table S3,† setup 2). The energy barrier we obtained was relatively higher (∼7.5 kcal mol⁻¹) and, at the same time the protonation of Tyr225 (anionic) found to be highly unfavorable (Fig. 3C, left), in agreement with our above conjecture. We noted that in all these calculations, Glu131 of the conserved Lys/Glu pair was modelled deprotonated (ESI, Table S3†). However, surprisingly, both the kinetics and thermodynamics of proton transfer between the tyrosine residues substantially improved when Glu131 was modelled protonated (Fig. 3C, right). This is a strong indication that long range electrostatic effects are at play and proton transfer in this domain can be triggered by protonation of Glu131 of Lys/Glu pair, as long as the proton hole resides on Lys241 (Fig. 3C, see also ESI, Fig. S7†). pK_a calculations performed on simulation snapshots obtained from classical molecular dynamics indeed support a higher proton affinity of Glu131 (see Fig. S2† panel D). We additionally tested this notion by modelling an H₃O⁺ ion (excess proton) between protonated Lys211 and neutral Tyr157. We found that when Glu131 was modelled protonated, a rapid proton transfer (in 50 fs) occurred to Lys241 in contrast to no such event when it was modelled anionic. Our simulation data directly show how protonation of conserved Lys/Glu ion pair can drive proton transfer in forward direction by reducing activation energy barriers.

The third pathway in the R region is located between Lys211 of the Lys/Glu pair and conserved Tyr157. There are several structurally-resolved water molecules in the pathway (Fig. 1B), and these are expected to support Grotthuss-like proton transfer. However, unbiased QM/MM MD simulations with Glu131 deprotonated showed no forward proton transfer from Lys211 to Tyr157. Instead, the proton was found to stabilize on Glu131, and was accompanied by the loss of the water path between Lys211 and Tyr157. Only when Glu131 was modelled neutral, we observed a water bridge (two H₂O molecules) between Lys211 and Tyr157, yielding a rather low energy barrier (∼4 kcal mol⁻¹) of protonation dynamics, with the proton position on Tyr157 much more favorable (Fig. 3B). We thus suggest that similar to the inter-tyrosine protonation dynamics (Fig. 3C, right panel), protonation of Tyr157 can be modulated by the charged state of Glu131 of Lys/Glu ion-pair that resides at a distance of ∼9 Å from the tyrosine.

Data from classical MD simulations complement the picture that emerged from QM/MM MD simulations and show that protonated Lys211 sidechain drifts towards neutral Tyr157 upon protonation of Glu131/Glu66 pair residing at the interface (ESI, Fig. S2,† see more below). This behavior possibly facilitates the formation of a water path between Lys211 to Tyr157, which can benefit the proton transfer dynamics in this region. To track the energetics of proton transfer from Lys211 all the way to Tyr225, we performed an additional large-scale QM/MM free energy simulation with Glu131 protonated (see also Fig. S2†). It yielded an activation energy barrier of ∼6 kcal mol⁻¹ with around 2 kcal mol⁻¹ more favorable proton location on Tyr225 (Fig. 4; see also ESI, Figs. S8 and S9†). The metastable state at d ∼ −0.7 Å corresponds to the deprotonated Tyr157 prior its subsequent re-protonation and is in good agreement with the data from small-scale QM/MM free energy runs (Fig. 3B and C).

Fig. 4 Proton transfer in the R region. Proton transfer occurs from Lys211 (protonated) → H₂O → H₂O→ Tyr157 (neutral) → H₂O → Tyr225 (deprotonated) as depicted in the top panel. The bottom panel shows PMF (potential of mean force) profile of proton transfer derived from our QM/MM umbrella sampling simulations. Green shaded area around the bold green line depicts the bootstrapping error range (see also methods). Occupational histograms for the reaction coordinates are given in ESI, Fig. S17F.† The inset marks the studied QM region with dark green circle. The core membrane subunits ND4, ND2, and ND4L are shown in blue, pink, and green, respectively.

Inter-subunit proton transmission in the membrane arm of complex I

Next, we decided to estimate the protonation energetics at the subunit–subunit interface (ND2/ND4 and ND2/ND4L interfaces, see Fig. 1B). The high-resolution structure of mitochondrial complex I¹⁰ shows continuous water pathways between titratable residues residing at the subunit–subunit interface (Fig. 1B, S1C and D, ESI†). High-resolution structural data on complex I from other species also show similar interfacial hydration (ESI, Fig. S10†), suggesting that proton transfer can occur (see also ref. 24). Therefore, we next investigated the proton transfer dynamics on a chain of three water molecules spanning Lys383 and ^ND4Glu142 via hydrogen bonds. PMF profiles obtained showed a low energy barrier of ∼3.5 kcal mol⁻¹ with a much more favorable proton position on ^ND4Glu142 (Fig. 5A and S11†). Also, for the ND4L-ND2 interface, we observed a very low barrier for proton transfer from ^ND4LGlu66 to Glu131 of ∼1.5 kcal mol⁻¹ with around 4 kcal mol⁻¹ more advantageous proton location on Glu131 (Fig. 5B and S12, ESI†). Remarkably, this driven proton transfer caused ^ND4LGlu30 (modelled neutral), which was outside the chosen US reaction coordinate (ESI, Fig. S1D and Movie S3†), to protonate and neutralize ^ND4LGlu66, highlighting the cooperative nature of proton transfer dynamics, as observed also in other biological systems.²⁵

Fig. 5 Free energy profiles of the proton transfer in the “interface regions”. (A) ND2-ND4 interface: Lys383 (protonated) → H₂O → H₂O → H₂O → ^ND4Gu142 (deprotonated) pathway. (B) ND4L-ND2 interface: ^ND4LGlu66 (neutral) – H₂O – Glu131 (deprotonated) pathway. The bold lines are PMF profiles with light shaded area highlighting the bootstrapping error range (see also methods). Occupational histograms for the reaction coordinates are given in ESI, Fig. S17G and H.† The inset depicts the core membrane subunits of respiratory complex I: ND4 (blue), ND2 (pink), and ND4L (green). The green (A) and blue (B) circles denotes the QM region where the proton transfer is investigated.

To conclude, our multiscale simulation data reveal that inter- and intra-subunit proton transfer can occur with low activation energy barriers, and favorable thermodynamics. This has important implications on the proton pumping mechanism of complex I, as discussed below.

Discussion

By applying QM/MM free energy simulations on a high-resolution structure of complex I, we have identified how a proton can migrate through the central hydrophilic axis of an antiporter-like subunit of complex I. According to our results, proton migration can occur with a hole-like mechanism with reasonable kinetics and thermodynamics. Four conserved titratable residues play a central role in catalyzing proton transfer in the antiporter-like subunits of complex I. However, our QM/MM simulations reveal that conserved tyrosine residues present in these membrane-bound subunits may play a key role in catalyzing proton transfer and are thus suggested here to be primary candidates for site directed mutagenesis studies to probe complex I proton pumping (see also ref. 10). The results also show a profound impact of the charge state of glutamic acid of conserved Lys/Glu pair. The protonation of Glu131 can reduce the energetic barriers of proton transfers between different residues, highlighting that long-range electrostatic effects are at play (see also ref. 6 and 16). Data from our QM/MM simulations also reveal that an explicit proton transfer can occur from conserved Lys/Glu pair to the central lysine, a notion that has not been investigated earlier. Furthermore, our calculations for the first time show that proton transfer can occur across interfaces between membrane-bound subunits with reasonable energetics. Based on these combined data, we discuss below a step-by-step mechanism of proton transfer in the ND2 subunit of complex I.

Proton transfer mechanism in the antiporter-like subunit of complex I

According to the proton pumping mechanism proposed by Parey et al.,¹⁰ a pumped proton is released towards the central hydrophilic axis via the E channel upon electrostatic unlatching by quinone reduction and protonation at the second quinone binding site. This proton electrostatically drives the pumping of protons loaded on the antiporter-like subunits towards the P side via a single exit route present in terminal ND5 subunit, which has been studied biochemically, structurally, and computationally.^10,14,26 Based on recent structural^6,10 and computational data,¹⁰ the two other antiporter-like subunits ND2 and ND4 have been suggested to be devoid of P side proton exit routes (but see also, ²⁷). However, they do uptake protons from the N side and E channel, thus catalyzing horizontal proton transfers parallel to the membrane.¹⁰

To expand these conclusions, we propose the following mechanism of the lateral proton transfer in respiratory complex I (Fig. 6). Each antiporter-like subunit has a maximal proton/charge capacity, which is when achieved, requires an external trigger (such as an electrostatic push by protons released upon Q redox reaction²⁸) causing the release of the proton to the neighboring subunit and so on so forth to the P side via ND5. Conversely, if the maximal proton/charge capacity is not reached, an antiporter-like subunit would accept protons from the bulk medium or neighboring subunits. As a starting point, we propose a proton transfer occurs from the N side of the membrane to the neutral Lys241 because of its high pK_a relative to the bulk (steps 1 to 2, changing charge of antiporter-like subunit from +1 to +2). The pathways observed in earlier MD simulations^14–16,24 and recent high-resolution structural data support this scenario.¹⁰ Protonation of Lys241 causes excess charge to build up in ND2, which drives the release of a proton present on terminal Lys383 across ND2/4 interface to anionic ^ND4Glu142 with a small barrier of ∼3.5 kcal mol⁻¹ (Fig. 5A). The energetically favorable proton transfer across interface creates a proton vacancy on Lys383 (step 2).

Fig. 6 Proton transfer in the antiporter-like subunit of complex I. The proposed mechanism of lateral proton transfer in ND2 subunit of mitochondrial respiratory complex I, which is a direct continuation of the scheme proposed in ref. 10. Within one cycle shown here, the enzyme takes up one proton from the N-side (stage 1), and one from the E-channel (stage 4), and releases two protons to the ND4 subunit (stages 5b and 11). ND2 subunit is shown in brown, ND4 is in grey, and ND4L is in green. Red, blue, and green circles denote acidic, basic, and polar amino acids, respectively. A whole number at the bottom of each inset indicates the total charge of the depicted residues in ND2 subunit. Blue arrows show proton transfer and black dotted lines highlight separation between subunits.

In the next step, positively charged Lys241 donates its proton to terminal Lys383 (Fig. 2) causing a series of intra-subunit protonic rearrangements (step 3). Structural and classical MD simulation data suggested that the neutral state of central lysine is coupled with the closure of N-side connectivity by the movement of TMH11a/b and a conserved phenylalanine gate.¹⁰ At this stage, Lys241 waits for the second pumped proton to arrive via the horizontal axis from quinone binding tunnel. However, this requires proton transfers involving conserved tyrosine residues and does not occur until the second pumped proton arrives to ^ND4LGlu66 (step 4). A cascade of forward proton transfers (steps 5–8) takes place causing total charge to change from +1 to 0 resulting in charge deficit in ND2 subunit. Inter-subunit proton transfer occurs from ^ND4LGlu66 to Glu131 in step 8, which perturbs the inter-tyrosine protonation dynamics by long-range electrostatic effects (steps 8 to 9).

Finally, protonation of the anionic tyrosine (Tyr157) from Lys211 of Lys/Glu pair restores the system to ground state for a second cycle. In the cycle shown in Fig. 6, two pumped protons are taken up, one from the N-side (step 1) and one from the E-channel (step 4). In the proposed mechanism, protons diffuse through the entire horizontal axis of the antiporter-like subunit under the action of proton deficit (created by pumping) and excess (created by pumped proton released as part of Q redox chemistry). As part of this cycle, central lysine residues in antiporter-like subunits pick protons from the N side and may work cooperatively. Indeed, structural data revealed a state in which N side connectivity in ND2 was open, whereas the same path was closed in neighboring ND4 subunit.¹⁰ Overall, our results demonstrate that a lateral proton transfer can occur across the entire width of antiporter-like subunit of respiratory complex I via a Grotthuss-type “hole” scheme. This long range proton transfer in the membrane arm of complex I can be driven by protons injected from E channel as part of Q redox chemistry.^10,28

Methods

Hybrid quantum mechanical/molecular mechanical (QM/MM) simulations were performed on the high-resolution cryo-EM structure of mitochondrial complex I from Yarrowia lipolytica (PDB 7O71). The model system was limited to 6 core subunits: ND1, ND2, ND3, ND4, ND4L, ND6. The missing residues of the ND3 subunit were created by homology modelling using the MODELLER tool.²⁹ All structural water molecules that were within 5 Å of the above-mentioned protein subunits were taken into account. To avoid the large-scale motions of the long terminus of ND3 subunit (residues 114 to 128), this region was excluded from the system. The QM/MM system encompassed almost 28 700 atoms.

To track the proton transfer energetics and dynamics, we performed biased and unbiased QM/MM MD simulations on the ND2 subunit that is shown to be most hydrated in the high-resolution structure of complex I.¹⁰ To enhance sampling and keep model systems computationally tractable, we split the horizontal axis of ND2 subunit (comprising both amino acids and structurally-resolved waters) into separate QM regions as part of QM/MM setup (see Table S1†). For biased QM/MM MD simulations, umbrella sampling (US) technique combined with weighted histogram analysis method (WHAM) were applied.³⁰ QM/MM MD simulations were carried out using the combination of NAMD v2.14 (refs. ³¹ and 32) and ORCA 5.0.3 (ref. 33) software packages. After obtaining the structural coordinates from PDB, hydrogen atoms were added using PSFGEN VMD plugin.³⁴ Amino acids were modelled in their standard protonation states unless otherwise mentioned (see below). Topotools VMD plugin was used for setting up the QM regions.³⁵ The obtained model system was minimized at MM level for 200 steps with the steepest-descent algorithm, followed by QM/MM minimization for 500 steps.

In all simulations the temperature was maintained at 310 K using Langevin thermostat.³⁶ Trajectories were obtained using Verlet integration method³⁷ with a 1-fs timestep. No constraints were applied on covalent bonds to hydrogen atoms. Non bonded interactions were calculated using the Verlet³⁸ scheme with a 12 Å cutoff, 10 Å switching, and 14 Å pairlist distance. QM part was treated with density functional theory (DFT) with hybrid B3LYP functional,³⁹ and def2-SVP basis set.⁴⁰ For additional testing and benchmarking purposes (see more below), CAM-B3LYP functional⁴¹ with def2-TZVP basis set⁴⁰ was also used. To account for the long-range London dispersion interactions, the DFT-D3 dispersion correction⁴² was applied. The SCF energy tolerance was set to 10 × 10⁻⁸ au. Coulomb interaction between QM and MM parts of the system was implemented using additive electrostatic embedding scheme.³² To describe the MM framework, CHARMM36 (ref. 43) force field was used.

To perform US free energy simulations, we used the linear combination of the pathway-forming atomic bonds as a reaction coordinate d (ESI, Fig. S13†):

d = d₁ − d₂ + d₃ − d₄ (1)

In the first simulation window, the reaction coordinate was restrained to its initial position (d₀) by a harmonic potential:

where the force constant k was chosen to be equal to 100 kcal mol⁻¹ Å⁻². Every next window was created by shifting d₀ by Δd = 0.2 Å with the velocity 10 Å ps⁻¹. This procedure was implemented with the “colvars” NAMD module.^32,44 After the generation procedure, each window was equilibrated for 1.8 ps, followed by a 2-ps production run further utilized for free energy analysis. The PMF profiles were calculated with WHAM analysis tool.⁴⁵ The bootstrapping errors were obtained using 1000 Monte-Carlo trials with a 50-fs correlation time.

The total QM/MM MD simulation time in this work (including unbiased runs, and all US windows) is ∼1.5 ns. The VMD software package was used for the trajectory analysis and visualization.³⁴

To investigate the effect of density functional and basis set on free energy profiles, we performed free energy calculations on the pathways Tyr157-H₂O-Tyr225 (R region), and ^ND4LGlu66-H₂O-Glu131 (ND4L–ND2 interface) using CAM-B3LYP density functional and def2-TZVP basis set (ESI, Table S1 and Fig. S14, S15†). Only a minor difference in activation energy barrier (∼1.5 kcal mol⁻¹ higher) was observed, which suggests the choice of B3LYP density functional and def2-SVP basis set is reasonable for studying proton transfer dynamics and energetics. Moreover, this level of theory was previously validated in extensive QM/MM studies of proton transfer in similar biological systems.^46–48

Classical atomistic MD simulations were also performed on the complete structure of complex I from Yarrowia lipolytica (PDB:7O71). A detailed description of the model system setup can be found in ref. 10. Briefly put, the simulations were performed in two charged states of the protein; PN1 – in which titratable amino acids (Lys, Glu, Asp and Arg) were modeled in their charged states, whereas in state PN2 – protonation states of these amino acids were defined based on pK_a calculations, which were performed with the PROPKA software.^19,49 In the current work, all three simulation replicas of PN1 and PN2 states were extended to 1 μs using the same simulations conditions and GROMACS v2020.3 package.⁵⁰ In another state (PN3), a ∼150 ns MD simulation was performed with ND2 residues Lys383 and Glu131 modelled neutral. Analysis of the simulation trajectories was performed using VMD software.³⁴ RMSF (Fig. S16†) of key residues constituting the proton transfer path remains low (especially in PN2 state) and residues do not show large scale conformational changes from the observed structural conformation in microseconds time scales (1 μs × 3 replicas). This also supports the usage of high-resolution structural data for QM/MM MD simulations performed here.

Data availability

All the simulation trajectories and analysis data are available from the corresponding author upon reasonable request.

Author contributions

OZ performed QM/MM simulations, analyzed data, prepared figures, and wrote the initial draft of the manuscript. AD analyzed data, prepared figures, movies, and contributed to the discussions. JL performed classical MD simulations and contributed to manuscript writing. VS designed and supervised the research, prepared figures, wrote and refined the manuscript.

Conflicts of interest

The authors declare no competing interests.

Acknowledgements

Academy of Finland, Sigrid Jusélius Foundation, Jane and Aatos Erkko Foundation, University of Helsinki, Magnus Ehrnrooth Foundation. VS thanks Volker Zickermann for helpful discussions, and Outi Haapanen for assistance in preparation of simulation movies. OZ thanks Shreyas Kaptan for fruitful discussions on computational methodology. The authors acknowledge CSC IT center for science (Helsinki, Finland) for the high-performance supercomputing time.

References

1. L. A. Sazanov, A giant molecular proton pump: structure and mechanism of respiratory complex I, Nat. Rev. Mol. Cell Biol., 2015, 16(6), 375–388.
2. A.-N. A. Agip, et al., Mammalian respiratory complex I through the lens of cryo-EM, Annu. Rev. Biophys., 2019, 48, 165–184.
3. I. Chung, et al., Cryo-EM structures define ubiquinone-10 binding to mitochondrial complex I and conformational transitions accompanying Q-site occupancy, Nat. Commun., 2022, 13(1), 1–13.
4. D. N. Grba and J. Hirst, Mitochondrial complex I structure reveals ordered water molecules for catalysis and proton translocation, Nat. Struct. Mol. Biol., 2020, 27(10), 892–900.
5. J. Gu, et al., The coupling mechanism of mammalian mitochondrial complex I, Nat. Struct. Mol. Biol., 2022, 29(2), 172–182.
6. D. Kampjut and L. A. Sazanov, The coupling mechanism of mammalian respiratory complex I, Science, 2020, 370(6516), eabc4209.
7. J. Gutiérrez-Fernández, et al., Key role of quinone in the mechanism of respiratory complex I, Nat. Commun., 2020, 11(1), 1–17.
8. K. Parey, et al., High-resolution cryo-EM structures of respiratory complex I: Mechanism, assembly, and disease, Sci. Adv., 2019, 5(12), eaax9484.
9. P. Kolata and R. G. Efremov, Structure of Escherichia coli respiratory complex I reconstituted into lipid nanodiscs reveals an uncoupled conformation, Elife, 2021, 10, e68710.
10. K. Parey, et al., High-resolution structure and dynamics of mitochondrial complex I—Insights into the proton pumping mechanism, Sci. Adv., 2021, 7(46), eabj3221.
11. A. Tursun, et al., Probing the proton channels in subunit N of Complex I from Escherichia coli through intra-subunit cross-linking, Biochim. Biophys. Acta, Bioenerg., 2016, 1857(12), 1840–1848.
12. M. Narayanan, et al., Roles of subunit NuoL in the proton pumping coupling mechanism of NADH: ubiquinone oxidoreductase (complex I) from Escherichia coli, J. Biochem., 2016, 160(4), 205–215.
13. F. Hoeser, et al., Respiratory complex I with charge symmetry in the membrane arm pumps protons, Proc. Natl. Acad. Sci., 2022, 119(27), e2123090119.
14. A. Djurabekova, et al., Proton motive function of the terminal antiporter-like subunit in respiratory complex I, Biochim. Biophys. Acta, Bioenerg., 2020, 1861(7), 148185.
15. M. E. Mühlbauer, et al., Water-gated proton transfer dynamics in respiratory complex I, J. Am. Chem. Soc., 2020, 142(32), 13718–13728.
16. M. Röpke, et al., Functional water wires catalyze long-range proton pumping in the mammalian respiratory complex I, J. Am. Chem. Soc., 2020, 142(52), 21758–21766.
17. V. Kravchuk, et al., A universal coupling mechanism of respiratory complex I, Nature, 2022, 609(7928), 808–814.
18. J. Michel, et al., Mutagenesis of the L, M, and N subunits of complex I from Escherichia coli indicates a common role in function, PLoS One, 2011, 6(2), e17420.
19. C. R. Søndergaard, et al., Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa values, J. Chem. Theory Comput., 2011, 7(7), 2284–2295.
20. C. J. T. de Grotthuss, Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique, Ann. Chim., 1806, 58, 54–73.
21. H. Ishikita, Origin of the pKa shift of the catalytic lysine in acetoacetate decarboxylase, FEBS Lett., 2010, 584(15), 3464–3468.
22. H. Ishikita and K. Saito, Proton transfer reactions and hydrogen-bond networks in protein environments, J. R. Soc., Interface, 2014, 11(91), 20130518.
23. G. Hong, et al., On understanding proton transfer to the biocatalytic [Fe–Fe] H sub-cluster in [Fe–Fe] H2ases: QM/MM MD simulations, Biochim. Biophys. Acta, Bioenerg., 2011, 1807(5), 510–517.
24. O. Haapanen and V. Sharma, Role of water and protein dynamics in proton pumping by respiratory complex I, Sci. Rep., 2017, 7(1), 1–12.
25. Yu Sugo, K. Saito and H. Ishikita, Conformational changes and H-bond rearrangements during quinone release in photosystem II, Biochemistry, 2022, 61(17), 1836–1843.
26. E. Nakamaru-Ogiso, et al., The membrane subunit NuoL (ND5) is involved in the indirect proton pumping mechanism of Escherichia coli complex I, J. Biol. Chem., 2010, 285(50), 39070–39078.
27. N. Klusch, et al., Cryo-EM structure of the respiratory I+ III2 supercomplex from Arabidopsis thaliana at 2 Å resolution, Nat. Plants, 2022, 1–15.
28. O. Haapanen and V. Sharma, A modeling and simulation perspective on the mechanism and function of respiratory complex I, Biochim. Biophys. Acta, Bioenerg., 2018, 1859(7), 510–523.
29. A. Fiser and A. Šali, Modeller: generation and refinement of homology-based protein structure models, Methods Enzymol., 2003, 374, 461–491.
30. S. Kumar, et al., The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method, J. Comput. Chem., 1992, 13(8), 1011–1021.
31. J. C. Phillips, et al., Scalable molecular dynamics on CPU and GPU architectures with NAMD, J. Chem. Phys., 2020, 153(4), 044130.
32. M. C. R. Melo, et al., NAMD goes quantum: an integrative suite for hybrid simulations, Nat. Methods, 2018, 15(5), 351–354.
33. F. Neese, et al., The ORCA quantum chemistry program package, J. Chem. Phys., 2020, 152(22), 224108.
34. W. Humphrey, A. Dalke and K. Schulten, VMD: visual molecular dynamics, J. Mol. Graphics, 1996, 14(1), 33–38.
35. A. Kohlmeyer, The TopoTools VMD plugin, 2017.
36. G. S. Grest and K. Kremer, Molecular dynamics simulation for polymers in the presence of a heat bath, Phys. Rev. A, 1986, 33(5), 3628.
37. M. P. Allen, and D. J. Tildesley, Computer simulation of liquids, Clarendon, Oxford, 1987.
38. L. Verlet, Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev., 1967, 159(1), 98.
39. A. D. Becke, Becke's three parameter hybrid method using the LYP correlation functional, J. Chem. Phys., 1993, 98(492), 5648–5652.
40. F. Weigend and R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, Phys. Chem. Chem. Phys., 2005, 7(18), 3297–3305.
41. T. Yanai, D. P. Tew and N. C. Handy, A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP), Chem. Phys. Lett., 2004, 393(1–3), 51–57.
42. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu, J. Chem. Phys., 2010, 132(15), 154104.
43. R. B. Best, et al., Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone ϕ, ψ and side-chain χ1 and χ2 dihedral angles, J. Chem. Theory Comput., 2012, 8(9), 3257–3273.
44. G. Fiorin, M. L. Klein and J. Hénin, Using collective variables to drive molecular dynamics simulations, Mol. Phys., 2013, 111(22–23), 3345–3362.
45. A. Grossfield, WHAM: the weighted histogram analysis method, version 2.0.10, http://membrane.urmc.rochester.edu/wordpress/?page_id=126.
46. S. L. Dürr, et al., The Role of Conserved Residues in the DEDDh Motif: the Proton-Transfer Mechanism of HIV-1 RNase H, ACS Catal., 2021, 11(13), 7915–7927.
47. W.-H. Deng, Y. Lu and R.-Z. Liao, Revealing the Mechanism of Isethionate Sulfite-Lyase by QM/MM Calculations, J. Chem. Inf. Model., 2021, 61(12), 5871–5882.
48. P. Saura, et al., Electric fields control water-gated proton transfer in cytochrome c oxidase, Proc. Natl. Acad. Sci., 2022, 119(38), e2207761119.
49. M. H. M. Olsson, et al., PROPKA3: consistent treatment of internal and surface residues in empirical p K a predictions, J. Chem. Theory Comput., 2011, 7(2), 525–537.
50. M. J. Abraham, et al., GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX, 2015, 1, 19–25.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sc01427d

---