Mechanistic Insights on Heme-to-Heme Transmembrane Electron Transfer Within NADPH Oxydases From Atomistic Simulations

Introduction

NADPH oxidases (NOX) encompass a family of membrane enzymes dedicated to the production of cellular reactive oxygen species (ROS) in a regulated manner (Geiszt and Leto, 2004; Bedard and Krause, 2007). NOX enzymes have conserved structures common to all family members. In Humans, seven isoenzymes are expressed: NOX1 to NOX5, Duox1 and Duox2. This family is subdivided into two groups; one according to the ability of the NOX to form a heterodimer with p22^phox (NOX1-4), the other according to the presence of calcium-binding EF-hand type motifs in its sequence (NOX5 and Duox1/2). These enzymes are expressed in many tissues, including kidney, fibroblasts, osteoclasts, and thyroid, where they produce ROS, albeit generally at diverse levels, in response to stimuli such as growth factors, cytokines or calcium. They are involved in a great number of physiological functions, such as host defense, post-translational proteins processing, inter- and intra-cellular signaling, regulation of gene expression and cell differentiation. Aberrant levels of NADPH-oxidase derived ROS, either too low or too high, can disturb the balance of cellular homeostasis, ultimately resulting in pathological states. NOX are widely distributed in different kingdoms of life including animals (in particular mammals), plants and lower organisms and their expression is specific to each kingdom of life (Bedard et al., 2007; Sumimoto, 2008; Hajjar et al., 2017).

Whereas the members of the NADPH oxidase family differ in their tissue distribution, mode of activation and physiological functions, they all share the capacity to generate ROS. Among them, superoxide is produced by NOX1-3 and NOX5, and H₂O₂ is produced by NOX4 and Duox enzymes. NOX are flavohemoproteins and electron transporters consisting in a N-terminus region of six transmembrane helices that bind two non-identical heme groups and a cytosolic C-terminus region (dehydrogenase domain) where an NAD(P)H and a Flavin Adenosine Diphosphate (FAD)-binding site are localized. This transmembrane protein permits to shuttle electrons brought by NADPH molecules from the cytosol toward the internal part of the phagosome (Figure 1). Directional electron transfers are catalyzed by a series of redox molecules encapsulated within the protein matrices that serve as stepping stones for the electrons. On the cytosolic side, FAD is the primary electron acceptor and gets a hydride anion, i.e. two electrons and one proton, from NADPH. It is hypothesized that electrons are then shuttled one-by-one to the first heme (Heme1), then to the second heme (Heme2), before being finally transferred to a dioxygen molecule to form the superoxide anion. Redox potentials of hemes in NOX2, which is the most studied member of the NADPH oxidase family and is considered as the prototype of NOX family, have been measured (Cross et al., 1995). To the best of our knowledge, there have been no reports of kinetic measurements of these intrinsic ET steps which are presumably faster than other rate-limiting steps of the overall NOX machinery (e.g., sub-units association or NADPH binding).

FIGURE 1

Figure 1. Electron transfer across the membrane catalyzes the production of superoxide anions. This image has been generated using the NOX5 model reported in Magnani et al. (2017). The yellow and green protein domains are the dehydrogenase and transmembrane domains respectively.

This sketchy description of transmembrane ET is supported by various biochemical data, notably by recent 3D structures. In 2017, Magnani et al. reported the first X-ray diffraction structures of the dehydrogenase and transmembrane domains of NOX5 from Cylindrospermum stagnale. This cyanobacteria NOX5 exhibits 40% sequence homology with human NOX5 (Magnani et al., 2017). Very recently, cryo-EM structures of mouse Duox1–DuoxA1 have been reported (Sun, 2020). NOX5 and Duox catalyze superoxide and hydrogen peroxide production, respectively. Their 3D-structures are well conserved and present two well-aligned hemes in the transmembrane domain (Figure 1). A phenylalanine residue (379 in NOX5 and 1097 in DUOX1-DUOXA1) is located midway between the two hemes suggesting this residue could play an important role to sustain inter-heme electron tunneling. These ET path schemes in NOX enzymes are valuable hypotheses but should be taken with caution as reasoning on single structures to predict tunneling pathways, ignoring structural fluctuations, has often proven to be misleading (Prytkova et al., 2007; de la Lande et al., 2010, 2013; Beratan et al., 2015). Putative binding pockets for O₂ molecules have been identified which could serve to maintain molecular oxygen close to Heme2, thus facilitating its reduction (Magnani et al., 2017; Sun, 2020).

Globally, we start to accumulate structural insights on NOX enzymes. However, little is known about the thermodynamics and kinetics of the succession of electron transfer steps within NOX. In this article, we report the first microscopic simulations of NOX5 within a realistic lipid membrane. Our aim is twofold. First, we aim at building a realistic model of NOX5 at the atomic level. We report molecular dynamics simulations of several hundreds of nanoseconds to assess its stability in a fully hydrated lipid bilayer environment. Second, we report a computational study of the inter-heme electron transfer step including evaluations of the ET free energy, of reorganization free energies and of tunneling pathways. We further reveal that O₂ moves back-and-forth between two binding pockets situated near Heme2. When Heme2 is reduced, O₂ dominantly populates the pocket which is the closest to the heme. The electron tunneling pathways afforded by this pocket are one order of magnitude stronger than those emerging when O₂ is present in the more remote pocket, hence providing a way to catalyze the final electron transfer to produce superoxide. Data and scripts to produce key figures and table of this article are available on the Zenodo data base (doi: 10.5281/zenodo.4424142).

Building of an Atomistic Model for NOX5

Our aim is to investigate inter-heme electron transfer. To this end, our approach is to develop an atomistic model of NOX5 inserted within a lipid bilayer based on the recently reported structures of helical transmembrane (TM) and C-terminal cytosolic dehydrogenase (DH) domains from Cylindrospermum stagnale NOX5. The initial coordinates of the TM and DH domains were taken from the Protein Data Bank structures with codes 5O0T (2.0 Å resolution) and 5O0X (2.2 Å resolution), respectively. We start by detailing the protocol we have followed to set-up the system and to carry out molecular dynamics simulations. In the following section we will focus on electron transfer modeling.

Molecular Dynamics Simulation of NOX Proteins

Preparation of the Simulation System

The MD simulations were carried out for an extended system comprising the NOX protein catalytic core embedded in a model phospholipid membrane and solvated with an ionic aqueous solution as shown in Figure 2. Magnani et al. (2017) identified a putative cavity to bind a dioxygen molecule in the vicinity of Heme2 (see the red bead on Figure 2, right). In our simulations this cavity is filled either by a dioxygen molecule or by a water molecule (see below).

FIGURE 2

Figure 2. Left, NOX5 catalytic core investigated in this study, composed of the transmembrane domain (in blue) and dehydrogenase domain (in orange) embedded in a lipid bilayer (in light blue and red). Water molecules are shown in gray. Right, insert shows the cofactor FAD (purple), Heme1 (green), Heme2 (yellow) and O2 (red).

Construction of a Full-Length NOX5 Model

