Cytochrome aa3 Oxygen Reductase Utilizes the Tunnel Observed in

Mar 16, 2018 - Cytochrome aa3 is the terminal respiratory enzyme of all eukaryotes and many bacteria and archaea, reducing O2 to water and harnessing ...
0 downloads 9 Views 19MB Size
Subscriber access provided by - Access paid by the | UCSB Libraries

Cytochrome aa3 oxygen reductase utilizes the tunnel observed in the crystal structures to deliver O for catalysis 2

Paween Mahinthichaichan, Robert B. Gennis, and Emad Tajkhorshid Biochemistry, Just Accepted Manuscript • DOI: 10.1021/acs.biochem.7b01194 • Publication Date (Web): 16 Mar 2018 Downloaded from http://pubs.acs.org on March 16, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Cytochrome aa3 oxygen reductase utilizes the tunnel observed in the crystal structures to deliver O2 for catalysis Paween Mahinthichaichan,‡ Robert B. Gennis,∗,§ and Emad Tajkhorshid∗,‡ ‡Department of Biochemistry, NIH Center for Macromolecular Modeling and Bioinformatics, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A. §Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A. E-mail: [email protected]; [email protected] Phone: +1 217 3339075; +1 217 2446914

Running header O2 delivery pathway in cytochrome aa 3

1 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Funding Information The authors acknowledge funding supports from the National Institutes of Health (NIH P41GM104601, U01-GM111251 and U54-GM087519 to E.T., and R01-HL16101 to R.B.G.) and the Office of Naval Research (ONR N00014-16-1-2535 to E.T.). P.M. gratefully acknowledges NIH support as a trainee of the Molecular Biophysics Training Program (5T32-GM008276) during his graduate study.

2 ACS Paragon Plus Environment

Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Abstract Cytochrome aa 3 is the terminal respiratory enzyme of all eukaryotes and many bacteria and archaea, reducing O2 to water and harnessing the free energy from the reaction to generate the transmembrane electrochemical potential. The diffusion of O2 to the heme-copper catalytic site, which is buried deep inside the enzyme, is the initiation step of the reaction chemistry. Our previous molecular dynamics (MD) study with cytochrome ba 3 , a homologous enzyme of cytochrome aa 3 in Thermus thermophilus, demonstrated that O2 diffuses from the lipid bilayer to its reduction site through a 25-˚ A long tunnel inferred by Xe-binding sites de1 tected by X-ray crystallography. Although a similar tunnel is observed in cytochrome aa 3 , this putative pathway appears partially occluded between the entrances and the reduction site. Also, the experimentally determined second-order rate constant for O2 delivery in cytochrome aa 3 (∼108 M−1 s−1 ) is 10 times slower than that in cytochrome ba 3 (∼109 M−1 s−1 ). A question to be addressed is whether cytochrome aa 3 utilizes this X-ray inferred tunnel as the primary pathway for O2 delivery. Using complimentary computational methods including multiple independent flooding MD simulations and implicit ligand sampling calculations, we probe the O2 delivery pathways in cytochrome aa 3 of Rhodobacter sphaeroides. All of the O2 molecules that arrived in the reduction site during the simulations were found to diffuse through the X-ray observed tunnel, despite its apparent constriction, supporting its role as the main O2 delivery pathway in cytochrome aa 3 . The rate constant for O2 delivery in cytochrome aa 3 , approximated using the simulation results, is 10 times slower than in cytochrome ba 3 , in agreement with the experimentally determined rate constants.

3 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction A-family heme-copper oxygen reductases (HCOs) including the aa 3 -type cytochrome c oxidases (cytochrome aa 3 ) are the terminal respiratory oxygen reductases of all eukaryotes and many bacteria and archaea. 2–5 These enzymes oxidize cytochrome c and reduce O2 to water, and use the free energy from the O2 reduction reaction to power the biosynthesis of ATP. 6–12 The atomic structures of cytochrome aa 3 from Rhodobacter sphaeroides (R.s), 13,14 Paracoccus denitrificans 15 and bovine heart mitochondria 16–19 have been determined by X-ray crystallography. The bacterial enzymes comprise 3-4 subunits (Fig. 1A), whereas the mitochondrial enzymes contain at least ten additional accessory subunits. 20,21 The smaller bacterial enzymes are widely used as models for the mitochondrial enzymes. A homologue of cytochrome aa 3 is cytochrome bo 3 of E. coli , 5,22,23 which uses ubiquinol instead of cytochrome c as the electron donor. The reduction of O2 takes place in Subunit I (SI), a 12-TM helical subunit, containing heme a and the O2 reduction site composed of heme a3 and CuB (Fig. 1A). The reaction requires four electrons and eight protons, of which four are for the reduction of O2 and four are translocated across the membrane. + + O2 + 8HN,chem,pump + 4e− → 2H2 O + 4HP,pump

The four electrons are provided by the sequential oxidation of four cytochrome c red molecules by the one-electron redox center CuA located in the periplasmic domain of Subunit II (SII). Electrons are sequentially shuttled to heme a and then to the reduction site. All of the protons come from the N (electrically negative) side of the membrane and are transferred via two proton-conducting input channels, the D and K channels. 8,20 The D channel transfers two of the chemical protons and all of the pumped protons. It is 25-30 ˚ A long and comprises a continuous hydrogen-bonded network formed by water molecules and conserved polar amino acids linking residue D132 (numbering for the R. sphaeroides enzyme) anchored on the surface of the N side and residue E286 located near heme a and the reduction site. The K channel, which transfers two chemical protons, begins at residue E101 in SII and includes residue K362 near the SI-SII interface. The two channels are not redundant, each providing chemical protons during different steps in the catalytic cycle. Starting with a fully reduced enzyme, the binding of O2 to heme a 3 initiates the O2 reduction chemistry and 4 ACS Paragon Plus Environment

Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

reaction steps associated with proton pumping. The O2 delivery pathway, however, has not been well characterized and is the topic of the current work. A previous MD study from this laboratory with a B-family HCO, cytochrome ba 3 from Thermus thermophilus 1 demonstrated that O2 diffuses from the lipid bilayer to the reduction site through a hydrophobic tunnel that is the same as that defined by Xe-binding sites observed in the crystal structure of the protein (Fig. 1D). 24,25 This tunnel is also apparent in the absence of bound Xe in all of the available crystal structures of cytochrome ba 3 . 26–28 The O2 delivery pathway in cytochrome ba 3 has two entrances enabling O2 access to the protein from the lipid bilayer. Both entrances begin at SI-lipid interfaces (Fig. 1D); one (termed TM1-3) is spanned by TM1, TM2 and TM3 helices, and other (termed TM4-5) is spanned by TM4 and TM5 helices. This structural feature is a “static” tunnel and is referred to as the “X-ray inferred” pathway. An equivalent tunnel with the same two entrances (Fig. 1A-C) is observed in the crystal structures of the A-family oxygen reductases including cytochrome aa 3 from different sources. 13,19,24,29–31 However, in the A-family enzymes, these tunnels appear to be partially occluded (Fig. 1B-C), consistent with the experimental apparent second order rate constant for O2 diffusion from the aqueous solution to the reduction site of cytochrome aa 3 (108 M−1 s−1 ) being 10 times slower than for cytochrome ba 3 . 32,33 The A-family oxygen reductases have a 7-TM helical subunit III (SIII) that is not present in cytochrome ba 3 or other B-family HCOs, and the TM4-5 entrance of the putative O2 pathway in SI is adjacent to SIII. However, there are no data that indicate a functional role of SIII in O2 delivery to the reduction site that might result from its proximity to the tunnel entrance. Although the elimination of SIII has deleterious effects, the A-family enzymes remain functional. 8,34–37 X-ray structures of A-family HCOs without SIII indicate no structural perturbations to SI or SII or to the O2 pathways 14,15,38–40 (Fig. 1B-C). Two MD studies have, however, reached opposite conclusions about the use of the X-ray inferred tunnel by O2 to reach the reduction site in the A-family HCOs. An early (1998) MD study by the Schulten group 41 applied the locally enhanced sampling (LES) technique 42 to simulate O2 diffusion in cytochrome aa 3 oxygen reductases from P. denitrificans (P.d ) and bovine. They observed an O2 molecule initially placed at the reduction site exiting via the entrance of the X-ray inferred tunnel in each case. However, these LES simulations 41 were performed in vacuum and lasted only several picoseconds. Because LES accelerates O2 5 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

diffusion by softening interactions between O2 and its surroundings, it might not correctly capture the dynamics that are functionally relevant for O2 delivery. More recently, an MD study of the A-family R.s. cytochrome aa 3 by Oliveira et al 43 probed the O2 delivery pathway by performing five 100-ns flooding simulations, in which multiple copies of O2 were explicitly included and simulated with the rest of the system (i.e. protein, lipids, water and ions). Unlike LES, flooding simulation is based on equilibrium, conventational simulations with no perturbation and rescaling of interactions between atoms in the system. However, no O2 entry was observed during these simulations. They performed energy estimations using the implicit ligand sampling (ILS) technique, a post-simulation analysis used in probing high affinity O2 sites, 44 and the results suggested two alternate dynamically formed tunnels as preferred O2 pathways over the X-ray inferred tunnel. Experimental mutagenesis data which support the use of the X-ray inferred tunnel as the O2 delivery pathway in the A-family HCOs are based on amino acid residues, such as V279 of P.d. cytochrome aa 3 45 and V287 of E. coli cytochrome bo 3 46 (equivalent to V287 R.s. cytochrome aa 3 ), that are located very close to the reduction site rather than within the tunnel. Also, although an X-ray crystallographic experiment detected two Xe binding sites in the R.s. enzyme, 13 those sites are located near the entrances of the pathway and therefore do not necessarily define a clear pathway for O2 delivery. Since the two previous MD simulations for O2 diffusion to the reduction site of A-family oxygen reductases arrived at different conclusions, the present study re-examines the issue by independently performing both ILS analyses and an extended set of flooding simulations of the diffusion of O2 to the reduction site of the R.s. cytochrome aa 3 . The X-ray model used for the simulations includes SI and SII (Fig. 2). To improve the statistics on the O2 delivery pathways compared to previous studies, twenty independent flooding simulations with different starting points were performed and each was extended to 150 ns. ILS analyses were performed on an independent 200-ns MD simulation. The results clearly show that the only pathway used by O2 to reach the reduction site is the X-ray inferred tunnel, similar to what has been shown for the B-family cytochrome ba 3 from T. thermophilus (T.t. 1 Constrictions along the pathway result in a rate of O2 delivery that is at least 10 times slower than the rate observed 32,33 and calculated 1 for cytochrome ba 3 .

6 ACS Paragon Plus Environment

Page 6 of 33

Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Cytochrome aa3 and its putative O2 delivery pathway

A

90° O 2

O2

B 4-subunit structure C 2-subunit structure D of cytochrome aa3

of cytochrome aa3

O2 pathway in cytochrome ba3

TM6

TM6

TM6 Xe1

TM5

TM5

Xe6

TM5

Xe5 Xe3 Xe4 Xe7

TM2

TM2 TM4 TM1

TM3

TM4 TM1

TM3

TM2 Xe2

TM1

TM4

TM3

Figure 1: Putative O2 delivery pathway in HCOs. A) Location of the pathway in cytochrome aa 3 . The R.s. enzyme comprises four subunits: SI (white), SII (brown), SIII (green) and SIV (not shown). SI and SII are the minimal functional units of HCOs. The O2 reduction site is highlighted using a dashed green box. The X-ray inferred pathway, located in SI, is shown in magenta surface drawn using CAVER 3.0. 47 This structural model is from PDB 1M56. 13 Red arrows on the right panel depict the entrances of the pathway. B-C) Projected from the top view, the putative O2 pathway in cytochrome aa 3 is illutrated by chains of magenta balls drawn by CAVER3.0. Unoccupied spaces within the protein were visualized using the molecular surface representation in VMD with the probe radius of 1.7 ˚ A, which is about the vdW radius of an oxygen atom. Comparison between the crystal structure with 4 subunits (B) and the one with only SI and SI subunits from PDB 2GSM 14 (C) indicates that SIII does not block the TM4-5 entrance. The simulations were performed using the 2-subuit structure. D) The equivalent pathway in cytochrome ba 3 was found to bind Xe (gray balls) by X-ray crystallographic experiments 24,25 and characterized as the O2 delivery pathway in MD simulations. 1 The structural model is from PDB 1XME 27 and the atom coordinates of Xe are from PDB 3BVD. 24 The Xe atoms are numbered according to Luna et al. 24

7 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Methods and Materials System preparation The membrane-embedded model of R.s. cytochrome aa 3 was prepared using the 2.0-˚ A crystal structure (PDB 2GSM) as the protein model. The structure comprises SI (catalytic subunit), which binds low-spin heme a, high-spin heme a 3 and CuB cofactors, SII, which binds the CuA complex, and 282 water molecules. Hydrogen atoms were added using PSFGEN in VMD. 48 Histidine residues except H102, H333, H334, H411, H419 and H421 of SI were in the HSE tautomeric form (N atom of the imidazole ring carrying proton). H102 and H421 are ligated to the Fe atom of heme a, H284, H333 and H334 are ligated to CuB , H419 is ligated to the Fe atom of heme a 3 , and H411 forms a hydrogen bond with the propionate A of heme a 3 . The carboxylate side chain of residue E286, for which pKa has been experimentally estimated to be >9, 49 was assigned to be protonated. The first principal axis of the protein was aligned with the z axis (membrane normal) using the OPM (Orientations of Proteins in Membranes) database. 50,51 The protein was then inserted into a patch of POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine) bilayer. Lipids that overlapped the protein were removed, keeping 132 lipids in the periplasmic and cytoplasmic leaflets. The membrane-embedded enzyme complex was then solvated with water. Water molecules that were in the membrane (-18 < z < 18), except those from the crystal structure , were removed, keeping 20,637 water molecules. Finally, 0.2 M NaCl (78 Na+ and 79 Cl− ions) was added to neutralize and ionize the system, resulting in a fully solvated model of 104,710 atoms. The CHARMM22 force field with φ/ψ corrections 52,53 was used to describe the protein and the heme cofactors, CHARMM36 54 for the lipids, and the TIP3P model 55 for water molecules. The cofactors were treated in the reduced state. The partial charge of 0.383 was used for CuB according to Hofacker and Schulten. 41 The vdW parameters of Cu with  = 0.19 kcal/mol and Rmin = 1.4 ˚ A were from Fuchs et al. 56 The force field parameters for the hydroxyethylfarnesyl side chain, attached to the porphyrin ring of hemes a and a 3 , were not available in the CHARMM force field library; they were constructed by analogy using parameterized alcohol, aldehyde, alkene, and alkane fragments; 57 the complete structure of the hemes with the atomic partial charges is shown in Fig. S1. For the CuA complex, 8 ACS Paragon Plus Environment

Page 8 of 33

Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

the partial charges of its Cu atoms and its ligated amino acids were assigned according to Hofacker and Schulten. 41

Simulation protocols MD simulations performed to prepare the systems consisted of the following steps: (1) 0.5ns melting of lipid tails during which only the lipid tails were allowed to move in order to achieve better packing of lipids around the inserted protein; (2) 0.5-ns simulation with restraints (k = 1 kcal/mol/˚ A2 ) applied to heavy atoms of the protein and cofactors (all lipid atoms and water moving) and with harmonic potentials (k = 0.1 kcal/mol/˚ A2 ) applied to keep water out of the membrane; (3) 0.5-ns simulations with only backbone atoms of the protein and heavy atoms of the cofactors restrained (k = 1 kcal/mol/˚ A2 ); (4) 1-ns simulation with only Cα atoms of the protein and heavy atoms of the cofactors restrained; and (5) 20-ns unrestrained relaxation. Energy minimization (1,000 steps) was performed at the beginning of Steps 1, 2 and 3 using the conjugate gradient algorithm. To maintain the ligation of CuB to H284, H333 and H334 and thereby the structure of the reduction site, the His-CuB , and heme a 3 Fe-CuB connections were described by bonded interactions (k = 200 kcal/mol/˚ A2 for bonds and k = 50 kcal/mol/rad2 for angles). All simulations were performed using NAMD2 58 with a time step of 2 fs and with the periodic boundary condition (PBC). All bonds involving hydrogen atoms were kept rigid using the SHAKE algorithm. 59 To evaluate long-range electrostatic interactions in PBC without truncation, the particle mesh Ewald (PME) method 60 with a grid density of 1/˚ A3 was used. The cutoff for van der Waals interactions was set at 12 ˚ A. All of simulation steps except the melting of lipid tails (Step 1) were performed in a flexible cell, which allows the system to change its dimensions independently as an NPT ensemble. The temperature was maintained at 310 K by Langevin dynamics 61 with a damping coefficient γ of 1 /ps. The Nos´e-Hoover Langevin piston method 61,62 with a piston period of 200 fs was used to maintain the pressure at 1 atm.