In the experiments that achieved crystallization of the DH domain, a short amino acid sequence PWLELAAA was added after the C-terminal Phe693 residue to enhance thermal stability and FAD retention (Magnani et al., 2017). We removed this artificial additional sequence in our model. On the other hand, 27 amino acid residues were missing in the N-terminal calcium-binding EF-hand domain (GAQQKSDMKSSTLFVAMDLMHQETKVD), for which no homology modeling template was available. We used the PSIPRED server (Buchan and Jones, 2019) to predict the secondary structure of this segment. The central part of this sequence is predicted to form an alpha helix. MODELLER (Šali and Blundell, 1993) was used to add this missing part, however the resulting models did not form an alpha-helix. This segment is only present in the NOX5 branch of the family, and is replaced by e.g., an 11-residue stretch in a NOX4. So, we constructed a NOX4-like variant, by deleting SSTLFVAMDLMHQET and only keeping GAQQKSDMKVD (Resid 605–615). This shorter sequence was modeled ab initio with MODELLER. 20 models were built and assessed based on their DOPE (Discrete Optimized Protein Energy) score (Shen and Sali, 2006). Differences between the 20 models were small: the RMSD difference between the model with best score and worst is 0.19 Å. We considered that 20 models are enough to build this missing part and we chose the model with the best score. We do not expect this assumption to influence the study of the ET process between two hemes, since the distance between the reconstructed loop and Heme1 is 35 Å.

In ref. Magnani et al. (2017), the X-ray structures of the two DH and TM domains were obtained separately, and docked together with some user-defined restraints using the Haddock server (Dominguez et al., 2003; van Zundert et al., 2016). We have followed the same docking strategies. Specifically, Asn288 (B-loop of TM), Lys361 (D-loop of TM), Thr520 (B-loop binding region of DH), and Phe677 (flavin-interacting C terminus of DH) were defined as “active” residues (preferentially placed at the interface by Haddock). Moreover, we constrained the edge-to-edge distance between FAD and Heme1 to be within 2.5 Å. We also constrained Phe677 to be close to the isoalloxazine ring of FAD with a face-to-face π-stacking interaction (Figure 2). The distance between Lys412 (C terminus of TM) and Glu413 (N terminus of DH) was constrained to 2 Å. The best docked complex (Haddock score of −2.5) was then fused by creating a covalent bond between Lys412 and Glu413. The geometry of the model was further validated by the Qmean server (Benkert et al., 2011) for model quality estimation. Our model obtained a QMEAN Z-score of −1.97. We can consider this score as satisfactory since the Duox1 cryo-EM structure has a slightly worse Z-score of −2.65.

The interface between the transmembrane and dehydrogenase domains was the main unknown when building our full-length model. Recent cryo-EM structures of the homologous Duox1 protein contain both domains, and offer a basis for comparison (Sun, 2020). We find that when aligning the transmembrane domain, our model of the dehydrogenase domain is rotated by 60° with respect to the equivalent domain in the 7D3E structure of the human DUOX1-DUOXA1 complex (Supplementary Figure 1). During the simulations, the relative orientation of the two domains fluctuates by up to 11° from the initial model.

Construction of the Extended System

The docking structure was used as initial coordinates of the catalytic core. All amino acid residues were assigned a standard protonation state at pH 7. All histidine residues were singly protonated at the Nδ sites based on the inspection of inter-residue interactions. The importance of the lipid metabolism was evidenced in NOX activities, although the regulatory role of the lipid membrane properties remains unclear. Membrane properties (charges, thickness, composition, lipid raft formation) have been shown to affect the NOX2 enzyme functioning (Shao et al., 2003; Souabni et al., 2014, 2017). Moreover, the phospholipid composition changes during NOX2 activation (Magalhaes and Glogauer, 2010; Bréchard et al., 2012; Joly et al., 2020). We made the choice of building a membrane with 1-Palmitoyl-2-oleoyl-sn-glycero-3-(phospho-rac-(1-glycerol)) (POPG) and 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine (POPE) in 4 to 1 proportions with CHARMMGUI (Jo et al., 2008; Brooks et al., 2009; Lee et al., 2016). Note that this membrane composition is one among many possibilities that could have been considered. The OPM (orientations of proteins in membranes) orientation of the TM domain has been used as starting point for membrane insertion. The protein membrane system was solvated with 38,341 water molecules in a hexagonal prism unit cell (a,b,c = 135Å,117Å,112 Å, alpha, beta, gamma = 90, 90, 120°). The system was neutralized electrically by adding 451 Na+ and 98 Cl− ions to obtain a 0.15 M NaCl solution.

Force Field Parameters

Standard protein residues and lipids were modeled using the CHARMM36 force field (Lee et al., 2016). We used the TIP3P water model (MacKerell et al., 1998). Force field parameters for the heme cofactor in ferrous and ferric states and for FADH• have been developed by our groups and can be found in Supplementary Materials. In this study, we consider two distinct redox states associated to the inter-hemes electron transfer (see next section): Heme1-Fe(II)/Heme2-Fe(III) and Heme1-Fe(III)/Heme2-Fe(II). In some simulations, a single dioxygen molecule was placed in the cavity shaped by the propionate groups of Heme2 and by residues Arg256, His313 and His317 as shown in Figure 2. This cavity has been previously proposed to host dioxygen prior to its reduction. Following the methodology proposed in Ref. Javanainen et al. (2017, 2020), we placed a charged dummy atom holding a 0.452 electric charge in between the oxygen nuclei, each holding a −0.226 electric charge in order to reproduce the quadrupole moment of the dioxygen molecule. A restraint energy term was added to the total potential energy with the Colvar module (Fiorin et al., 2013) of NAMD (Phillips et al., 2005, 2020) to maintain the dioxygen molecule within the pocket during the MD simulations (details are given later in the text).

Simulation Protocol

Initial equilibration and production phase were carried out using the NAMD program (Phillips et al., 2020). We carried out MD simulations for two different systems. One system includes the presence of O₂ in the aforementioned cavity, while the other is without dioxygen. We conducted in each case MD simulations for the system in the initial and final redox states, represented by different partial charge distributions on heme and iron moieties. The protein membrane system was equilibrated with a standard procedure as suggested from CHARMM-GUI (Jo et al., 2008; Lee et al., 2016). The equilibration phase has been carried out in the isothermal–isobaric ensemble (NPT) under periodic boundary conditions with six cycles of 25, 25, 25, 200, 200, 200 ps by gradually releasing harmonic restraints around the initial positions of protein backbone and membrane lipid atoms at 310 K and 1 bar. For protein, the harmonic restraints force constants of each cycle were set to 10.0, 5.0, 2.5, 1.5, 1.0, and 0.5 kcal/mol/Ų, while membrane lipid restraints were applied with force constants of 5.0, 5.0, 2.0, 1.0, 0.2, and 0 kcal kcal/mol/Ų. The integration time step was set to 1 fs for the first three cycles and then increased to 2 fs for the last three cycles. All bonds involving H are fixed in length. To maintain planarity of the membrane, the center of mass of the lipid head groups in the upper and lower layers had their positions restrained in the z direction to +19 and −19 Å relative to the center of mass of the bilayer by application of a constraining potential with a force constant of 5.0, 5.0, 2.0, 1.0, 0.2, and 0 kcal/mol/Ų. The particle-mesh Ewald (PME) method was used for the calculation of electrostatic interactions. Non-bonding interactions were treated using a cutoff of 12.0 Å. Pressure was controlled by the Nose-Hoover Langevin piston method (Phillips et al., 2005) while the temperature was controlled by Langevin dynamics.

Production runs in the canonical ensemble (NVT) were carried out for 300 ns for each redox state. Geometries were saved every 2 ps to sample the vertical energy gap. We also conducted two 100 ns MD simulations for each state with the same equilibrated structure but different initial velocities to assess the variability of the computed ET parameters with respect to the initial conditions of the simulation. Altogether, we generated a set of 12 MD simulations, 6 for each state (with and without oxygen), with 3 replicas of 300, 100, and 100 ns for each state, respectively, and two charge distributions each representing the initial and final redox states.