9 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Production runs The production runs consist of a 200-ns apo simulation performed in the absence of O2 molecules (used for ILS analysis and starting flooding simulations), and a set of twenty 150-ns flooding simulations in which a large number of O2 molecules were added. The apo simulation started from the 20-ns time point of the relaxation simulation. The set of twenty flooding simulations started from snapshot taken at different time points from the 200-ns apo simulation described above: 4 snapshots (at 25-ns, 100-ns, 150-ns, and 200-ns time points, respectively) were selected from the apo simulation and each one were used to seed 5 indepedent simulations, resulting in an ensemble of twenty simulations. Flooding simulations Flooding simulations were carried out to probe potential O2 delivery pathways and the dynamics associated with the O2 delivery process. To maximize the sampling of O2 delivery pathways within a limited timescale (150 ns), 130 O2 molecules, corresponding to a concentration of 210 mM with respect to the volume of the simulation box, were added to the equilibrated structure of the membrane-embedded cytochrome aa 3 . At the start of the simulation, 70 O2 molecules were placed in the membrane and 60 molecules in the aqueous phas generating a concentration of 185 mM in the aqueous phase with all O2 molecules occupying space outside the protein (Fig. 2). Twenty 150-ns simulations, adding upto a total of 3,000 ns, were carried out to probe for potential O2 pathways. The simulated O2 molecules are described by the standard CHARMM force field. 52 The partial charges for the oxygen atoms of O2 are +0.021 and -0.021, respectively. Its intramolecular interactions are described by the bond distance of 1.23 ˚ A and the spring constant of 600 kcal/mol/˚ A2 . The vdW parameters of the oxygen atoms are with  = -0.12 kcal/mol and Rmin = 1.7 ˚ A. We note that the O2 model with the partial charges of ±0.021 and the one with the partial charges of 0 used in the ILS calculations, in which no electrostatic terms are included, show negligible differences in O2 solvation in the aqueous solution and O2 partitioning in the membrane (Fig. S3-S4 and Table S1). Moreover, it is important to note that although the phrase “O2 ligand for heme” is tagged in the CHARMM topology file containing both O2 and hemes, these charges of ±0.021 are far away from the strongly polarized O2 molecule when ligated to the heme iron. A calculation by Daigle et al 63 showed that oxygen atoms of an O2 molecule ligated to heme are strongly polarized with the partial charges of -0.18 and -0.32, clearly indicating that the slightly charged model (±0.021) used in the flooding simulations is not representing 10 ACS Paragon Plus Environment

Page 10 of 33

Page 11 of 33

a heme-bound O2 . The very small charges of ±0.021 used in the present study (and in other simulation studies 64,65 ) are used to take into account a portion of the polarization that O2 would experience when it approaches to strongly charged portions of the protein. In such cases, the rotation of O2 can orient the small introduced dipole with the surrounding field, thereby giving a small portion of polarization effects. 120

A

O2 in the membrane

100

2

80 NO

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

60

O2 in the aqueous solution

40 20 0

20

40

60

80

100

120

140

time (ns)

B t ~ 20ns

t=0

Figure 2: Flooding simulations of O2 diffusion in the membrane-embedded cytochrome aa 3 . A) An example taken from one of the simulation illustrating the number of O2 molecules (NO2 ) localized in the membrane and in the aqueous solution. At the beginning of the simulation, 70 O2 molecules were placed in the membrane and 60 in the aqueous solution. After reaching the equilibrium, ∼110 O2 molecules localized in the membrane and ∼20 localized in the aqueous solution. B) Snapshots taken from one of the simulations illustrating the starting (left) and equilibrium points. Helices shown in white belong to SI, while ones shown in brown belong to SII. Lipid molecules are shown in orange sticks. To identify whether O2 delivery events took place during the simulations, we used the 6-˚ A distance cutoff from CuB as the criterion. Then, we obtained the identity of the delivered O2 molecules and examined the simulation trajectories whether they were indeed localized in the reduction site.

11 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Implicit ligand sampling (ILS) Complementary to the flooding simulations in which ligand diffusion is explicitly probed, ILS was employed to identify potential regions and pathways for O2 insertion that may not be sufficiently sampled by flooding simulaitons. ILS calculates ligand-interaction energies (Ei ) in any position inside the protein over an ensemble of protein conformations and ligand orientations, 44 which estimate a 3D free energy map of inserting an O2 molecule at position i (∆Gi ). < e−Ei /RT > pi = −RT ln ∆Gi = −RT ln p0 p0 where p0 (in vacuum) = 1 and pi is the probability of inserting an O2 molecule at position i. Following the assumption that small, hydrophobic gases weakly interact with proteins and therefore do not affect the protein structure and dynamics, the 200-ns apo trajectory (10,000 frames) was analyzed for O2 pathways. O2 molecules were sampled in a 55 × 50 × 75 ˚ A3 grid with spacing of 1 ˚ A, covering the entire cytochrome aa 3 enzyme. Ten orientations of O2 were sampled in each subgrid, which contained 3×3×3 interaction sites. The solvation free energy of O2 (∆Gsol ) was used as the reference for calculating the partitioning free energy of O2 (∆Gi,sol ). ∆Gi,sol = ∆Gi − ∆Gsol ˚3 using of NaCl solution ILS and ∆Gsol was independently calculated over a 30×30×30 A free-energy perturbation (FEP). 66 Both techniques yield ∆Gsol values of 2.1 kcal/mol, which is consistent to a previous calculation by Cohen et al. 44

12 ACS Paragon Plus Environment

Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Results and Discussion O2 delivery pathway to the reduction site To unequivocally probe for pathways used by O2 to diffuse to the reduction site of cytochrome aa 3 , flooding simulations with 130 O2 molecules were performed. Twenty simulations were performed using different starting points taken from the equilibrated system of the apo simulation: at t = 25 ns, 100 ns, 150 ns or 200 ns. The O2 molecules were initially placed outside the protein, 70 molecules in the membrane and 60 molecules in the aqueous solution. Each simulation lasted 150 ns, but it took only ∼20 ns to achieve steady distributions of O2 in both the membrane and aqueous phases (Fig. 2), where ∼110 O2 molecules resided in the membrane and ∼20 resided in the aqueous solution. In 15 of 20 simulations, O2 molecules were observed to enter the reduction site; the number of this event ranged from 1 to 7 with an average of 2 events per simulation, corresponding to an average of O2 entry in every 75 ns. Based on the 6-˚ A distance cutoff from the center of Fe of heme a 3 and CuB , the residence times of O2 in the reduction site ranged from ∼1 ns to ∼90 ns (Fig. 3). O2 molecules that reached the reduction site were identified and their diffusion dynamics was visually examined in order to locate the O2 delivery pathways. Their trajectories are shown in Fig. 4; each colored line represents an individual O2 molecule. The results of flooding simulations showed that the delivered O2 molecules enter the reduction site via a pathway, which contains three entry branches (entrances) accessible from the membrane. One of the entrances begins at a lipid-protein interfacial region formed by TM1, TM2 and TM3 (TM1-3) helices, one begins at the TM4-TM5 (TM4-5) interface, and one begins at the TM5-TM6 (TM5-6) interface. In some of the simulations, O2 is observed to enter via all three entrances, while in others, the entries occur via only one or two entrances (Fig. 4). The location of the pathway in cytochrome aa 3 is similar to the one previously defined for O2 delivery in cytochrome ba 3 , which resembles a Y-shaped tunnel (shown in Fig. 1C). The O2 delivery pathway in cytochrome ba 3 1 corresponds to a two-branched hydrophobic tunnel as was also determined to bind Xe by X-ray crystallography. 24,25 The TM1-3 entrance of cytochrome aa 3 is equivalent to Branch A of cytochrome ba 3 , so it is referred to as Branch A. The TM4-5 entrance is equivalent to Branch B, so it is referred to as Branch B. The TM5-6 entrance was not found in cytochrome ba 3 and is denoted as Branch C.

13 ACS Paragon Plus Environment

Biochemistry

Distance between O2 and the reduction site (Å)

simulation time (ns)

A

simulation time (ns)

F

K simulation time (ns)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 33

20

15

10

5

0

0

20

B

15

10

5

0

0

20

C

15

10

5

0

0

D

20

15

10

5

0

0

E

20

20

20

20

20

40

40

40

40

40

60

60

60

60

60

80

80

80

80

80

100

100

100

100

100

120

120

120

120

120

140

140

140

140

140

20

15

10

5

0

0

G

20

15

10

5

0

0

20

H

15

10

5

0

0

20

I

15

10

5

0

0

J

20

20

20

20

40

40

40

40

40

60

60

60

60

60

80

80

80

80

80

100

100

100

100

100

120

120

120

120

120

140

140

140

140

140

20

15

10

5

0

L

20

15

10

5

0

0

M

20

15

10

5

0

0

20

N

15

10

0

5

0

O

20

20

20

20

40

40

40

40

40

60

60

60

60

60

80

80

80

80

80

100

100

100

100

100

120

120

120

120

120

140

140

Bottleneck BNC

Bottleneck

BNC

BNC

5

0

20

15

10

5

0

20

15

10

5

0

140

140

Bottleneck

10

0

20

140

15

0

20

0

20 0

Bottleneck

BNC

Bottleneck

BNC

Figure 3: Access of O2 to the reduction site during 150-ns timespan of the simulations quantified by the distance between an O2 molecule and the binuclear center (denoted as BNC) of the reduction site (O2 -cat) as the function of simulation time. Each panel represents one of the 15 simulations in which O2 delivery took place. In relation to Table 1, panels A-D denote Sims. 1A-D, panels E-H denote Sims. 2A-D, panels I-K denote Sims. 3A and 3D-E, and panels L-O denote Sims. 4A-B and 4D-E. Each line corresponds to the trajectory of an individual O2 molecule diffusing into the reduction site (shaded green), defined by the proximity of < 6 ˚ A from the centers of mass between the Fe atom of heme a3 and CuB . The dashed line in each panel indicates the location of the main bottleneck, corresponding to the kinetic barrier region separating the reduction site from the rest of the O2 delivery pathway, where O2 resided very transiently while diffusing towards the reduction site.

14 ACS Paragon Plus Environment

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

A

B

C

TM6 C

TM4

TM1 TM3

TM3

A

TM3

A

F

TM1

B

TM4

TM3 A

I

TM6

B

TM4

TM1

A

H

TM6

TM2

TM3

A

G

TM5

TM5

TM2

TM4

TM1

TM6 C

TM5

TM2

B

TM4

B TM1

TM6 C

TM5

TM2

E

TM6 C

TM5

TM2

D

TM6

J

TM6

TM6

TM6

C TM5

TM5

TM2

TM4

TM1

TM2

B

TM4

TM1

TM3

B

TM6

TM4

B

TM1

TM3

TM4

N

O

TM6

TM5 B

TM1

TM3

TM4

TM6

TM4

TM1

B

TM6 C

C

TM5

TM2

TM3 A

TM4

TM1

TM3

C

TM2

TM5

TM2

A

TM6

C

B

TM5 B

TM3

M

TM5

TM1

TM2

A

L

TM2

TM4

TM1

TM3 A

K

TM5

TM2

TM5

TM2 TM1

TM3 A

A

TM4 TM3

TM5

TM2 TM1

B TM4 TM3

Figure 4: O2 access routes. The trajectories of the delivered O2 molecules into the reduction site during each of the 15 simulations are illustrated by line traces with different colors. O2 molecules enter the X-ray inferred pathway via three membrane-accessible branches, denoted as “A”, “B” and “C”. Branch A originates at TM1, TM2 and TM3. Branch B originates at TM4 and TM5 which is adjacent to SIII. Branch C, which was previously unknown, originates at TM5 and TM6. In some of the simulations, O2 were found entering via all 3 branches, whereas in others, they entered via only 1 or 2 of the branches were visited by O2 . In relation to Table 1, panels A-D denote Sims. 1A-D, panels E-H denote Sims. 2A-D, panels I-K denote Sims. 3A and 3D-E, and panels L-O denote Sims. 4A-B and 4D-E.

15 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Comparisons of O2 delivery pathways in cytochrome aa 3 and cytochrome ba 3 The slower O2 delivery rate in cytochrome aa 3 (1×108 M−1 s−1 ) compared to the one in cytochrome ba 3 (1×109 M−1 s−1 ) 32,33 is probably due to the presence of diffusions barriers within the protein. However, the inability to observe O2 delivery in some of the flooding simulations makes it difficult to compare our current results with cytochrome aa 3 and our previous ones with cytochrome ba 3 using the results of flooding simulations. Therefore, ILS analyses were performed to calculate thermodynamically favorable O2 regions within cytochrome aa 3 . The results of ILS were used to map out potential O2 delivery pathways, which were then compared to the pathways obtained from the flooding simulations and to the results of cytochrome ba 3 obtained from our previous study. 1 The comparison between cytochrome aa 3 and cytochrome ba 3 is shown in Fig. 5A-B as 3D free energy isosurfaces, in which the colored surfaces represent free energy states of O2 insertion. For cytochrome ba 3 , the entire O2 delivery pathway is energetically favorable as indicated in Fig 5B by the Y-shaped red surface corresponding to ∆G of -3.3 kcal/mol. For cytochrome aa 3 , the pathway identified by the flooding simulations are found to be less favorable for O2 insertion relative to the one in cytochrome ba 3 . The pathway contains several higher energy regions between the entrances and the reduction site illustrated in Fig 5A by the discontinuity of the magenta surfaces corresponding to ∆G of -1.8 kcal/mol forming diffusion barriers. Here, O2 molecules diffusing from any of the three entrances to the reduction site encounter a barrier located in the vicinity of residues F172 and G283. In both cytochrome aa 3 and cytochrome ba 3 , the results of flooding simulations are correlated to the results of ILS calculations. In our previous study of O2 delivery in cytochrome ba 3 , although flooding simulations were performed for only 50 ns using the O2 concentration of 210 mM, we were able to observe O2 delivery on an average of 20 events. 1 For cytochrome aa 3 , the simulations lasted much longer (150 ns), but only two delivery events occurred on the average, which represent a 30-fold decrease compared to cytochrome ba 3 . Fewer O2 delivery events indicate that the pathway in cytochrome aa 3 is less accessible to O2 than the one in cytochrome ba 3 . The study by Oliveira et al 43 with the same R.s. cytochrome aa 3 presented five 100-ns flooding simulations but found no O2 entry. These results are consistent with the present 16 ACS Paragon Plus Environment

Page 16 of 33

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

A

Cytochrome aa3

Cytochrome ba3

B TM6

TM6

C

G283

G232

W172

TM5

Y133

TM5

TM2

TM2

B

B

TM4

TM4

A

A

C TM9

TM8

TM7 TM10 TM6

TM1

TM5

TM2 TM4 TM3

Figure 5: Potential O2 pathways predicted by ILS. A) Equivalent regions to the O2 delivery pathway of cytochrome aa 3 . Red isosurfaces correspond to ∆G of ∼-3.3 kcal/mol. Magenta isosurfaces correspond to ∆G ∼-1.8 kcal/mol. Discontinuous gaps between the isosurfaces indicate unfavorable O2 insertion regions, correlated with hindrances along the pathway. Constricting residues are labeled; E286 (not shown) is located below F282 and G283. B) The entire pathway of cytochrome ba 3 is highly favorable for O2 . The ILS results of cytochrome ba 3 were obtained by analyzing the simulation trajectories for our previous study. 1 C) Other potential O2 delivery pathways in cytochrome aa 3 . Pink surfaces correspond to ∆G of +1 kcal/mol. The O2 pathway, highlighted in dashed red box, is the only region extending all the way from the membrane to the reduction site. The alternate pathways proposed by Oliveira et al 43 start in between TM7 and TM8 and in between SII (brown), TM8 and TM9.

17 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

study insofar we also found free energy barriers that limit the rate of O2 diffusion to the reduction site. Five of the 20 simulations performed in the present study exhibited no O2 entry events while the rest observed only few events (Table 1). Based on their ILS analyses, however, Oliveira et al 43 concluded that O2 prefers to use alternate pathways rather than the X-ray inferred pathway to reach the reduction site 43 although no passage of O2 through such alternative pathways was observed. Different free energy isosurfaces, especially the ones at higher free energy contours (e.g., at ∆G = 1 kcal/mol, shown in Fig. 5C), were used to test for the existence of potential O2 pathways besides the X-ray inferred pathway. However, despite the free energy barriers within the X-ray inferred pathway, the results of flooding simulations show that all of the 39 O2 molecules reaching the reduction site in cytochrome aa 3 during the flooding simulations from the solutions (Table 1) use this X-ray inferred pathway, suggesting that the X-ray inferred pathway is the primary O2 delivery pathway in cytochrome aa 3 .