Analysis of Protein Stability

The protein backbone root-mean-square deviation (RMSD) time series calculated during MD simulations are shown in Supplementary Figure 2. The docking structure is taken as reference for the RMSD calculation. Panels (a) to (c) correspond to the MD with oxygen present, while panels (d) to (f) are without O₂. Black (resp. red) curves correspond to dynamics performed in the initial (resp. final) redox states. For all simulations, the RMSD stays under around 4 Å, assessing a slight distortion with respect to the initial structural model. Only in the 300 ns-long MD simulation without O₂ in the initial redox state a small drift toward higher values of RMSD is seen at the end of the simulation.

Analysis of NOX5 in the Membrane

We have assessed the stability of the insertion of the transmembrane domain of the protein in the lipid bilayer by computing the membrane thickness and the relative position of the protein in the bilayer throughout the simulations. The membrane thickness is estimated by the distance between the center-of-mass of upper and lower lipids. It is stable during all the MD simulations with a fluctuation of less than 2 Å (Supplementary Figure 3). We use the distance between the center of mass of the membrane and the center of mass of the transmembrane part of the protein to measure the relative position of the protein inside the bilayer. Its evolution in all the MD simulations is shown in Supplementary Figure 4. It is found to be stable with fluctuations less than 6 Å and no apparent drift on the hundreds of ns timescale.

Heme-to-heme Electron Transfer

Having assessed the stability of the in-silico model of NOX5 inserted in a lipid membrane, we focus in this section on the heme-to-heme electron transfer. We follow the framework of the Marcus Theory of electron transfer to evaluate the free energy of the reaction (ΔA°) and the reorganization energy (λ) (Marcus and Sutin, 1985). Warshel et al. showed that both quantities can be obtained from microscopic simulations under the Linear Response Approximation as (King and Warshel, 1990):

$\begin{array}{ll}\Delta A^0=\frac{1}{2}({\langle \Delta E\rangle }_i{+\langle \Delta E\rangle }_f)& (1)\end{array}$

$\begin{array}{ll}{\lambda }^{St}=\frac{1}{2}({\langle \Delta E\rangle }_i{-\langle \Delta E\rangle }_f)& (2)\end{array}$

where ΔE is the diabatic energy gap, that is the difference in potential energy of the system in the initial and final redox state (ΔE = E_f−E_i), or alternatively said when the transferred electron sits on Heme1 or Heme2, respectively. Note that ΔE is calculated for a given set of nuclear coordinates and is therefore referred to as a vertical energy gap. 〈…〉_x refers to an average performed over a canonical ensemble of configurations of the system in redox state x (i or f). ΔE can be decomposed into a contribution coming from the heme cofactors (inner-sphere contribution), from the environment (outer-sphere contribution) and from a coupling term between them. The latter refers to the mutual polarization (E^mp) that differs for the two redox states: ΔE = ΔE^is+ΔE^os+ΔE^mp. Neglecting ΔE^mp, a separation of the total free energy of the reaction as ΔA = ΔA^is+ΔA^os is obtained. This approximation has been tested in other heme proteins and turned out to be reasonable (Blumberger, 2008). We adopt it here, too. Furthermore, as ET is taking place between two chemically identical hemes, ΔA^is = 0, therefore ΔA≈ΔA^os.

The reorganization energy calculated with Equation 2 is referred to as the Stokes reorganization energy, hence the upper script St. The fluctuation of the energy gap can also be used to define the reorganization energy:

$\begin{array}{ll}{\lambda }_x^{var}=\frac{var(\Delta E_x)}{2k_{B}T}& (3)\end{array}$

In Equation 3, ΔE_x is the energy gap computed over a trajectory performed in the state x (i or f). From MD simulations performed on each electronic state involved in the ET, one obtains two different values: ${\lambda }_i^{var}$ and ${\lambda }_f^{var}$. If the LRA and the ergodic principle apply, the three definitions of the reorganization energy coincide: ${\lambda }^{St}={\lambda }_i^{var}={\lambda }_f^{var}$. A simple measure of non-ergodicity can be obtained as:

$\begin{array}{ll}{\chi }_G=\frac{{\lambda }_i^{var}+{\lambda }_f^{var}}{2{\lambda }^{St}}& (4)\end{array}$

Values of χ_G close to 1 correspond to ergodic systems.

Neglecting ΔE^mp as above, reorganization energy can be decomposed into two contributions: λ = λ^is+λ^os. λ^os has been evaluated by Eqs. 2 and 3 using force field-based interaction energies computed between the two hemes and the environment. λ^is has been evaluated by Density Functional Theory calculation (see Supplementary Material for details). The geometries of the heme cofactors including the apical and axial histidine ligands (modeled as methyl-imidazoles) have been optimized and λ^is has been calculated as:

$\begin{array}{ll}{\lambda }^{is}=2.\left(\frac{E_{gR}^O-E_{gR}^R+E_{gO}^R-E_{gO}^O}{2}\right)& (5)\end{array}$

$E_{gX}^X$ stands for the energy of the heme in the redox state X(O∨R) in the equilibrium geometry gX. We obtained a value of 0.10 eV for λ^is. For each system (with or without O₂) we have computed ET thermodynamic parameters using the full 300 ns-long MD simulations, discarding the first 20 ns. In order to assess the variability of ET parameters with respect to the sampling, we have split the 300 ns-long simulations into three segments of 100 ns. Adding the two small 100 ns-long simulations, we have access to 5 sets of data that we consider to be independent, hereafter referred to as Set1 to Set 5. From these datasets, we have computed ET thermodynamic parameters, discarding each time the first 20ns, considered as an equilibration time. A block analysis over the 5 obtained sets of values has finally been performed. The results obtained for Heme1 to Heme2 ET are summarized in Figures 3, 4 (all the individual numerical values are given in Supplementary Table 1).

FIGURE 3

Figure 3. Evolution of the vertical energy gap ΔE (E_f−E_i) along MD simulations with O₂ (A) and without O₂ (B). Each line section represents 80ns of MD simulation. The first three sections correspond to parts of the 300ns-long simulation (Set 1–3) while the two last sections come from the two independent 100ns-long replica simulations (Set 4 and 5). Energy gaps computed based on simulations in the initial (resp. final) redox state are shown in red (resp.green). Data shown correspond to running averages of ΔE for clarity. Horizontal lines are guides to the eye.

FIGURE 4

Figure 4. ET reaction free energies and outer-sphere reorganization energies of the ET between hemes with O₂ (A) and without O₂ (B). The gray bars represent values obtained for the five 100ns-long data sets. The green bars are the results of a block average using the data of these five data sets. The red bars correspond to values obtained using the full 300ns-long trajectory. Error bars correspond to twice the uncertainty.

Figure 3 presents the evolution of the energy gaps ΔE in the five data sets Set1 to Set5 with O₂ (panel a) or without O₂ (panel b), computed in the initial (in red) and final (in green) redox states. We begin by describing the results in presence of dioxygen. The data clearly shows that fluctuations of the energy gap are greater in the initial redox state (in red). What is more, a clear bimodality appears for the energy gaps in this redox state, with ΔE oscillating around two values, which are indicated with horizontal dotted lines in Figure 3. This oscillation seems to take place at the scale of tens of ns. No clear-cut localized change in structure of the protein and/or its environment could be found to account for this bimodality. That said, we decomposed ΔE into different contributions, namely TM domain, DH domain, environment (water, membrane and counter ions) and cofactors (FAD, hemes, O₂). The result of this decomposition is shown in Supplementary Figures 5, 6 for simulations in presence and absence of O₂ respectively. For the initial redox state we recover a bimodality when considering the contributions from the TM domain and from the environment. Longer simulations could help to decipher the origin of the bimodality of ΔE. It could come from slow concerted motion involving the TM and the lipids composing the membrane and its solvation layer. On the other hand, the other components (FAD, Hemes, O₂, DH domain) don't exhibit any bimodality in the energy gap distributions. In addition, these results are independent on the presence or not of dioxygen in the cavity.