Delivery rate of O2 To directly connect the results of the present study to the ones experimentally determined from time-resolved absorption spectroscopy, 32,33 we provide an approximation of the second order rate of O2 delivery (kobs ) to cytochrome aa 3 . Since flooding simulations provide dynamic details of O2 diffusion, the obtained data can be used to semi-quantitatively describe steps associated with O2 delivery. The delivery of O2 to the reduction site involves two major steps: 1) the diffusion of O2 from the solution to the entrance(s) of the pathway and 2) the migration of O2 from the entrance(s) to the reduction site. Assuming that the consumption of O2 is 100% efficient once in the reduction site, kobs can be calculated by using the following steady-state kinetics model: k1

k

2 O2 + E E(O2 ) − → E(O2 )cat

k−1

where E(O2 ) is the species in which O2 has arrived at one of the entrances and E(O2 )cat is when O2 is in the reduction site. kobs is defined as: kobs =

[O2 ]k1 k2 . [O2 ]k1 + k−1 + k2

18 ACS Paragon Plus Environment

Page 18 of 33

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

k1 is the rate constant of O2 reaching the entrance(s) of the pathway describing the diffusion step of O2 from the solution to the reduction site. It was calculated as the reciprocal of the product of the time taken to observe the first O2 molecule diffusing from the solutions to the entrance(s) of the pathway during the simulation (tent ) and the aqueous concentration of O2 ([O2 ]). The average tent calculated from all 20 simulations is 2.7 ns and the average [O2 ] in the aqueous solution is 67 mM (Table 1), so k1 is ∼5.5×109 M−1 s−1 . k −1 is the dissociation rate constant of O2 from the entrance of the pathway to the membrane. It is the product of k 1 and the standard concentration (1 M) over the partitioning ratio of O2 at the entrance and the membrane defined as Pent,mem ). Pent,mem is inversely proportional to exponent of the substraction of ∆Gent from ∆Gmem . ∆Gmem is ∆G of O2 in the membrane with respect to the aqueous solution and is -2 to -1.5 kcal/mol. 1,67 ∆Gent is ∆G of O2 at the entrances of the pathway with respect to the aqueous solution. Because the ∆G contour with ∆G = -3.3 kcal/mol is found at all of the entrances according to ILS calculations (Fig. 5A), ∆Gent is approximated to be -3.3 kcal/mol, indicating that O2 is 1.3 to 1.8 kcal/mol or 10-20 folds more favored to partition at the entrance of the delivery pathway than in the membrane. Hence, k−1 is approximated to be ∼2.8-5.5×108 s−1 . k2 is the rate constant of O2 migration from the entrance(s) to the reduction site and is the reciprocal of the time of O2 to diffuse into the reduction site after reaching the entrance(s) of the pathway (tcat−ent ). tcat−ent is obtain by subtracting the time taken to observe the first event of O2 to the catalytic site or tcat , which is ∼57 ns, from tent . Although the average tcat−ent is 54 ns, O2 delivery occurred only in 15 out of 20 simulations or 75% of the total number of simulations. To account for this, the average value of tcat−ent was scaled to 72 ns. Thus, k2 is ∼1.39×108 s−1 . Under physiologically relevant O2 concentrations, it is reasonable to assume that [O2 ] k1  k−1 and k2 , so [O2 ]k1 k2 . kobs ≈ k−1 + k2 This leads to the estimated kobs /[O2 ] of 1.1-1.8×109 M−1 s−1 , which is 8-13 fold slower than the one estimated by our previous MD study with cytochrome ba 3 (15×1010 M−1 s−1 ). 1 Although this estimated kobs /[O2 ] is 10-fold faster than the experimentally determined second-order rate constant of 1×108 M−1 s−1 , 32,33 the experimental rate constant was determined from the time required to convert the enzyme from the fully reduced non-O2 bound state to the ferrous-oxy (O2 bound) state. Hence, the experimental measurement includes the chemical 19 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ligation of O2 to the heme a 3 , the reaction that is beyond to scope of classical MD simulation. The experimental rate constant for O2 to form the heme Fe adduct in T.t. cytochrome ba 3 is 1×109 M−1 s−1 , 32,33 which is 10-fold faster than that determined for R.s. cytochrome aa 3 . The computationally estimated rate constant for O2 delivery to the reduction site of cytochrome aa 3 is 1.1×109 M−1 s−1 , which means that even when the O2 concentration is as low as 10 µM, the rate of O2 diffusion into the reduction site (∼104 s−1 ) will be considerably faster than the rate of O2 catalysis (kcat ∼102 s−1 ). Because of the specific pathway allowing rapid diffusion of O2 from the membrane to the reduction site of the enzyme, the rate of catalysis will not be limited by the ambient concentration of O2 until that concentration in solution is in the low micromolar range.

Constrictions along the O2 delivery pathway in cytochrome aa3 Although the O2 delivery pathway in cytochrome aa 3 contains an additional entrance (Branch C) when compared to cytochrome ba 3 , its presence probably has little or no effect on O2 migration to the reduction site. This is concluded based on previous simulations on cytochrome ba 3 with in silico mutants designed to block the two O2 entrances (Branches A and B in cytochrome ba 3 ). Those mutants, however, did not appreciably impair the passage of O2 . 1 ˚ distant from the reduction site. Three O2 entrances of cytochrome aa 3 merge at 12-15 A At this location, the flooding simulations identified a low O2 sampling region, which can be seen in the left panel of Fig. 6A. The presence of this barrier region is correlated with a broad range of residence times (1-90 ns) for O2 in the reduction site (Fig. 3). This region was also characterized by ILS calculations as a diffusion barrier for O2 . The previous study by Oliveira et al 43 also observed this kinetic barrier and, using a similar approach, calculated the energetic cost of ∼9.5 kcal/mol for O2 to pass through this region, thereby excluding the X-ray inferred pathway as the primary O2 delivery pathway. 43 In the present study, the height of this barrier is suggested to be much smaller since O2 molecules crossed this region in 15 of 20 150-ns simulations. These results support the role of the X-ray-inferred tunnel as the main, if not the only, pathway used for the substrate O2 to diffuse to the reduction site. The region that coincides with the kinetic barrier is surrounded by bulky amino acids M107, W172, F282 and E286 (Fig. 6B, left). This structural feature is common in A-family

20 ACS Paragon Plus Environment

Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Cytochrome aa3

Cytochrome ba3

A TM6

TM6

C

TM5 TM5 TM2

B

TM2

B

TM4

A

TM4

A

B W229

W280 V287

E286

TM5

M107

TM2

V236

F282

G283 W172

TM6

TM6

A77

L246

F108

TM2

V194

G232

I235 Y133 I78

TM5

T231

L200 A149

TM4

TM4

Figure 6: Constricting residues along the O2 delivery pathway of cytochrome aa 3 . A) Occupancy of O2 in cytochrome aa 3 (left) vs. cytochrome ba 3 (right). Red lines are the collections of O2 molecules partitioning in the protein during a flooding simulation (150 ns for aa 3 and 50 ns for ba 3 ). The results of cytochrome ba 3 were taken from the simulation trajectories performed in our previous study. 1 B) Comparison of amino acid residues lining the pathway in cytochrome aa 3 (left) and cytochrome ba 3 (right). Constricting residues, which surround the kinetic barrier region of the O2 delivery pathway in cytochrome aa 3 are labeled in bold; their equivalences in cytochrome aa 3 are also labeled in bold.

21 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