ET reaction free energies obtained with Equation 1 are negative whatever data set is used to compute them, with values ranging roughly from −0.5 to −0.1 eV (see Figure 4 and Supplementary Table 1), indicating that the ET from Heme1 to Heme2 is favorable. The values obtained from the full 300 ns-long MD simulations or from the five 100 ns-long datasets are consistent (green and red bars in Figure 4). Because of the larger fluctuations in ΔE, values of ΔA⁰ obtained without O₂ are more dispersed (see e.g., Figure 4). However, the average values obtained for the two systems are similar within the uncertainties of the calculations (−0.28 eV and −0.24 eV respectively with and without O₂). Outer-sphere reorganization energies are presented in Figure 4, calculated with Equation 2 (λ^St) or Equation 3 (${\lambda }_i^{var}$ and ${\lambda }_f^{var}$). λ^Stvalues lie between 1 and 1.4 eV, with averages of 1.04 and 1.24 eV for λ^St respectively with and without O₂. Note that these reorganization energies are computed with a non-polarizable forcefield. Accounting for electrostatic induction is expected to decrease the values of the reorganization energies by 30% approximately (Blumberger, 2008). Once again one observes greater fluctuations in the values when O₂ is not bound to Heme2. This is especially the case when one considers ${\lambda }_x^{var}$ values computed over the 300 ns-long MD simulations in the initial electronic state. This observation should be taken with caution as the duration of the simulation is likely not to be long enough with respect to the slow degree of freedom that is responsible for the bimodality in ΔE. Whichever way they are computed, reorganization energies are greater in absence of dioxygen (see Supplementary Table 1). This dependency may arise because the presence of O₂ “rigidifies” the environment around Heme2, lowering the structural differences between both electronic states, noting that the uncertainties of the mean values are large, respectively 0.07 and 0.24 eV. Further investigations with longer simulations or more replicas should help to ascertain these conclusions more strongly. Values of ${\lambda }_x^{var}$ are found to be greater than λ^St for every dataset, leading to χ_G values between 1.1 and 2.0, most of the time lying between 1.3 and 1.5. For simulations without O₂, the long 300 ns MD simulations gives rise to a much larger ${\lambda }_i^{var}$(ca. 3.5 eV, red bar in Figure 4, right). This high value, associated to a large statistical uncertainly (1.17 eV, Supplementary Table 1) results from the bimodality in the energy gap distribution which is apparent in Figure 3. Together with the low values obtained for λ^os, this finding reveals that the structural reorganization around the two hemes is somewhat moderate, since much larger χ_G values have been reported (Matyushov, 2015).

We have decomposed ΔA⁰ and λ^St into contributions arising from the molecular components of the system, namely the TM and DH domains, FAD, the Hemes, the counter-ions, the membrane and the solvent (Supplementary Table 2). It turns out that the negative ΔA⁰ results from a balance between positive and negative contributions. With the current membrane composition (POPG/POPE in 4:1 proportion), the TM domain and membrane largely disfavor inter-heme ET (1.52 ± 0.18 and 0.60 ± 0.14 eV, respectively, considering the average over the five sets of trajectories), while water and counter-ions provide a net driving force for ET (−0.27 ± 0.12 and −1.42 ± 0.68, respectively). We can thus extrapolate that a change of membrane composition, especially if it incorporates more of less charged lipids would directly impact inter-heme ET. Interestingly, the DH domain has little influence on ΔA⁰ and marginally reorganizes upon electron transfer. As in our model the DH domain might have a non-physiological orientation with respect to the TM domain in view of the recent cryo-EM structure of DUOX1 (see above, and Supplementary Figure 1), this result is reassuring regarding our estimates of inter-heme ET thermodynamics, i.e., the precise orientation of the DH domain marginally determines these parameters, likely because the DH is far away from the inter-heme region.

Electron Tunneling Pathways

We have carried out tunneling pathway analyses along MD trajectories with the pathway model of Beratan and co-workers. This model is built on the assumption that the electron tunnels from the electron donor to the acceptor along a pathway that is defined as a succession of covalent or hydrogen bonds and of through space contacts (Beratan et al., 1987). The decay of the electronic coupling associated with a given pathway is expressed as the product of multiplication of a constant contact coupling ($H_{if}^{contact}$) by a factor (ε_tot) that reflects the attenuation of the coupling caused by the presence of intervening medium between the donor and the acceptor. ε_tot depends on the number of covalent or hydrogen bonds (denoted N_c and N_hb, respectively) and of through-space jumps (N_ts) composing the electron transfer pathway. This factor is calculated according to a set of mathematical expressions empirically calibrated by the authors of the model (Onuchic and Beratan, 1990; Beratan et al., 1991).

$\begin{array}{ll}{H}_{if}=H_{if}^{contact}{\epsilon }_{tot}& (6)\end{array}$

$\begin{array}{ll}{\epsilon }_{tot}={\displaystyle \prod _{N_c}}{\epsilon }_c\times {\displaystyle \prod _{N_{hb}}}{\epsilon }_{hb}\times {\displaystyle \prod _{N_{ts}}}{\epsilon }_{ts}& (7)\end{array}$

$\begin{array}{ll}{\epsilon }_c=0.6& (8)\end{array}$

$\begin{array}{ll}{\epsilon }_{hb}=0.36\times exp\left[-{\beta }_S(R_H-2.8)\right]& (9)\end{array}$

$\begin{array}{ll}{\epsilon }_{ts}=0.6\times exp\left[-{\beta }_S(R_S-1.4)\right]& (10)\end{array}$

R_S is the atom–atom distance for through-space jumps, R_H is the hydrogen bond length (taking the inter-heavy atom distance) and β_S is a characteristic decay factor. The latter was originally parametrized in the 1990's to a value of 1.7Å⁻¹, but later investigations suggested that a value of 1.1Å⁻¹ enables better agreement with quantum chemistry calculations (Prytkova et al., 2005). This is the value we have used here. Covalent bonds within the porphyrin aromatic rings were made fully conductive (i.e., setting ε_c = 1 between two such atoms) in our analysis, so that the PM results do not depend on the particular choice of the donor/acceptor atom within the two hemes. The best pathway was identified for each analyzed structure using the Dijkstra algorithm (Dijkstra, 1959) with an in-house program (de la Lande et al., 2010; El Hammi et al., 2012). Our program (i) seeks for the pathway producing the most efficient electronic coupling with the Dijkstra algorithm, (ii) decomposes the total pathway over covalent, through-hydrogen-bond and through-space interactions (iii) indicates which amino-acid residues or other molecules are involved in tunneling pathways, (iv) provides adequate files for visual representation of the pathways with the VMD program. Structures were extracted from MD simulation trajectories every 20 ps. To make the pathway search tractable, the structure was preliminarily pruned to include only the atoms susceptible to take part in the tunneling process, that is residues not lying too far from the direct inter-heme axis. The pruned sub-system comprises the two hemes and their axial ligands (His299, His313, His372, His385) as well as resides Leu268, Ile269, Lys300, Leu301, Val302, Gly303, Gln304, Val305, Met306, Phe307, Ala308, Leu309, Ala310, Ile311, Val312, Leu345, Leu346, Val347, Phe348, Ile349, Ile350, Met351, Trp352, Trp378, Phe379, Trp392 (Figure 5). As seen on the picture, the hemes are separated by a layer of five aromatic amino acids, among which Trp378, as highlighted in Magnani et al. (2017) but also Phe307, Phe348, Trp352, Phe379, and Trp392.

FIGURE 5

Figure 5. Left, restricted molecular system considered in ET pathway analyses. The hemes are shown in red. Right, for clarity only the closest residues from the Donor-Acceptor axis are represented (tryptophan residues in gray and phenylalanine residues in violet).

Results are collected in Table 1 and on Figure 6. Two types of pathways are identified over the duration of the MD simulations. The first one is running through Trp378, while the other is running through Phe348 as illustrated by the green and red tubes in Figure 6. They connect the hemes by the edge of the porphyrin plane. Our simulations thus confirm the expectation that Trp378 plays a key role in mediating inter-heme tunneling, but also reveal the occurrence of pathways running through Phe348. As seen on the distributions depicted in Figure 6, the tunneling pathways running through Phe348 are in general slightly more efficient than those running through Trp378. The respective distributions among the two kinds of pathways seem to depend on the redox state but it is difficult to find a rational in terms of structure of the intervening medium, and both kinds of pathways give rise to rather similar decay factors. Overall, the average decay factor equals to 2.11 10⁻⁵ with a standard deviation of 9.98 10⁻⁶ when analyzing MD trajectories of the Fe²⁺/Fe³⁺ redox state, and to 3.03 10⁻⁵ and 1.52 10⁻⁵ for the Fe³⁺/Fe²⁺ redox state.

TABLE 1

Table 1. Semi-empirical pathway analyses of tunneling matrix elements mediating inter-heme electron transfer. $\langle H_{if}^{PM}\rangle$ are given in eV².

FIGURE 6

Figure 6. Left, illustration of the two types of pathways encountered for heme-to-heme ET in NOX5 running either through Trp378 (in green) or through Phe348 (in red). Right, distributions of tunneling decay factors for MD simulations in the Fe²⁺/Fe³⁺ and Fe³⁺/Fe²⁺ redox states (in dashed and plane lines respectively), deconvoluted by pathway types.

As $\langle {{\epsilon }_{tot}}^2\rangle ={\langle {\epsilon }_{tot}\rangle }^2+{\sigma }^2$ with ${\sigma }^2=\langle {({\epsilon }_{tot}-\langle {\epsilon }_{tot}\rangle )}^2\rangle$, the ratio $R_{coh}^{cl}$ defined as $R_{coh}^{cl}={\langle {\epsilon }_{tot}\rangle }^2/\langle {{\epsilon }_{tot}}^2\rangle$ permits to assess if tunneling is governed by fluctuations of the electronic coupling or by its average (Balabin and Onuchic, 2000; Balabin et al., 2012). If $R_{coh}^{cl}$ approaches unity it means that ${\langle {\epsilon }_{tot}\rangle }^2\gg {\sigma }^2$, tunneling is governed by the average value. On the other hand, if $R_{coh}^{cl}$ approaches 0.5, it means that ${\langle {\epsilon }_{tot}\rangle }^2\ll {\sigma }^2$, i.e. tunneling is likely to be dominated by the fluctuations. We obtain values larger than 0.8 suggesting tunneling is actually dominated by a single kind of pathway in both cases. This observation is consistent with the fact that almost all pathways run through Trp378 or Phe348. Note that $R_{coh}^{cl}$ is a classical coherence-like parameter in the sense that it is calculated from ε_tot provided by the semi-empirical pathway model. $R_{coh}^{cl}$ does not explicitly capture fine quantum effects like interferences among competing tunneling pathways running through the intervening medium (Prytkova et al., 2007). The pathway does not give access to H_if, but to the decay ε_tot caused by the presence of an intervening medium (see the first equation). A crude estimation of the absolute coupling may still be obtained by setting $H_{if}^{contact}$ to 4.3 10¹³ Hz (0.177829 eV) (as in Prytkova et al., 2005), yielding an average value of 3.75 10⁻⁶ eV for $\langle H_{if}^{PM}\rangle$ from MD trajectories of the Fe²⁺/Fe³⁺ state.

Probing the Putative O₂-accessible Cavity and Tunneling Pathways

The authors of Magnani et al. (2017) identified a putative cavity to bind the dioxygen molecule in the vicinity of Heme2 (right insert in Figure 2). In the simulation without O₂, we saw that a water molecule occupied this cavity, that is shaped by residues Arg256, His313, and His317 and Heme2 (Figure 7). We replaced this water molecule with a dioxygen molecule to perform MD simulations with O₂. First attempts to run simulations without constraint on O₂ position lead to its escape from the protein matrix. This may be due to the fact that there is only one dioxygen molecule in the system and/or to the fact that the protein has not had enough time to accommodate to the O₂ molecule. Since we are interested in the impact of O₂ binding on the dynamics of the protein and on the ET parameters, we decided to constrain O₂ to remain in the cavity. To this end we used Colvar options in NAMD. We defined a distance between the center of mass (COM) of the cavity and O₂ [hereafter referred to as d(O₂-COM)]. We added a boundary semi-harmonic potential restraint applying on this distance when it is more than 5 Å and with a force constant of 20 kcal/mol/Ų.

FIGURE 7

Figure 7. (A) and (B) Images showing the two favored O₂ positions with respect to the binding cavity. The protein backbone is represented in gray ribbon and the cavity is highlighted in all-atom representation: Arg256, His313, and His317 and Heme2, O₂ (in red). The cavity in gray bubble. (C) Histograms of the distance between O₂ and the center of mass (COM) of the cavity in the initial (black) and final (red) redox states. Results from the 300ns-long MD runs.

Figure 7C shows the histograms of this distance in the initial and final electronic states for the 300 ns-long MD simulations. A bimodal structure appears with two different positions of O₂ around the cavity. The oscillations of O₂ between these two positions occur at the timescale of few ns. The first position is O₂ inserted inside the cavity with d(O₂-COM) of around 2.5 Å (Figure 7A) and a distance between O₂ and the Fe atom of Heme2 (d(O₂-Fe)) around 5.5 Å. The other one corresponds to O₂ lying on the edge of the cavity with d(O₂-COM) of around 4.3 Å and d(O₂-Fe) around 8.5 Å (Figure 7B). The relative proportion of occupancy of the two sites is reversed in the two electronic states. After the electron transfer (in red), O₂ is more often close to Heme2 than in the initial state (black).