oxygen reductases. Residue M107 is within TM2, W172 is in the loop connecting TM3 and TM4, and F282 and E286 are both within TM6. W172 and F282 have been suggested to restrict the access of O2 . 24 In cytochrome ba 3 , the equivalent region contains amino acids with smaller side chains (Fig. 6B, right), and the barrier height has been estimated to be much less, around 1.5 kcal/mol (Fig. 5C). 1 For example, M107Rs and F282Rs are equivalent to A77 and T231 of cytochrome ba 3 , respectively. W172Rs is equivalent to Y133 of cytochrome ba 3 , and replaceing this tyrosine by a tryptophan by site-directed mutagenesis resulted in a 5-fold decrease in the O2 delivery rate. 68 E286Rs is the terminal amino acid of the D channel, a proton-delivery pathway that is present in the A-family HCOs but is absent in B-family HCOs; 26,69 in T.t cytochrome ba 3 , it is equivalent to I235 (Fig. 6B). The X-ray inferred pathway in cytochrome aa 3 is not a static pathway, and requires dynamics to change from a closed to an open form to allow for O2 to migrate across the barrier (Fig. 7 and Fig. S2). Conformational dynamics of the bulky residues of cytochrome aa 3 were assessed by measuring minimum pairwise distances of F282-E286, F282-F108, W172-E286 and M107-E286 (Fig. 7A). The influences of this dynamics on O2 migration were examined by ILS calculations over selected periods (Fig. 7B). When the side chains of F282, E286 and M107 protrude into the pathway, they approach each other as well as F108 and W172, significantly interfering with O2 migration. For example, during the 171-175 ns period of the 200-ns apo simulation in which all of the four distances came close to 4 ˚ A, the free energy for O2 to migrate through the kinetic barrier can be greater than 4 kcal/mol (Fig. 7B). As indicated in Fig. 7A, the occurrence of these conformational changes, however, is relatively infrequent and might not be observed using shorter simulation times. Finally, it is noted that, in the A-family HCOs, both pumped and chemical protons pass from the D channel to the periplasmic surface and the reduction site via hydrogen-bonded water chains. The terminal region of the O2 delivery pathway is also overlapped to that of the D channel. Since chemical protons and O2 use the same delivery route to the reduction site, the presence of transiently localized water molecules may interfere with the passage of O2 . However, in the simulations performed in the present study as well as those reported in the recent study by Oliveira et al, 43 the region connecting the terminus of the D channel and the reduction site remained dehydrated. Therefore, we can only conclude that the observed restriction of O2 migration is related to the presence of amino acids with bulky side chain. The transient presence of water in this region could also affect the rate of O2 delivery in the A-family HCOs which possess the D channel. 22 ACS Paragon Plus Environment

Page 22 of 33

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

A

F282--E286

F282--F108

W172--E286

M107--E286

B Averaged over 171-175 ns period

C

Averaged over 200 ns

Figure 7: Conformational dynamics of constricting residues and O2 accessibility. A) Minimum distances between lining amino acid residues during the 200-ns apo simulation. Dynamics of the residues affects O2 passage. For example, from 171 to 175 ns (highlighted in red), the side chain of E286 comes in contact with those of W172, M107 and F108, while that of F282 comes very close to that of F108. These close contacts restrict O2 passage. In the crystal structures, N-Cα -Cβ -Cγ of F282 is ∼180◦ which constricted the pathway; its transition to 60◦ increases its contact distances to E286 and F108 from 4 to 6 ˚ A, relieving the hindrance. Residue F108 is located opposite to F282 and in cytochrome ba 3 , it is equivalent to residue I78 (Fig. 6B) B) ∆G isosurface of +1 kcal/mol (pink surfaces) calculated using ILS over the period from 171 to 175 ns. The kinetic barrier is encircled in red, indicating a high energy of O2 insertion and the blockage of O2 passage. C) ∆G isosurface of +1 kcal/mol calculated over the entire 200-ns trajectory. Referring to Fig. 5A, ∆G of the two regions adjacent to the barrier are lower than -3.3 kcal/mol indicated in red surfaces. Hence, the ∆G barrier is greater than 4 kcal/mol when the constricting residues simultaneously form very close contacts. 23 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Sim. 1A 1B 1C 1D 1E 2A 2B 2C 2D 2E 3A 3B 3C 3D 3E 4A 4B 4C 4D 4E Avg.

Page 24 of 33

Table 1: Access of O2 to the reduction site of cytochrome aa 3 No. entry events t cat (ns) Bulk t ent (ns) t cat−ent (ns) 3 4 1 6 0 2 2 2 2 0 2 0 0 1 2 2 7 0 2 2 2

92.6 57.15 36.45 27.02 NA 11.45 39.95 107.38 53.48 NA 91.95 NA NA 78.01 9.97 33.7 58.74 NA 105.03 50.94 57

57 mM 63 mM 60 mM 71 mM 63 mM 70 mM 69 mM 74 mM 70 mM 73 mM 68 mM 74 mM 70 mM 70 mM 60 mM 62 mM 62 mM 64 mM 65 mM 81 mM 67 mM

2.11 6.84 2.59 3.28 3.91 1.93 0.34 0.63 4.38 2.23 1.72 4.04 0.62 4.38 1.94 2.5 0.82 3.85 2.26 2.77 2.7

90.49 50.31 33.86 23.74 NA 9.52 39.61 106.75 49.1 NA 90.23 NA NA 73.63 8.03 31.2 57.92 NA 102.77 48.17 54

t cat is the time taken to observe the first entry event of O2 to the reduction site. t ent is the time taken to observe the first O2 molecule reaching the entrance(s) of the O2 delivery pathway. Sims 1X, 2X, 3X and 4X were prepared using the 25-ns, 100-ns, 150-ns and 200-ns equilibrated complex of the 200-ns apo simulation as the starting model.

24 ACS Paragon Plus Environment

Page 25 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

Conclusion This MD study characterizes the X-ray inferred tunnel as the primary delivery pathway for O2 in R.s. cytochrome aa 3 . This conclusion is opposite to the one made by Oliviera et al, likely related to shorter simulation times used in the previous study 43 due to slow conformational dynamics of amino acids lining the pathway. Because there are multiple constriction regions along the pathway indicated by both flooding simulations and ILS calculations, the simulation time needs to be sufficiently long to observe the migration of O2 to the reduction site of R.s. cytochrome aa 3 . We conclude that both R.s. cytochrome aa 3 and T.t. cytochrome ba 3 use similar pathways to effectively deliver O2 , and this pathway also appears to be conserved in other HCOs.

25 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acknowledgements The authors acknowledge funding supports from the National Institutes of Health (NIH P41GM104601, U01-GM111251 and U54-GM087519 to E.T., and R01-HL16101 to R.B.G.) and the Office of Naval Research (ONR N00014-16-1-2535 to E.T.). P.M. gratefully acknowledges NIH support as a trainee of the Molecular Biophysics Training Program (5T32-GM008276) during his graduate study. Computational resources were provided by XSEDE (XSEDE MCA06N060) and Blue Waters (ACI-1440026).

Supporting Information Available Included in the Supporting Information are the distributions of atomic partial charges assigned to the heme cofactors (Fig. S1), the observed conformations of key amino acid residues at the main barrier of the O2 delivery pathway (Fig. S2), and the comparisons of O2 models used in the flooding simulations and ILS calculations (Fig. S1-S4 and Table 1). The comparisons of the slightly charged O2 model, in which the oxygen atoms contain the partial charges of ±0.021, and apolar model, in which the atoms contain the partial charges of 0, were done by calculating the solvation free energies (Table S1), the interactions between O2 and water molecules depited by the pairwise distribution or g(r) (Fig. S3) and the partitioning profiles of O2 in the membrane and in the aqueous solution (Fig. S4).

26 ACS Paragon Plus Environment

Page 26 of 33

Page 27 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