To explore possible implications of the existence of these two connected cavities for electron tunneling pathways connecting Heme2 to O₂, we have carried out pathway analyses as in the previous section, but enlarging the pruned system to encompass amino acid residues from the binding pocket. We have analyzed all the MD simulations (in the initial and final electronic states of the inter-heme ET). Two kinds of pathways are identified. The “short pathway” type implies a direct through-space jump from the porphyrin ring to O₂. It is illustrated in Figure 8 (middle). The “long pathway” involves an intermediate passage through His317 before jumping to O₂. The short pathways are associated to an average decay factor of 0.065 while the long pathways are associated to an average decay factor of 0.004, which is directly correlated to the respective lengths of the pathways. Interestingly, the short pathways are predominant when Heme2 is reduced, while the long pathways are dominant when Heme2 is oxidized. This could be a piece of the mechanism ensuring superoxide production catalysis. When Heme2 is reduced, thus ready to transfer an electron to a dioxygen molecule, the latter binds closer to the porphyrin ring, hence enabling very fast tunneling. The molecular origin of the dynamical control of the position shifting is not completely clear. From our simulations, it appears that the positively charged arginine Arg256 and histidine His317 play a role in this process. The conformations accessible to Arg256 in both redox states are similar. There is no frank residue repositioning upon Heme2 reduction. On the other hand, while Arg256 is almost exclusively hydrogen bonded to the propionate moieties in the ferrous state, it also often makes alternative hydrogen bonds to His317 in the ferric state (Supplementary Figure 8). In such situations, the cavity for O₂ binding close to Heme2 is inaccessible. Further simulations will be needed to fully explain this behavior.

FIGURE 8

Figure 8. Electron tunneling pathway analysis between Heme2 and the dioxygen molecule. Top, normalized distributions of decay factors collected along the MD simulations in initial and final redox states of the inter-heme ET [(Heme2(III)) and Heme2(II) respectively]. Middle and bottom, representative snapshots of the two dominant pathways encountered during MD simulations.

Discussion and Conclusion

This study reports the first molecular simulations of a NOX system inserted into a lipid bilayer membrane. Starting from separate experimental structures of the DH and TM domains of NOX5, we have built a fused model of the full protein. This model is stable in MD simulations at the timescale of hundreds of nanoseconds. The interface between the two domains involves the same residues as in the very recent cryo-EM structures of the homologous protein DUOX1, however the DH domain has a different orientation in our model. We have focused on electron transfer between Heme1 and Heme2 within the TM domain, and on the impact of O₂ on ET thermodynamics parameters. Dioxygen binding in the vicinity of Heme 2 appears to reduce the fluctuations of the energy difference between the two redox states involved in the ET, leading to more stable values of ET parameters. In both cases, O₂ binding doesn't significantly alter the thermodynamics of inter-herme electron transfer. Negative values of the redox free energies on the order of few tenths of eV have been obtained, corresponding to spontaneous ET. Outer-sphere reorganization energies of about 1.1–1.3 eV (including contributions of the inner and outer-spheres) have been observed. Decomposition of these thermodynamics parameters over the molecular components of the system has provided additional insight. First, the DH has little impact on the inter-heme ET, suggesting that the relative orientation of the DH with respect to the cryo-EM structure is not a major concern for the present study of inter-heme ET. Second, we have highlighted the subtle balance existing between the protein matrix, the membrane, the water and the counter-ions. Our results should stimulate other studies aiming at investigating for instance the effect of membrane composition on inter-heme ET. Considering that outer-sphere reorganization energies have been calculated with a non-polarizable forcefield and are thus overestimated by ca. 30%, the values obtained fall within the range of values (0.7–1.0 eV) obtained recently for inter-heme ET within multi-hemic protein (Breuer et al., 2012). We have finally evaluated the tunneling coupling matrix element to fall within the μeV based on the tunneling pathway model (Beratan et al., 1991; Prytkova et al., 2005). We have been able to confirm the important role of Trp378 and Phe348 to sustain electron tunneling between the two hemes.

Dioxygen reduction by electron transfer is another essential part of NOX5 function that remains unclear though of prominent importance to understand the NOX machinery. It was suggested from Electron Paramagnetic Resonance and Raman spectroscopies (Yamaguchi et al., 1989; Hurst et al., 1991; Ueno et al., 1991; Fujii et al., 1999) that since Heme2 is hexacoordinated Heme2-to-O₂ ET is likely to proceed via an outer-sphere mechanism. This hypothesis was corroborated by stopped-flow and rapid-scanning spectroscopy oxidation-reduction kinetics on porcine neutrophils (Isogai et al., 1995). The identification of a putative binding cavity in the 3D structure of NOX5 further supports this hypothesis (Magnani et al., 2017). In NOX5, the cavity is formed by Arg256, His313, and His317. We have investigated the localization of dioxygen near Heme2 by means of molecular dynamics simulations. Our simulations revealed that the cavity deduced from the X-ray structure is actually flexible, offering two possible binding sites, one being closer to the heme than the other, while the other is more accessible to solvent. The occupancy of these two sites depends on the redox states of the two hemes. When Heme2 is in a reduced state, the close position of O₂ is largely favored. Importantly, this occupancy shift is coupled to the emergence of ten-times more efficient electron tunneling pathways than when O₂ is positioned in the remote cavity. This 10-fold increase of tunneling matrix element would induce a 100-fold increase of ET rate form Heme2 to O₂, therefore providing a powerful mean to catalyze superoxide production. We speculate that Arg256 is playing a central role in this process, since this positively charged residue necessarily rearranges when Heme2 is reduced. On the other hand, our MD simulations don't reveal any dramatic cavity shape modifications, but more subtle rearrangements. The stability of the hydrogen bond network among Arg256, His317, and the Heme2's propionate groups that shape the cavity appears to be sensitive to the charge present on Heme2. Further investigations will be needed to fully decipher this mechanism. Our study thus provides a possible mechanistic explanation to the experimental results reported in Magnani et al. (2017) indicating that the oxidation rate in R256S and H317R mutants were at least fivefold lower than observed for wild type. We now need to investigate the thermodynamics of this ET, including the possibility of spin-inversion on Heme2 cofactor as suggested by some authors (Doussière et al., 1996; Maturana et al., 2001), as well as the mechanism of primary binding of dioxygen into the remote cavity when Heme2 is oxidized.

This work is the first step toward a global investigation of the mechanism of superoxide production by the NOX machinery using molecular simulation. Many aspects deserve further analysis (computation of kinetic parameters, study of the electron transfers between FAD and Heme1, influence of lipid composition and/or mutation of key residues sustaining ETs in NOX…). The refinement of computational methods (such as the use of polarizable forcefields) will allow us to ascertain the conclusions drawn in this study, but the results obtained here open the path to a better understanding of NOX.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: doi: 10.5281/zenodo.4424142; https://zenodo.org/record/4424142#.X_ciUyXjLeM/b.

Author Contributions

MB, JH, and AL designed research. XW performed research. XW, JH, FC and AL carried out complementary analyses. All authors discussed and analyzed research data. All authors wrote the manuscript.

Funding

Computational work was performed using HPC resources from GENCI (Grant Number A0070701714) and from LBT-HPC thanks to support from Geoffrey Letessier. This work was supported by the Initiative d'Excellence program from the French State (Grant DYNAMO, ANR-11-LABX-0011 and Grant CACSICE, ANR-11-EQPX-0008). MB thanks Sesame Ile-de-France for co-funding the display wall used for data analysis.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Prof. Mattevi who provided us with early pdb file containing atomic coordinates of the NOX5 3D model they built and that turned very useful to initiate our work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.650651/full#supplementary-material

References

Balabin, I. A., Hu, X., and Beratan, D. N. (2012). Exploring biological electron transfer pathway dynamics with the pathways plugin for VMD. J. Comput. Chem. 33, 906–910. doi: 10.1002/jcc.22927

PubMed Abstract | CrossRef Full Text | Google Scholar