References (1) Mahinthichaichan, P., Gennis, R., and Tajkhorshid, E. (2016) All the O2 consumed by Thermus thermophilus cytochrome ba 3 is delivered to the active site through a long, open hydrophobic tunnel with entrances within the lipid bilayer. Biochemistry 55, 1265–1278. (2) Pereira, M., Santana, M., and Teixeira, M. (2001) A novel scenario for the evolution of haem-copper oxygen reductases. Biochim. Biophys. Acta – Bioener. 1505, 185–208. (3) Sousa, F. L., Alves, R. J., Ribeiro, M. A., Pereira-Leal, J. B., Teixeira, M., and Pereira, M. M. (2012) The superfamily of heme-copper oxygen reductases: Types and evolutionary considerations. Biochim. Biophys. Acta – Bioener. 1817, 629–637. (4) Morris, R. L., and Schmidt, T. (2013) Shallow breathing: bacteria life at low O2 . Nat. Rev. Microbiol. 11, 205–212. (5) Wikstr¨om, M., Sharma, V., Kaila, V. R., Hosler, J. P., and Hummer, G. (2015) New Perspective on Proton Pumping in Cellular Respiration. Chem. Rev. 115, 2196–2221. (6) Babcock, G., and Wikstr¨om, M. (1992) Oxygen activation and the conservation of energy in cell respiration. Nature 256, 301–309. (7) Ferguson-Miller, S., and Babcock, G. T. (1996) Heme/Copper Terminal Oxidases. Chem. Rev. 96, 2889. (8) Hosler, J. P., Ferguson-Miller, S., and Mills, D. A. (2006) Energy Transduction: Proton Transfer Through the Respiratory Complexes. Annu. Rev. Biochem. 75, 165–187. (9) Kaila, V. R., Verkhovsky, M. I., and Wikstr¨om, M. (2010) Proton-coupled electron transfer in cytochrome oxidase. Chem. Rev. 110, 7062–7082. (10) Brzezinski, P., and Gennis, R. B. (2008) Cytochrome c oxidase: exciting progress and remaining mysteries. J. Bioenerg. Biomembr. 40, 521–531. (11) Rich, P. R., and Marechal, A. (2013) Functions of the hydrophilic channels in protonmotive cytochrome c oxidase. J. R. Soc. Interface 10, 20130183. (12) Ducluzeau, A. L., Schoepp-Cothenet, B., van Lis, R., Baymann, F., Russell, M. J., and Nitschke, W. (2014) The evolution of respiratory O2 /NO reductases: an out-of-thephylogenetic-box perspective. J. R. Soc. Interface 11, 20140196. (13) Svensson-Ek, M., Abramson, J., Larsson, G., T¨ornroth, S., Brzezinski, P., and Iwata, S. (2002) The X-ray Crystal Structures of Wild-type and EQ(I-286) Mutant Cytochrome c Oxidases from Rhodobacter sphaeroides. J. Mol. Biol. 321, 329–339.

27 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(14) Qin, L., Hiser, C., Mulichak, A., Garavito, R. M., and Ferguson-Miller, S. (2006) Identification of conserved lipid/detergent-binding sites in a high-resolution structure of the membrane protein cytochrome c oxidase. Proc. Natl. Acad. Sci. USA 103, 16117–16122. (15) Ostermeier, C., Harrenga, A., Ermler, U., and Michel, H. (1997) Structure at 2.7 ˚ A resolution of the Paracoccus denitrificans two-subunit cytochrome c oxidase complexed with an antibody FV fragment. Proc. Natl. Acad. Sci. USA 94, 10547. (16) Tsukihara, T., Aoyama, H., Yamashita, E., Tomizaki, T., Yamaguchi, H., ShinnzawhaItoh, K., Nakashima, R., Yaono, R., and Yoshikawa, S. (1996) The Whole Structure of the 13-Subunit Oxidized Cytochrome c Oxidase at 2.8 ˚ A. Science 272, 1136–1144. (17) Yoshikawa, S., Shinzawaitoh, K., Nakashima, R., Yaono, R., Yamashita, E., Inoue, N., Yao, M., Fei, M. J., Libeu, C. P., Mizushima, T., Yamaguchi, H., Tomizaki, T., and Tsukihara, T. (1998) Redox-coupled crystal structural changes in bovine heart cytochrome c oxidase. Science 280, 1723–1729. (18) Tsukihara, T., Shimokata, K., Katayama, Y., Shimada, H., Muramoto, K., Aoyama, H., Mochizuki, M., Shizawa-Itoh, K., Yamashita, E., Yao, M., Ishimura, Y., and Yoshikawa, S. (2003) The low-spin heme of cytochrome c oxidase as the driving element of the proton-pumping process. Proc. Natl. Acad. Sci. USA 100, 15304–15309. (19) Shinzawa-Itoh, K., Aoyama, H., Muramoto, K., Terada, H., Kurauchi, T., Tadehara, Y., Yamasaki, A., Sugimura, T., Kurono, S., Tsujimoto, K., Mizushima, T., Yamashita, E., Tsukihara, T., and Yoshikawa, S. (2007) Structures and physiological roles of 13 integral lipids of bovine heart cytochrome c oxidase. EMBO J. 26, 1713–1725. (20) Abramson, J., Svensson-Ek, M., Byrne, B., and Iwata, S. (2001) Structure of cytochrome c oxidase: a comparison of the bacterial and mitochondrial enzymes. Biochim. Biophys. Acta 1544, 1–9. (21) Yoshikawa, S., and Shimada, A. (2015) Reaction Mechanism of Cytochrome c Oxidase. Chem. Rev. 115, 1936–1989. (22) Abramson, J., Riistama, S., Larsson, G., Jasaitis, A., Svensson-Ek, M., Laakkonen, L., Puustinen, A., Iwata, S., and Wikstr¨om, M. (2000) The structure of the ubiquinol oxidase from Escherichia coli and its ubiquinone binding site. Nat. Struct. Biol. 7, 910–917. (23) Musser, S. M., Stowell, M. H., and Chan, S. I. (1993) Comparison of ubiquinol and cytochrome c terminal oxidases: An alternative view. FEBS Lett. 327, 131–136. (24) Luna, V. M., Chen, Y., Fee, J. A., and Stout, C. D. (2008) Crystallographic Studies of Xe and Kr Binding within the Large Internal Cavity of Cytochrome ba 3 from Thermus thermophilus: Structural Analysis and Role of Oxygen Transport Channels in the Heme-Cu Oxidases. Biochemistry 47, 4657–4665. 28 ACS Paragon Plus Environment

Page 28 of 33

Page 29 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