Balabin, I. A., and Onuchic, J. N. (2000). Dynamically controlled protein tunneling paths in photosynthetic reaction centers. Science 290, 114–117. doi: 10.1126/science.290.5489.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Bedard, K., and Krause, K.-H. (2007). The NOX family of ROS-generating NADPH oxidases: physiology and pathophysiology. Physiol. Rev. 87, 245–313. doi: 10.1152/physrev.00044.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bedard, K., Lardy, B., and Krause, K.-H. (2007). NOX family NADPH oxidases: not just in mammals. Biochimca 89, 1107–1112. doi: 10.1016/j.biochi.2007.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Benkert, P., Biasini, M., and Schwede, T. (2011). Toward the estimation of the absolute quality of individual protein structure models. Bioinformtics, 27, 343–350. doi: 10.1093/bioinformatics/btq662

PubMed Abstract | CrossRef Full Text | Google Scholar

Beratan, D. N., Betts, J. N., and Onuchic, J. N. (1991). Protein electron transfer rates set by the bridging secondary and tertiary structure. Science 252, 1285–1288. doi: 10.1126/science.1656523

PubMed Abstract | CrossRef Full Text | Google Scholar

Beratan, D. N., Liu, C., Migliore, A., Polizzi, N. F., Skourtis, S. S., Zhang, P., et al. (2015). Charge transfer in dynamical biosystems, or the treachery of (static) images. Acc. Chem. Res. 48, 474–481. doi: 10.1021/ar500271d

PubMed Abstract | CrossRef Full Text | Google Scholar

Beratan, D. N., Onuchic, J. N., and Hopfield, J. J. (1987). Electron tunneling through covalent and noncovalent pathways in proteins. J. Chem. Phys. 86, 4488–4498. doi: 10.1063/1.452723

CrossRef Full Text | Google Scholar

Blumberger, J. (2008). Free energies for biological electron transfer from QM/MM calculation: method, application and critical assessment. Phys. Chem. Chem. Phys.10, 5651–5667. doi: 10.1039/b807444e

PubMed Abstract | CrossRef Full Text | Google Scholar

Bréchard, S., Plançon, S., and Tschirhart, E. J. (2012). New insights into the regulation of neutrophil NADPH oxidase activity in the phagosome: a focus on the role of lipid and Ca2+ signaling. Antioxid. Redox Signal. 18, 661–676. doi: 10.1089/ars.2012.4773

CrossRef Full Text | Google Scholar

Breuer, M., Zarzycki, P., Shi, L., Clarke, T. A., Edwards, M. J., Butt, J. N., et al. (2012). Molecular structure and free energy landscape for electron transport in the decahaem cytochrome MtrF. Biochem. Soc. Trans. 40, 1198–1203. doi: 10.1042/BST20120139

PubMed Abstract | CrossRef Full Text | Google Scholar

Brooks, C. L., Brooks, A. D. III., Mackerell, L. Jr, Nilsson, R. J., Petrella, B., Roux, Y., et al. (2009). The biomolecular simulation program. J. Comput. Chem. 30, 1545–1614. doi: 10.1002/jcc.21287

CrossRef Full Text | Google Scholar

Buchan, D. W. A., and Jones, D. T. (2019). The PSIPRED protein analysis workbench: 20 years on. Nucl. Acid Res. 47, W402–W407. doi: 10.1093/nar/gkz297

PubMed Abstract | CrossRef Full Text | Google Scholar

Cross, A. R., Rae, J., and Curnutte, J. T. (1995). Cytochrome b of the neutrophil superoxide-generating system contains two nonidentical hemes: POTENTIOMETRIC STUDIES OF A MUTANT FORM OF Gp91. J. Biol. Chem. 270, 17075–17077. doi: 10.1074/jbc.270.29.17075

CrossRef Full Text | Google Scholar

de la Lande, A., Babcock, N. S., Rezáč, J., Sanders, B. C., and Salahub, D. R. (2010). Surface residues dynamically organize water bridges to enhance electron transfer between proteins. Proc. Natl. Acad. Soc. 107, 11799–11804. doi: 10.1073/pnas.0914457107

PubMed Abstract | CrossRef Full Text | Google Scholar

de la Lande, A., Babcock, N. S., Rezáč, J., Sanders, B. C., and Salahub, D. R. (2013). Surface residues dynamically organize water bridges to enhance electron transfer between proteins. Proc. Natl. Acad. Soc. 110, 1136–1137. doi: 10.1073/pnas.1220833110

PubMed Abstract | CrossRef Full Text | Google Scholar

Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numer. Math. No. 1, 269:BF01386390. doi: 10.1007/BF01386390

CrossRef Full Text | Google Scholar

Dominguez, C., Boelens, R., and Bonvin, A. M. J. J. (2003). HADDOCK: a protein–protein docking approach based on biochemical or biophysical information. J. Am. Chem. Soc. 125, 1731–1737. doi: 10.1021/ja026939x

PubMed Abstract | CrossRef Full Text | Google Scholar

Doussière, J., Gaillard, J., and Vignais, P. V. (1996). Electron transfer across the O²- generating flavocytochrome b of neutrophils. Evidence for a transition from a low-spin state to a high-spin state of the heme iron component. Biochemistry 35, 13400–13410. doi: 10.1021/bi960916b

PubMed Abstract | CrossRef Full Text | Google Scholar

El Hammi, E., Houee-Levin, C., Rezac, J., Levy, B., Demachy, I., Baciou, L., et al. (2012). New insights into the mechanism of electron transfer within flavohemoglobins: tunnelling pathways, packing density, thermodynamic and kinetic analyses. Phys. Chem. Chem. Phys. 14, 13872–13880. doi: 10.1039/c2cp41261f

PubMed Abstract | CrossRef Full Text | Google Scholar

Fiorin, G., Klein, M. L., and Hénin, J. (2013). Using collective variables to drive molecular dynamics simulations. null 111, 3345–3362. doi: 10.1080/00268976.2013.813594

CrossRef Full Text | Google Scholar

Fujii, H., Finnegan, M. G., and Johnson, M. K. (1999). The active form of the ferric heme in neutrophil cytochrome B558 is low-spin in the reconstituted cell-free system in the presence of amphophil1. J. Biochem. 126, 708–714. doi: 10.1093/oxfordjournals.jbchem.a022507

CrossRef Full Text | Google Scholar

Geiszt, M., and Leto, T. L. (2004). The NOX family of NAD(P)H oxidases: host defense and beyond. J. Biol. Chem. 279, 51715–51718. doi: 10.1074/jbc.R400024200

PubMed Abstract | CrossRef Full Text | Google Scholar

Hajjar, C., Cherrier, M. V., Dias Mirandela, G., Petit-Hartlein, I., Stasia, M. J., Fontecilla-Camps, J. C., et al. (2017). The NOX family of proteins is also present in bacteria. mBio, 8:e01487–e01417. doi: 10.1128/mBio.01487-17

CrossRef Full Text | Google Scholar

Hurst, J. K., Loehr, T. M., Curnutte, J. T., and Rosen, H. (1991). Resonance raman and electron paramagnetic resonance structural investigations of neutrophil cytochrome B558. J. Biol. Chem. 266, 1627–1634. doi: 10.1016/S0021-9258(18)52340-1

CrossRef Full Text | Google Scholar

Isogai, Y., Iizuka, T., and Shiro, Y. (1995). The mechanism of electron donation to molecular oxygen by phagocytic cytochrome B558. J. Biol. Chem. 270, 7853–7857. doi: 10.1074/jbc.270.14.7853

CrossRef Full Text | Google Scholar

Javanainen, M., Vattulainen, I., and Monticelli, L. (2017). On atomistic models for molecular oxygen. J. Phys. Chem. B, 121, 518–528. doi: 10.1021/acs.jpcb.6b11183

CrossRef Full Text | Google Scholar

Javanainen, M., Vattulainen, I., and Monticelli, L. (2020). Correction to “on atomistic models for molecular oxygen”. J. Phys. Chem. B 124, 6943–6946. doi: 10.1021/acs.jpcb.0c06376

PubMed Abstract | CrossRef Full Text | Google Scholar

Jo, S., Kim, T., Iyer, V. G., and Im, W. (2008). CHARMM-GUI: a web-based graphical user interface for CHARMM. J. Comput. Chem. 29, 1859–1865. doi: 10.1002/jcc.20945

PubMed Abstract | CrossRef Full Text | Google Scholar

Joly, J., Hudik, E., Lecart, S., Roos, D., Verkuijlen, P., Wrona, D., et al. (2020). Membrane dynamics and organization of the phagocyte NADPH oxidase in PLB-985 cells. Front. Cell Dev. Biol. 8:1295. doi: 10.3389/fcell.2020.608600

PubMed Abstract | CrossRef Full Text | Google Scholar

King, G., and Warshel, A. (1990). Investigation of the free energy functions for electron transfer reactions. J. Chem. Phys. 93, 8682–8692. doi: 10.1063/1.459255

CrossRef Full Text | Google Scholar

Lee, J., Cheng, X., Swails, J. M., Yeom, M. S., Eastman, P. K., Lemkul, J. A., et al. (2016). CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J. Chem. Theory Comput. 12, 405–413. doi: 10.1021/acs.jctc.5b,00935

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKerell, A. D., Bashford, D., Bellott, M., Dunbrack, R. L., Evanseck, J. D., Field, M. J., et al. (1998). All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 102, 3586–3616. doi: 10.1021/jp973084f

PubMed Abstract | CrossRef Full Text | Google Scholar

Magalhaes, M. A. O., and Glogauer, M. (2010). Pivotal advance: phospholipids determine net membrane surface charge resulting in differential localization of active Rac1 and Rac2. J. Leuk. Biol. 87, 545–555. doi: 10.1189/jlb.0609390

PubMed Abstract | CrossRef Full Text | Google Scholar

Magnani, F., Nenci, S., Millana Fananas, E., Ceccon, M., Romero, E., Fraaije, M. W., et al. (2017). Crystal structures and atomic model of NADPH oxidase. Proc. Natl. Acad. Soc. 114, 6764–6769. doi: 10.1073/pnas.1702293114

PubMed Abstract | CrossRef Full Text | Google Scholar

Marcus, R. A., and Sutin, N. (1985). Electron transfers in chemistry and biology. Biochim. Biophys. Acta 811, 265–322. doi: 10.1016/0304-4173(85)90014-X

CrossRef Full Text | Google Scholar

Maturana, A., Arnaudeau, S., Ryser, S., Banfi, B., Hossle, J. P., Schlegel, W., et al. (2001). Heme histidine ligands within Gp91 Phox modulate proton conduction by the phagocyte NADPH oxidase. J. Biol. Chem. 276, 30277–30284. doi: 10.1074/jbc.M010438200

PubMed Abstract | CrossRef Full Text | Google Scholar

Matyushov, D. (2015). Protein electron transfer: is biology (thermo) dynamic? J. Phys. Condens. Matter, 27:473001. doi: 10.1088/0953-8984/27/47/473001

PubMed Abstract | CrossRef Full Text | Google Scholar

Onuchic, J. N., and Beratan, D. N. (1990). A predictive theoretical model for electron tunneling pathways in proteins. J. Chem. Phys. 92, 722–733. doi: 10.1063/1.458426

CrossRef Full Text | Google Scholar

Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., et al. (2005). Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 1781–1802. doi: 10.1002/jcc.20289

CrossRef Full Text | Google Scholar

Phillips, J. C., Hardy, D. J., Maia, J. D. C., Stone, J. E., Ribeiro, J. V., Bernardi, R. C., et al. (2020). Scalable molecular dynamics on CPU and GPU architectures with NAMD. J. Chem. Phys. 153:044130. doi: 10.1063/5.0014475

PubMed Abstract | CrossRef Full Text | Google Scholar

Prytkova, T. R., Kurnikov, I. V., and Beratan, D. N. (2005). Ab initio based calculations of electron-transfer rates in metalloproteins. J. Phys. Chem. B 109, 1618–1625. doi: 10.1021/jp0457491

PubMed Abstract | CrossRef Full Text | Google Scholar

Prytkova, T. R., Kurnikov, I. V., and Beratan, D. N. (2007). Coupling coherence distinguishes structure sensitivity in protein electron transfer. Science 315, 622–625. doi: 10.1126/science.1134862

PubMed Abstract | CrossRef Full Text | Google Scholar

Šali, A., and Blundell, T. L. (1993). Comparative protein modelling by satisfaction of spatial restraints. J. Mol. Biol. 234, 779–815. doi: 10.1006/jmbi.1993.1626

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, D., Segal, A. W., and Dekker, L. V. (2003). Lipid rafts determine efficiency of NADPH oxidase activation in neutrophils. FEBS Lett. 550, 101–106. doi: 10.1016/S0014-5793(03)00845-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, M., and Sali, A. (2006). Statistical potential for assessment and prediction of protein structures. Prot. Sci. 15, 2507–2524. doi: 10.1110/ps.062416606

PubMed Abstract | CrossRef Full Text | Google Scholar

Souabni, H., Machillot, P., and Baciou, L. (2014). Contribution of lipid environment to NADPH oxidase activity: influence of sterol. Biochemistry 107, 33–42. doi: 10.1016/j.biochi.2014.10.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Souabni, H., Wien, F., Bizouarn, T., Houée-Levin, C., Réfrégiers, M., and Baciou, L. (2017). The physicochemical properties of membranes correlate with the NADPH oxidase activity. Biochim. Biophys. Acta. 1861, 3520–3530. doi: 10.1016/j.bbagen.2016.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Sumimoto, H. (2008). Structure regulation and evolution of NOX-family NADPH oxidases that produce reactive oxygen species. FEBS J. 275, 3249–3277. doi: 10.1111/j.1742-4658.2008.06488.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, J. (2020). Structures of mouse DUOX1–DUOXA1 provide mechanistic insights into enzyme activation and regulation. Nat. Struct. Mol. Biol. 27, 1086–1093. doi: 10.1038/s41594-020-0501-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ueno, I., Fujii, S., Ohya-Nishiguchi, H., Iizuka, T., and Kanegasaki, S. (1991). Characterization of neutrophil B-type cytochrome in situ by electron paramagnetic resonance spectroscopy. FEBS Lett. 281, 130–132. doi: 10.1016/0014-5793(91)80375-D

PubMed Abstract | CrossRef Full Text | Google Scholar

van Zundert, G. C. P., Rodrigues, J. P. G. L. M., Trellet, M., Schmitz, C., Kastritis, P. L., Karaca, E., et al. (2016). The HADDOCK2.2 web server: user-friendly integrative modeling of biomolecular complexes. J. Mol. Biol. 428, 720–725. doi: 10.1016/j.jmb.2015.09.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamaguchi, T., Hayakawa, T., Kaneda, M., Kakinuma, K., and Yoshikawa, A. (1989). Purification and some properties of the small subunit of cytochrome B558 from human neutrophils. J. Biol. Chem. 264, 112–118. doi: 10.1016/S0021-9258(17)31230-9

PubMed Abstract | CrossRef Full Text | Google Scholar