(25) Luna, V. M., Fee, J. A., Deniz, A. A., and Stout, C. D. (2012) Mobility of Xe Atoms within the Oxygen Diffusion Channel of Cytochrome ba 3 Oxidase. Biochemistry 51, 4669–4676. (26) Soulimane, T., Buse, G., Bourenkov, G. P., Bartunik, H. D., Huber, R., and Than, M. E. (2000) Structure and mechanism of the abberrant ba 3 -cytochrome c oxidase from Thermus thermophilus. EMBO J. 19, 1766–1776. (27) Hunsicker-Wang, L., Pacoma, R., Chen, Y., Fee, J., and Stout, C. (2005) A novel cryoprotection scheme for enhancing the diffraction of crystals of recombinants cytochrome ba 3 oxidase from Thermus thermophilus. Acta Cryst. D 61, 340–343. (28) Tiefenbrunn, T., Liu, W., Chen, Y., Katritch, V., Stout, C. D., Fee, J. A., and Cherezov, V. (2011) High Resolution Structure of the ba 3 Cytochrome c Oxidase from Thermus thermophilus in a Lipidic Environment. PLoS One 6, e22348. (29) Iwata, S., Ostermeier, C., Ludwig, B., and Michel, H. (1995) Structure at 2.8 ˚ A Resolution of Cytochrome C Oxidase From Paracoccus denitrificans. Nature 376, 660–669. (30) Tsukihara, T., Aoyama, H., Yamashita, E., Tomizaki, T., Yamaguchi, H., ShinnzawhaItoh, K., Nakashima, R., Yaono, R., and Yoshikawa, S. (1995) Structures of Metal Sites of Oxidized Bovine Heart Cytochrome C Oxidase at 2.8 ˚ A. Science 269, 1069–1074. (31) Lyons, J. A., Aragao, D., Slattery, O., Pisliakov, A. V., Soulimane, T., and Caffrey, M. (2012) Structural insights into electron transfer in caa3 -type cytochrome oxidase. Nature 487, 514–518. (32) Szundi, I., Funatogawa, C., Fee, J. A., Soulimane, T., and Einarsd´ottir, O. (2010) CO impedes superfast O2 binding in ba 3 cytochrome oxidase from Thermus thermophilus. Proc. Natl. Acad. Sci. USA 107, 21010–21015. (33) Einarsd´ottir, O., Funatogawa, C., Soulimane, T., and Szundi, I. (2012) Kinetic studies of the reactions of O2 and NO with reduced Thermus thermophilus ba 3 and bovine aa 3 using photolabile carriers. Biochim. Biophys. Acta – Bioener. 1817, 672–679. (34) Hosler, J. P. (2004) The influence of subunit III of cytochrome c oxidase on the D pathway, the proton exit pathway and mechanism-based inactivation in subunit I. Biochim. Biophys. Acta 1655, 332–339. (35) Varansi, L., and Hosler, J. P. (2012) Subunit III-depleted cytochrome c oxidase provides insight into the process of proton uptake by proteins. Biochim. Biophys. Acta – Bioener. 1817, 545–551. (36) Mills, D. A., and Hosler, J. P. (2005) Slow Proton Transfer through the Pathways for Pumped Protons in Cytochrome c Oxidase Induces Suicide Inactivation of the Enzyme. Biochemistry 44, 4656–4666. 29 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(37) Mills, D. A., Tan, Z., Ferguson-Miller, S., and Hosler, J. (2003) A Role for Subunit III in Proton Uptake into the D Pathway and a Possible Proton Exit Pathway in Rhodobacter sphaeroides Cytochrome c Oxidase. Biochemistry 42, 7410–7417. (38) Qin, L., Sharpe, M. A., Garavito, R. M., and Ferguson-Miller, S. (2007) Conserved lipid-binding sites in membrane proteins: a focus on cytochrome c oxidase. Curr. Opin. Struct. Biol. 17, 444–450. (39) Qin, L., Mills, D. A., Buhrow, L., Hiser, C., and Ferguson-Miller, (2008) A Conserved Steroid Binding Site in Cytochrome c Oxidase. Biochemistry 47, 9931–9933. (40) Qin, L., Liu, J., Mills, D. A., Proshlyakov, D. A., Hiser, C., and Ferguson-Miller, S. (2009) Redox-Dependent Conformational Changes in Cytochrome c Oxidase Suggest a Gating Mechanism for Proton Uptake. Biochemistry 48, 5121–5130. (41) Hofacker, I., and Schulten, K. (1998) Oxygen and Proton Pathways in Cytochrome c Oxidase. Proteins: Struct., Func., Gen. 30, 100–107. (42) Elber, R., and Karplus, M. (1990) Enhanced Sampling in Molecular Dynamics: Use of the Time-Dependent Hartree Approximation for a Simulation of Carbon Monoxide Diffusion through Myoglobin. J. Am. Chem. Soc. 112, 9161–9175. (43) Oliveira, A. S., Damas, J. M., Baptista, A., and Soares, C. M. (2014) Exploring O2 diffusion in A-type cytochrome c oxidases: molecular dynamics simulations uncover two alternative channels towards the binuclear site. PLoS Comput. Biol. 10, e1004010. (44) Cohen, J., Arkhipov, A., Braun, R., and Schulten, K. (2006) Imaging the migration pathways for O2 , CO, NO, and Xe inside myoglobin. Biophys. J. 91, 1844–1857. (45) Riistama, S., Puustinen, A., Verkhovsky, M. I., Morgan, J. E., and Wikstr¨om, M. (2000) Binding of O2 and Its Reduction Are Both Retarded by Replacement of Valine 279 by Isoleucine in Cytochrome c Oxidase from Paracoccus denitrificans. Biochemistry 39, 6365–6372. (46) Riistama, S., Puustinen, A., Garc´ıa-Horsman, A., Iwata, S., Michel, H., and Wikstr¨om, M. (1996) Channelling of dioxygen into the respiratory enzyme. Biochim. Biophys. Acta 1275, 1–4. (47) Chovancova, E., Pavelka, A., Benes, P., Strnad, O., Brezovsky, J., Kozlikova, B., Gora, A., Sustr, V., Klvana, M., Medek, P., Biedermannova, L., Sochor, J., and Damborsky, J. (2012) CAVER 3.0: A Tool for the Analysis of Transport Pathways in Dynamic Protein Structures. PLoS Comput. Biol. 8, e1002708. (48) Humphrey, W., Dalke, A., and Schulten, K. (1996) VMD – Visual Molecular Dynamics. J. Mol. Graphics 14, 33–38.

30 ACS Paragon Plus Environment

Page 30 of 33

Page 31 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

(49) Namslauer, A., Aagaard, A., Katsonouri, A., and Brzezinski, P. (2003) Intromolecular Proton-Transfer Reactions in a Membrane-Bound Proton Pump: The Effect of pH on the Peroxy to Ferryl Transition in Cytochrome c Oxidase. Biochemistry 42, 1488–1498. (50) Lomize, M. A., Lomize, A. L., Pogozheva, L. D., and Mosberg, H. I. (2006) OPM: Orientations of Proteins in Membranes database. Bioinformatics 22, 623–625. (51) Lomize, M. A., Pogozheva, I. D., Joo, H., Mosberg, H. I., and Lomize, A. L. (2012) OPM database and PPM web server: resources for positioning of proteins in membranes. Nucleic Acids Res. 40, D370–376. (52) MacKerell, Jr., A. D. et al. (1998) All-atom empirical potential for molecular modeling and dynamics studies of proteins. J. Phys. Chem. B 102, 3586–3616. (53) MacKerell, Jr., A. D., Feig, M., and Brooks, III, C. L. (2004) Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comp. Chem. 25, 1400–1415. (54) Klauda, J. B., Venable, R. M., Freites, J. A., O’Connor, J. W., Tobias, D. J., Mondragon-Ramirez, C., Vorobyov, I., MacKerell Jr., A. D., and Pastor, R. W. (2010) Update of the CHARMM all-atom additive force field for lipids: Validation on six lipid types. J. Phys. Chem. B 114, 7830–7843. (55) Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983) Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. (56) Fuchs, J.-F., Nedev, H., Poger, D., Ferrand, M., Brenner, V., Dognon, J.-P., and Crouzy, S. (2006) New model potentials for sulfur-copper(I) and sulfur-mercury(II) interactions in proteins: From ab initio to molecular dynamics. J. Comp. Chem. 27, 837–856. (57) Vanommeslaeghe, K., Hatcher, E., Acharya, C., Kundu, S., Zhong, S., Shim, J., Darian, E., Guvench, O., Lopes, P., Vorobyov, I., and MacKerell, Jr., A. D. (2010) CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comp. Chem. 31, 671–690. (58) Phillips, J. C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R. D., Kale, L., and Schulten, K. (2005) Scalable Molecular Dynamics with NAMD. J. Comp. Chem. 26, 1781–1802. (59) Ryckaert, J.-P., Ciccotti, G., and Berendsen, H. J. C. (1977) Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comp. Phys. 23, 327–341. 31 ACS Paragon Plus Environment

Biochemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(60) Darden, T., York, D., and Pedersen, L. G. (1993) Particle mesh Ewald: An N ·log(N ) method for Ewald sums in large systems. J. Chem. Phys. 98, 10089–10092. (61) Martyna, G. J., Tobias, D. J., and Klein, M. L. (1994) Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 101, 4177–4189. (62) Feller, S. E., Zhang, Y., Pastor, R. W., and Brooks, B. R. (1995) Constant pressure molecular dynamics simulation: The Langevin piston method. J. Chem. Phys. 103, 4613–4621. (63) Daigle, R., Guertin, M., and Lague, P. (2009) Structural characterization of the tunnels of Mycobacterium tuberculosis truncated hemoglobin N from molecular dynamics simulations. Proteins: Struct., Func., Bioinf. 75, 735–747. (64) Wang, Y., Cohen, J., Boron, W. F., Schulten, K., and Tajkhorshid, E. (2007) Exploring Gas Permeability of Cellular Membranes and Membrane Channels with Molecular Dynamics. J. Struct. Biol. In press. (65) Wang, Y., and Tajkhorshid, E. (2010) Nitric oxide conduction by the brain aquaporin AQP4. Proteins: Struct., Func., Bioinf. 78, 661–670. (66) Frenkel, D., and Smit, B. Understanding Molecular Simulation From Algorithms to Applications; Academic Press: California, 2002. (67) Wang, Y., and Tajkhorshid, E. (2007) Molecular Mechanisms of Conduction and Selectivity in Aquaporin Water Channels. J. Nutr. 137, 1509S–1515S. (68) McDonald, W., Funatogawa, C., Li, Y., Szundi, I., Fee, Y. C. J. A., Stout, C. D., and Einarsd´ottir, O. (2013) Ligand access to the active site in Thermus thermophilus ba 3 and bovine heart aa 3 cytochrome oxidases. Biochemistry 52, 640–652. (69) Chang, H. Y., Hemp, J., Chen, Y., Fee, J. A., and Gennis, R. B. (2009) The cytochrome ba 3 oxygen reductase from Thermus thermophilus uses a single input channel for proton delivery to the active site and for proton pumping. Proc. Natl. Acad. Sci. USA 106, 16169–16173.

32 ACS Paragon Plus Environment

Page 32 of 33

Page 33 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Biochemistry

For Table of Contents use only

33 ACS Paragon Plus Environment