Filename: b313220j

Back | Show only these items: | Modeling protonated water networks in bacteriorhodopsin

1
Modeling protonated water networks in bacteriorhodopsin

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj1

A model is constructed to investigate in atomistic detail the structure and dynamics of protonated water networks on the extracellular side of the transmembrane proton pump bacteriorhodopsin.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj2

The protein is embedded in a solvated lipid bilayer membrane described by established force fields.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod1

Most importantly, the protonated water network is treated in the framework of a mixed density functional electronic structure theory/molecular mechanics approach.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod2

This QM/MM Car–Parrinello molecular dynamics approach employed allows for stable dynamics on a picosecond time scale.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod3

The structural building process, force field parameterizations and subsequent equilibration of the total system consisting of the protein, the lipid membrane and the hydration layer is described in detail followed by a QM/MM simulation of the protonated water network.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj3

It is found that hydrogen-bonded networks around both H₃O⁺ and H₅O₂⁺ cores can be stabilized in the protein matrix, leading to so-called Eigen and solvated Zundel complexes, H₃O⁺· (H₂O)₃ and H₅O₂⁺·(H₂O)₄, respectively.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res1

It turns out that both complexes behave qualitatively similarly to the gas phase, implying that the H₅O₂⁺ core displays an essentially symmetric hydrogen bond with the excess proton being equally shared between two water molecules.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res2

The dynamics of this hydrogen bond is found to be complex featuring slow large-amplitude motion of the central proton as well as complex dynamics of the protonated water cluster hydrogen bonds formed with the protein matrix.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res3

These findings are consistent with the proposal that an essential component of the so-called proton release group “XH” could consist of a protonated water network stabilized by polar aminoacids.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con1

Motivation: Vectorial proton transport through membranes

Bacteriorhodopsin (BR) is a photoactive transmembrane protein residing in the purple membrane of the extreme halophile Halobacterium salinarium.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac1

Being a light-driven proton pump, BR establishes actively vectorial proton transport from the cytoplasmatic medium to the extracellular side.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac2

The protein consists of seven α-helices, which contain the co-factor retinal acting as the chromophore during the photo-excitation process.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac3

Upon photon absorption, the all-trans 15-syn retinal, which is covalently bound to the Lys216 residue through a protonated Schiff base linkage, undergoes an ultrafast isomerization into a 13-cis 15-anti conformation.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac4

This photo-induced conformational change in the chromophore triggers a cascade of conformational changes in the protein, the so-called BR photocycle, during which a proton is translocated vectorially from cytoplasm to the extracellular side.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac4

The photocycle consists of various major intermediates, labeled K, L, M, N, and O, which are characterized by their absorption maxima and can be probed individually by a variety of experimental techniques.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac5

For recent articles that comprehensively summarize the tremendous progress made by experiments on understanding the delicate structure/function relationship evident in BR see e.g^{refs. 1–6}. and references therein; a very concise exposition of the essentials is given in .^{ref. 7}

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac6

Although the molecular structure of BR is in parts known even on a time-resolved level in terms of the above-mentioned intermediates, the detailed microscopic nature of the proton translocation process still remains unclear.

Type: Motivation | Advantage: None | Novelty: None | ConceptID: Mot1

However, recent experiments and computer simulation studies suggest that BR uses quasi one-dimensional hydrogen-bonded water chains, so-called water wires, during the proton relay process, as suggested many years ago for similar situations based on theoretical considerations.^8,9

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac7

There is also now clear evidence, e.g. from high resolution structures of various intermediates, that water molecules are actively involved in the proton translocation process in BR as opposed to just being “solvent”.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac8

A great deal of further insight could be obtained from time-resolved Fourier-transform infrared spectroscopy (FTIR) which allows for the experimental monitoring of the formation of various protonated species throughout the photocycle.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac9

Despite these advances in our understanding mainly driven by experiment, the detailed role and degree of involvement of these water molecules has still not been resolved.

Type: Motivation | Advantage: None | Novelty: None | ConceptID: Mot2

These are exactly the kind of problems where modern ab initio molecular dynamics (MD) simulation techniques can step into place in order to gain further insight into the microscopic details underlying the proton translocation process in BR.

Type: Motivation | Advantage: None | Novelty: None | ConceptID: Mot3

For a realistic modeling of the overall proton translocation, which very likely has to proceed along a multitude of consecutive hydrogen-bond (H-bond) formation, breaking, and proton transfer steps, an accurate description of the underlying reactive potential energy surface has to be employed and accompanied by a dynamical treatment of the atomic movements.

Type: Motivation | Advantage: None | Novelty: None | ConceptID: Mot4

These two prerequisites can be achieved by treating the electronic degrees of freedom explicitly relying on electronic structure calculations, while evolving concurrently the nuclear degrees of freedom using forces obtained “on the fly”.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met1

Such ab initio, and in particular Car–Parrinello molecular dynamics approaches¹⁰ have proven to be a key tool in studying a large number of chemical reactions ranging from proton transfer in liquid water, systems under extreme conditions, or reactions at surfaces.¹¹

Type: Method | Advantage: Yes | Novelty: Old | ConceptID: Met1

Unfortunately, the size of the systems amenable to this elegant approach is limited by the currently available computer resources and will not exceed 10³ atoms in the near future when it comes to molecular dynamics simulations.

Type: Method | Advantage: No | Novelty: Old | ConceptID: Met1

Therefore, in order to be able to study such complex biomolecular systems on a microscopic level a so-called quantum-mechanical/molecular-mechanical (QM/MM) approximation has to be employed; note that the system introduced below, i.e. one BR molecule in a double-layer lipid membrane which is solvated by explicit water, includes close to 30 000 atoms in total.

Type: Method | Advantage: None | Novelty: New | ConceptID: Met2

Within the QM/MM framework, the total system is partitioned into a fairly small subsystem, the chemical active region, being treated in a complete electronic structure calculation approach, and the “rest”, which is described by established empirical force fields.

Type: Method | Advantage: None | Novelty: New | ConceptID: Met3

This partitioning scheme has the advantage that the computational effort is concentrated on the part of the system where it is mandatory to perform an electronic structure calculation, whereas the effects of the surrounding classical MM part is taken into account by more approximate means.

Type: Method | Advantage: Yes | Novelty: New | ConceptID: Met3

Due to the saving in computational time the QM/MM approach opens the doorway to study interesting reactive chemical processes in extended biorelevant systems.

Type: Method | Advantage: Yes | Novelty: New | ConceptID: Met2

In the long run our work aims at helping to elucidate the various stages of the proton transfer processes that occur during the BR photocycle after retinal isomerization.

Type: Goal | Advantage: None | Novelty: None | ConceptID: Goa1

To this end, a recently developed QM/MM methodology¹² based on a Car–Parrinello type approach^10,11 for the QM part in the framework of density functional theory (DFT) coupled to a MM part treated within the Gromos force field and software^13–15 is used; see ^{ref. 16} for an overview of this novel methodology including several ongoing applications.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met4

The first goal certainly is to build a detailed microscopic all-MM model that resembles the situation found for BR in the experiments as closely as possible, see Fig. 1.

Type: Goal | Advantage: None | Novelty: None | ConceptID: Goa2

This implies that a BR monomer¹⁷ has to be embedded in a membrane, which itself needs to be hydrated by water.

Type: Method | Advantage: None | Novelty: New | ConceptID: Met5

Secondly, the methodology should provide a potential for the proton transfer processes taking place in the QM part with sufficient accuracy and reliability.

Type: Goal | Advantage: None | Novelty: None | ConceptID: Goa3

Finally, the employed computational framework should result in stable classical dynamics for all nuclei in order to be able to extract meaningful dynamical information, in particular for the QM subsystem.

Type: Goal | Advantage: None | Novelty: None | ConceptID: Goa4

This paper reports our initial results on the structural model building of the total system, the partial force field parametrization done, the dynamical relaxation of the total QM/MM system, and structural as well as dynamical data concerning the embedded protonated water network in two different protonation states.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj4

Hereby we focus on the so-called Eigen and (solvated) Zundel complexes, protonated water clusters based on four and two (six) water molecules, respectively.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj5

They are located in a hydrophilic pocket on the extracellular side of BR, see Fig. 2, which is known to contain about three to four crystallographically resolved water molecules in the proximity of residues Glu194 and Glu204, see ^{ref. 18 and 19}.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac10

In particular, close to the L → M transition during the photocycle a proton is released from the so-called “proton release group” (PRG or “XH”) to the extracellular surface of the membrane.²⁰

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac11

Most interestingly, FTIR experiments²¹ provide evidence that the PRG is unlikely to be Glu204 as suggested earlier, rather IR continuum absorbance changes indicate intramolecular proton transfer to the extracellular side via a H-bonded network, and that this network is stabilized by Glu204 and Arg82.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac12

In addition, theoretical studies²² of pK_a values suggest that the proton is stored in terms of a protonated water cluster, in particular H₅O₂⁺, rather than at the close by glutamic acid residues Glu194 and Glu204.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac13

Thus, it is plausible that this H-bonded water network in the hydrophilic pocket “stores” the proton and as such might constitute the PRG, possibly together with not yet identified neighboring groups.

Type: Hypothesis | Advantage: None | Novelty: None | ConceptID: Hyp1

We report herein results from fully atomistic QM/MM molecular dynamics simulations concerning both the structure and dynamics of different protonated water clusters in this pocket of the BR transmembrane proton pump.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj6

Building, relaxing and simulating the structural model

Starting structures

As described in Section 1 we are interested in the proton release after the photo-induced isomerization as is happening close to the L → M stage of the photocycle.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj7

To this end we have assembled a model starting from a monomeric protein unit¹⁷ from the crystal structure reported in ^{ref. 23} as obtained from the Brookhaven Protein Data Bank entry 1E0P (chain B therein).

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod4

The partially resolved residues (Met163, Arg227, Glu232) were constructed by finding approximate starting locations for the missing atoms based upon known bond lengths and angles for those amino acids.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod5

The missing residues (#1–4 and #233–248), not present in the crystal structure due to the strong disorder in the C and N termini, were neglected with the exception of residues #2–4 which were approximated from the PDB data bank entry 1QM8 representing the structure of wild type BR in the ground state at 100 K. As this structure only contains four water oxygen atom positions with three of them being on the extracellular part of BR, the coordinates of 21 already known internal water molecules have been taken from the wild type ground state reported in crystal structure 1E0P (chain A).

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod6

In order to mimic the cell wall in which the BR is embedded we have chosen to employ a POPC membrane (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine).

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod7

To construct the lipid bilayer we employed a relaxed starting structure as reported in the literature.^24,25

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod8

This structure was modified to yield a cubic supercell containing 125 POPC molecules in total with a central cavity large enough to accommodate the BR protein assembled as described above.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod9

The entire lipid bilayer including the BR pump was then solvated by 6504 randomly oriented water molecules, represented in terms of three “explicit atoms”, to provide the initial starting configuration for classical molecular dynamics simulations.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod10

The position of the proton within the protein was chosen to be in the hydrophilic pocket near Glu204 following the arguments outlined in Section 1.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod11

In order to allow for long classical relaxation and equilibration runs before switching to the QM/MM interface, where the timescale is severely limited by the Car–Parrinello QM component, the crystallographic water molecules Wat403 and Wat404 were replaced by a H₂O₅⁺ model complex as suggested earlier²².

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod12

Force field parameters

The retinal chromophore is covalently bound to the Lys216 residue via a N = C bond resulting from the formation (and subsequent deprotonation) of a Schiff base between C–OH of retinal and the –NH₂ group of Lysine.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac14

In the part of the photocycle of current interest the retinal nitrogen atom of the Schiff base is not protonated as it is in the ground state of BR.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac15

This necessitated the parameterization of all-atom Charmm^26,27 as well as Gromos¹³ force fields, to be used for the classical simulations as outlined in Section 2.3, in order to represent a deprotonated “Lysine plus retinal” combined residue which we denote here by Lyr, see e.g^{ref. 28}. for a parameterization of the usually needed protonated Schiff base.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod13

For both the Charmm and Gromos force fields the required charges for conjugated C and N atoms and the N = C stretching force constant were obtained from isolated system DFT electronic structure calculations employing the BLYP gradient corrected density functional^29,30 with the standard 6-31G** basis set using the Gaussian98³¹ quantum chemistry package.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod14

The structure of the gas phase lysine/retinal adduct “Lyr” was hereby fully relaxed and the N = C force constant was taken from the diagonal element of the force constant matrix as projected onto this local stretching mode.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod15

Atomic charges were obtained by performing an Natural Atomic Orbital (NAO) analysis.³²

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod16

The NAO approach gives atomic charges that are known to exibit only a weak basis set dependency, unlike the more ubiquitously employed Mulliken scheme, and has the useful feature that for a given charge group as assigned in either of the two force fields one obtains the same net charge.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac16

Note that different, though probably very similiar, values for the partial atomic charges could have been obtained from other well-established charge fitting approaches, like e.g. ESP fitting³³ or Bader analysis.³⁴

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac17

In order to bring the computed NAO charges in accord with the existing fragments of the Charmm force field they have been scaled by a common factor of 0.7, see Table 1.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod17

For the Gromos force field all other parameters required to construct the Lyr fragment were taken from the existing parameters for the all-trans retinal.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod18

In the case of the Charmm force field torsional motion about the conjugated C = C and C = N bonds also required parameterization.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod19

The barrier heights required in the Charmm force field were approximated from the force constant obtained as a diagonal element from the force constant matrix for a given torsional mode assuming a perfect cosinusoidal torsional potential as required for the Charmm force field.^26,27

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod19

The derived force field parameters are listed in Table 2 where labels follow the nomenclature outlined in Fig. 3.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod20

Non-bonded Lennard-Jones (LJ) parameters for conjugated C atoms of retinal in Charmm were taken from unsaturated C atoms such as found in benzene or ethylene.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod21

Both force fields provided reasonable optimized equilibrium structures, with average absolute deviations in the conjugated C = N and C = C bond lengths of about 0.01–0.03 Å in comparison to the DFT reference data and ca. 0.03–0.04 Å with respect to the experimental starting structure, see Table 3.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs1

As a more stringent test of the Charmm force field the neutral Lyr molecule was run as an isolated molecule for 1 ns by molecular dynamics in the gas phase.

Type: Experiment | Advantage: None | Novelty: None | ConceptID: Exp1

During this time no cis/trans isomerization about the π-conjugated retinal backbone was observed indicating that the derived force field provides a stable local minimum of the deprotonated 13-cis 15-anti retinal.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res4

For the classical relaxation and equilibration runs the H₅O₂⁺ Zundel cation was defined as a new “residue” (keeping its optimized structure frozen) in the Charmm and Gromos force fields; note that this is only needed for relaxation purposes and that in the QM/MM simulations this species and larger protonated water clusters are treated quantum-mechanically within DFT.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod22

Thus, the derived force field need only qualitatively reflect the proposed species.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod22

As such LJ parameters were taken to correspond with the appropriate TIP3P water model in Charmm or the SPC water as used in the Gromos force field.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod23

The atomic charges on the Zundel ion (O: −0.82, H_terminal: 0.53, and H_bridge: 0.52) and the O–H_bridge stretching and O–H_bridge–O bending force constants (147 kcal mol⁻¹ Å⁻² and 32 kcal mol⁻¹ rad⁻², respectively), have been obtained in the same fashion as for the retinal/lysine complex Lyr from DFT calculations of isolated H₅O₂⁺.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod24

The position of this ion within the protein was chosen to coincide with crystallographic waters Wat403 and Wat404 which were replaced with the parameterized Zundel complex H₅O₂⁺.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod25

The starting configuration was first relaxed via local geometry optimization with all the protein atoms frozen after which the Zundel oxygen atoms were kept fixed in space for the remainder of all the classical simulations.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod26

Simulation of the POPC lipids within the Gromos force field was made possible by the availability of the set of LJ parameters within the Gromacs simulation package^35,36.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod27

Classical relaxation and equilibration

The equilibration of a complex and inhomogeneous system containing about 30 000 atoms needs to be done in several states which we outline in the following protocol.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj8

Our classical MD run, which targets to generate starting phase space points for the subsequent QM/MM studies, is comprised of two distinct phases.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met6

Phase A involves the use of the Charmm all-atom force field²⁷ and the NAMD simulation package^37,38 in order to prepare efficiently a stable structural model at ambient conditions, including the optimization of the lattice constant of the cubic supercell.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met6

Phase B involves the use of the Gromos force field and software^13–15 for preparing equilibrated starting configurations being used in the QM/MM production runs.¹²

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met6

In phase A (Charmm force field and NAMD code), the following steps have been performed:

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod28

(1) The BR structure (as assembled above) was first relaxed via geometry optimization in order to relax the added fragments, while keeping the crystallographically known heavy atoms including the internal water molecules inside the pump fixed.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod29

(2) Next, this partially relaxed BR complex was placed into the membrane cavity (as prepared above) and the resulting BR/membrane assembly was relaxed further by geometry optimization keeping the protein backbone, the internal water molecules, and the phospholipid headgroup atoms fixed.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod30

(3) This relaxed BR/membrane complex was than solvated by 6504 TIP3P water molecules and only these water molecules were allowed to relax (keeping the BR/membrane subsystem fixed including the Zundel complex and the internal water), followed by a 50 ps NVT MD simulation.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod31

The MD was accomplished with a time step of 0.15 fs, Berendsen temperature thermostatting at a temperature of 300 K, and the Shake algorithm for O–H bonds of water.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod32

(4) Then, the protein side chains, lipid aliphatic tails, internal waters, H-atoms of the Zundel complex, and finally the phospholipid head groups were subsequently released in a stepwise fashion over a period of about 300 ps of simulation time and the system was allowed to equilibrate at constant volume for a further 200 ps, while preserving the other constraints.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod33

(5) After a reasonably equilibrated low energy state had been reached within these constant volume simulations, a MD run within the NpT ensemble employing a variable cubic cell and the Berendsen pressure coupling scheme was performed in order to optimize the lattice constant of the supercell.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod34

The average pressure of the system was slowly equilibrated at 1 bar over a period of 500 ps.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod35

The resulting structure of phase A is the initial configuration for phase B, making use of the Gromos force field and program as required for the QM/MM interface, which proceeded as in the final two steps as in phase A; these simulations were run with a conservative time step of 0.1 fs, the Shake algorithm for all “heavy atom to H-bonds”, and Berendsen thermostats at 300 K and barostats at 300 K and 1 bar.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod36

The system was first equilibrated with constant volume over a period of 100 ps followed by a further 100 ps constant pressure simulation.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod37

Finally the structure was further equilibrated in the NVT ensemble using the optimized fixed cubic cell with axis length of 72.35 Å (corresponding to the average lattice constant at p ≈ 1 bar) for an additional 1 ns prior to collecting equilibrium all-MM classical data for about 3 ns.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod38

This trajectory is needed for an analysis of the long-term stability of the classical biomatrix around the QM part as reported in Section 3.2.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod38

In this phase B, in addition to fixing the BR protein backbone atoms and the Zundel ion oxygen atoms (formerly being water Wat403 and Wat404) the oxygen atom of the crystallographically reported²³ water molecule Wat405, being also located in the hydrophilic pocket close to Glu204, was fixed.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod39

This water molecule remained in close proximity with the Zundel ion throughout phase A and was therefore fixed in phase B in order to provide up to three water molecules that will allow the formation of a rather extended protonated water network once the QM/MM interface is switched on.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod40

As an aside it is noted that yet another internal water molecule (denoted²³) also migrated into this pocket during phase B of our equilibration and continued to reside there throughout the run.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs2

This is in agreement with the propositions advanced in ^{ref. 21, 22 and 39} that the cavity should contain about three to four water molecules.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res5

100

It is most interesting to note that upon switching to the quantum treatment of these water molecules including the proton, this network spontaneously rearranged and formed an Eigen cation H₃O⁺·(H₂O)₃ with several H-bonds to close by residues from the pocket, see Section 3.3 for a detailed analysis.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs3

QM/MM interfacing

101

The equilibrated final structural model from above was used as input for the QM/MM simulations runs relying on the recently introduced approach by Laio, VandeVondele, and Röthlisberger.¹²

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj9

102

This approach is the basis for an efficient and dynamically stable way of performing mixed QM/MM Car–Parrinello molecular dynamics in biochemical systems.

Type: Object | Advantage: Yes | Novelty: New | ConceptID: Obj9

103

In a nutshell, it is based on a DFT description of the electronic structure of the QM subsystem in terms of a plane wave/pseudopotential approach using the well-established CPMD code.^40,11

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met7

104

The QM/MM interface¹² combines this code to the Gromos software^14,15 that represents the MM part in terms of the Gromos force field.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met7

105

The interface approach uses an efficient Hamiltonian electrostatic coupling scheme to include the electrostatic effects due to the classical environment.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met8

106

This is accomplished by a hierarchical multipole method which is used to describe the interaction of the QM and MM regions (in addition to keeping the standard van der Waals interaction terms between the QM and MM atoms or sites as parameterized in the MM force field used).

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met8

107

The QM/MM electrostatic energy term is hereby decomposed into a short-range part in which the interactions are accurately described by directly evaluating the Coulombic terms resulting from the electronic charge distribution located on all real space grid points with all MM atoms.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met9

108

The long-range electrostatic part (for distance larger then about 10 Å in the current simulation) is treated by a multipole expansion of the QM electronic charge distribution.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met9

109

Furthermore, in order to prevent the so-called electron spill-out problem, which is especially severe if very flexible basis sets such as plane waves are used, the bare Coulomb potential is renormalized in terms of a Padé representation such that a finite value is approached as r → 0.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met10

110

The proper 1/r behavior is recovered far away from the nuclei, rr_c, where r_c is close to the covalent radius of the corresponding atom.

Type: Method | Advantage: None | Novelty: Old | ConceptID: Met10

111

This modification prevents that the MM atoms, which are lacking the Pauli repulsion from the missing electron cloud, act as traps for electron density leading to an artificial over-polarization and electron flow of the electron density.

Type: Method | Advantage: Yes | Novelty: Old | ConceptID: Met10

112

For further details we refer to the original literature.¹²

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac18

113

In particular, the QM subsystem is treated here using Troullier–Martins pseudopotentials⁴¹ and the wave function is expanded in a plane wave basis set with 70 Ry cutoff inside a cubic box with length 30 au subject to non-periodic boundary conditions.¹¹

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod41

114

The BLYP density functional,^29,30 which has proven to give reliable structure, dynamics and energetics for a wide range H-bonded system, is employed.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod42

115

A conservative time step of 0.05 fs was chosen and a fictitious electron mass of 500 au was employed in all QM/MM simulations, whereas a time step of 0.10 fs and the same electron mass was found to guarantee stable numerical integration of the Car–Parrinello equations of motion for pure QM trajectories; note that the proper hydrogen mass is used instead of substituting by deuterium as often done.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod43

116

All dynamical DFT simulations were performed in the NVT ensemble employing Nosé–Hoover chain thermostats¹¹ with a target temperature of 300 K. Despite the fact that these simulations were started using the carefully pre-equilibrated configurations from above, the QM subsystem still needed to be relaxed initially in order to adjust as smoothly as possible to switching from the all-MM to the QM/MM interface treatment.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod44

117

The cubic box containing all 28 365 (QM and MM) atoms (not counting the MM hydrogens treated in the united-atom approximation) was 72.35 Å, see Fig. 1, and trajectories in the typical time-regime accessible to Car–Parrinello dynamics,¹¹i.e. ∼10 ps, were generated within the QM/MM framework.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs4

118

After switching from all-MM to QM/MM simulations no constraints were applied in the simulation of the Eigen complex, whereas some distance constraints had to be imposed in order to stabilize the solvated Zundel complex, see Section 3.3 for more details.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod45

119

Finally, all-QM reference calculations of protonated water clusters in the gas phase at 300 K were performed to generate consistent data sets for the sake of comparison.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod46

Simulation of protonated water networks in BR

120

Protonated water clusters H⁺·(H₂O)_n are quite stable molecular ions in the gas phase⁴² In particular, it is known that the so-called “Zundel” and “Eigen complexes”, obtained by protonating clusters consisting of two and four water molecules, respectively, are two archetypal “solvating schemes”.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac19

121

In the Zundel complex, H₅O₂⁺, the most probable position of the proton is midway between the two solvating water molecules, schematically denoted as [H₂O⋯H⋯OH₂]⁺.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac20

122

This situation is very different for the Eigen complex, where the proton and one water molecule form a hydronium core H₃O⁺, which itself is microsolvated by the three more molecules, i.e. H₃O⁺·(H₂O)₃.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac21

123

These two basic structures can be thought of as building blocks of ever larger clusters.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac22

124

Of particular importance in the following discussion is what we call the “solvated Zundel complex” obtained by saturating the two terminal water molecules with four additional H-bonds leading to H₅O₂⁺·(H₂O)₄.

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac23

125

Standard structural parameters to investigate H-bonds in general, and in particular in protonated complexes, are the so-called asymmetric stretch coordinate δ = R_OH–R_HO, defined in terms of the covalent O–H and H–bonding H⋯O distances of the shared proton with respect to the acceptor and donor heavy atoms, in conjunction with the donor–acceptor distance R_OO.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod47

126

Thus, distribution functions in terms of both the δ and R_OO coordinates, such as the one shown in Fig. 4, are utmost useful in characterizing H-bonds qualitatively (and of course quantitatively).

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod47

127

Symmetrized H-bonds, where the proton sits right in between donor and acceptor such as in idealized Zundel ions, are characterized by |δ| ≈ 0 Å.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod48

128

However, the coordinates δ and R_OO do not vary independently in real H-bonds, which can be gleaned from the anisotropic shape distributions such as those in Fig. 4.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res6

129

On the other hand strongly asymmetric H-bonds like those occurring in Eigen-type complexes have |δ| ≫ 0 Å; note that the sign of δ depends on the labeling convention of donor vs. acceptor.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod49

130

In conclusion, δ can also be considered a proton transfer coordinate that measures, as a function of time, as protons migrate along H-bonds; note that a more complex order parameter is needed in order to describe proton diffusion beyond a single H-bond.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod50

Evaluating the QM/MM interfacing: Water dimer and Eigen complex

131

As a first assessment of the quality of the interfacing between the QM and MM parts as established via H-bonds, the dipole moment of a QM water molecule as solvated by one MM (SPC) water molecule was compared to the proper all-QM treatment of the water dimer in the gas phase at 100 K; here the QM molecule is the donor water, i.e. schematically [HOH]_QM⋯[OH₂]_MM where ⋯ denotes the H-bond.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj10

132

In the all-QM case, the dipole moments of both water molecules can be obtained separately via the exact Wannier representation of the total dipole moment of the system, which can be decomposed approximately into linearly additive subsystem contributions stemming from the two molecular fragments, see ^{ref. 43 and 44} for background information.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod51

133

In the QM/MM case the dipole moment of the donor, i.e. the QM fragment, is obtained from the total dipole moment inside the QM box, whereas the one of the MM acceptor molecule is computed using the fixed SPC charges only.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod52

134

As seen from Fig. 5(a), the significant polarization of the donor water molecule upon H-bonding, leading to a significant increase of its dipole moment in comparison to the gas phase value of an isolated QM water monomer, is nicely captured by the interface.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs5

135

Actually, using SPC water as acceptor (instead of the “QM acceptor” as in the all-QM reference simulation) seems to even slightly overpolarize the QM subsystem.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs6

136

This would be consistent with the fact that SPC water is parameterized to represent bulk liquid water at ambient conditions, where the average dipole moment of the water molecules is significantly increased from a gas phase value, ≈1.85 D, to about 3 D due to H-bonding, see .^{ref. 44}

Type: Result | Advantage: None | Novelty: None | ConceptID: Res7

137

This consequence of parameterizing a non-polarizable water model to represent bulk water becomes evident by comparing the much larger average MM water monomer dipole moment in panel (b) to the QM monomer distribution in (a).

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj11

138

Furthermore, the MM acceptor in the QM/MM dimer acquires a much larger dipole moment as appropriate according to the reference data from the acceptor fragment in the all-QM treatment.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs7

139

This delicate point of coupling non-polarizable bulk water MM models via H-bonds to QM subsystems needs further study.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con2

140

We mention in passing that the coupling of the protonated water network to the BR backbone, see Section 3.4, is far more distant from the relevant H-bond(s) than in this test case where the “crucial” H-bond is cut by the interface.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs8

141

As a further check on the quality of the interface the microsolvated H₃O⁺ ion, i.e. H₃O⁺· (H₂O)₃, was treated in an all-QM approach and using a QM/MM approximation where the QM subsystem was the H₃O⁺ core only.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj12

142

The resulting two-dimensional H-bond distribution functions, see Fig. 6, look qualitatively similar, in particular the delicate correlation that shorter H-bonds are also more symmetric, i.e. the smaller R_OO the smaller |δ|, is nicely captured.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs9

143

The near to quantitative agreement is brought out more clearly in the corresponding one-dimensional representations along δ and R_OO as plotted in Fig. 7(a) and (b), respectively.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs10

144

Closer inspection shows that the coupling to SPC water leads to systematically shorter R_OO distances and thus to “stronger” H-bonds, consistent with the analysis of the dipole moment of the water dimer from above.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res8

145

Nevertheless, the “asymmetry” of the H-bonds as measured by δ is well represented in the QM/MM approximation to the all-QM simulation, see Fig. 7(a).

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs11

Evaluating the stability of the biomatrix

146

After having assessed the quality of the QM/MM interface as such, the next checks concern the stability of the biomatrix that hosts the QM subsystem, see Fig. 1.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj13

147

A critical aspect of this concern is the integrity of the interface between the lipid membrane and the solvating water layer which must be stable to provide a proper environment for the BR monomer.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj13

148

In order to assess this it is necessary to examine the structural stability of the system on time scales much larger then the “few picoseconds” typically accessible in ab initio simulations.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj13

149

To this end, we perform an analysis on an equilibrated 3 ns classical MM trajectory of the combined biomatrix consisting of the BR protein in a solvated POPC membrane; note that in this all-MM run the BR backbone was constrained as described in Section 2.3.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod53

150

The total mass density perpendicular to the membrane is decomposed in Fig. 8 into contributions coming from its major components.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs12

151

This analysis finds a stable layering in terms of lipid tails providing a hydrophobic membrane interior of approximately 20 Å thickness and an external layer of lipid headgroups which form a well-defined interface with the solvating water phase.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs13

152

From this graph one can see that the upper and lower portions of the (periodic) simulation cell are occupied solely by water in a layer which is approximately 6 Å thick with an additional 10 Å which is composed of over 95% water by weight.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs14

153

Most importantly, the density ≈0.97–0.98 kg/l that is observed in the region comprised entirely of water comes fairly close to the density of 0.988 kg l⁻¹ reported in ^{ref. 45} (Table 2) for bulk SPC water at 300 K and 1 bar.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res9

154

Finally, Fig. 8 also shows that the BR protein is well separated from its periodic image by over 14 Å, which corresponds to approximately four water solvation layers.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res10

155

Considering the lipid membrane in atomic level detail we first examine the lipid/water interface.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj14

156

There are on average about 77 water molecules within the first solvation shell (within a distance cutoff of 4.5 Å) of the 20 atoms representing the lipid head group with six to seven of these within the typical 2.7 Å H-bonding distance of the carbonyl and phosphate oxygens.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs15

157

This is in contrast to the 32 atoms representing the hydrophobic lipid chains where there is less than 12 water molecules within the first solvation shell, the majority of which are associated with the C atoms close to the head groups.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs16

158

Thus, the head groups of the lipids are well solvated whereas those of the side chains see very little of the solvent.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res11

159

The internal structure of the membrane can be analyzed in terms of orientational order parameters that describe the degree of order of the lipid tails along the chains.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj15

160

In particular, the order parameter SCCk = 〈3 cos (θ_k) − 1〉/2 is evaluated where θ_k is the angle of the kth carbon–carbon bond with respect to the normal of the lipid membrane surface, which we denote as the absolute z-axis.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod54

161

This is done for 14 carbon–carbon bonds starting from the head group down to the –CH₃ terminus of the hydrophobic lipid chain, see Fig. 9(a).

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod55

162

This analysis is done separately for the two chains in each POPC molecule, one of them, chain A, consisting of only C–C single bonds, whereas chain B contains one CC double bond at position k = 9.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod56

163

Both chains display their most pronounced order close to the head group, assuming a standard value of around S^CC ≈ 0.3; note that both layers contributing to the membrane were averaged together taking into account the underlying symmetry.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs17

164

As expected, the degree of ordering decays strongly towards the center of the lipid double layer for chain A and B. A major difference is that chain B does not decay smoothly, but shows a pronounced zig-zag anomaly centered around k = 9, i.e. the location of the CC unit in that chain.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res12

165

Inspection of the angular distribution functions p(cos (θ_k)) can be used to get a feel for the spread of the various bond angles along the chain.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod57

166

These distributions are shown for k = 2 and 14 in Fig. 9(b) for chain B. In general, for the bonds close to the lipid headgroup, i.e. small k, the distribution p(cos (θ_k)) is unimodal with a broad but well defined peak.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs18

167

This peak becomes broader as k increases, that is moving away from the headgroups into the hydrophobic lipid interior, and eventually smooths out with all angles having a sizeable probability in the limit of large k.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs19

168

The presence of the CC bond in chain B causes this broadening effect to occur at smaller k than in the overall more ordered chain A, which is in qualitative accord with description of the double bond perturbing the local order of this chain.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res13

169

Overall, it can be concluded that the POPC bilayer membrane (including the BR monomer) has a stable structure where the head groups are well solvated, the initial part of the lipid chains are orientationally ordered, and where the lipid tails created a hydrophobic interior that is well ordered.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con3

170

Most importantly, however, is that this structure shows no indication of degradation over the entire 3 ns of the all-MM run (after an additional 1 ns of initial equilibration).

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con4

171

Thus we have confidence that the biomatrix as set-up and relaxed provides a decent representation of BR embedded in a solvated lipid bilayer, which is stable on time scales 2–3 orders of magnitude larger than the QM/MM simulations.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con5

Embedding Eigen and Zundel complexes in BR

172

As outlined in the first section, the protonated water networks in BR are thought to be embedded via H-bonds to side groups of the protein.

Type: Hypothesis | Advantage: None | Novelty: None | ConceptID: Hyp2

173

In the present case, the three external water molecules of the Eigen complex H₃O⁺·(H₂O)₃ indeed connected via H-bonds to the groups Glu194, Glu204, Thr205, Tyr79 and Tyr57 of the BR backbone after switching from the all-MM run to the QM/MM simulation.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs20

174

The time evolution of these contacts in terms of the H-bond length R_OO is plotted in Fig. 10.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs21

175

These bond distances are found to fluctuate about meaningful average values without showing any tendency of systematic drifts or H-bond breaking, with the exception of the donor H-bonds of one of the terminal water molecule in the Eigen complex to an oxygen atom of the Glu194 residue.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res14

176

In the latter case several attempts to break this H-bond take place up to about 9.2 ps, where a successful breaking event occurs leading to a new arrangement for about 0.65 ps.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs22

177

Closer analysis of the H-bond angles of the respective H-bonds, see Fig. 11 shows that a water molecule is found to rotate about its local C_2v axis, thereby trading in one H-bond for another one.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res15

178

This is evident from the “flip” in both H-bond angles of a tagged water hydrogen atom to both Glu194 residue oxygen atoms (labeled O1 and O2 in Fig. 11).

Type: Result | Advantage: None | Novelty: None | ConceptID: Res16

179

Note that after about 0.65 ps a recrossing event occurs during which the system falls back into its previous H-bonding pattern.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs23

180

Another interesting fluctuation can be seen in Fig. 12 where the H-bond dynamics of another terminal water molecule in the Eigen complex is shown.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs24

181

This water molecule forms a donor H-bond with one of its hydrogen atoms to the Tyr57 residue, whereas the other hydrogen atom is bound to Tyr205 (see Fig. 10 for the O–O distances as a function of time in these H-bonds).

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs24

182

This three-fold coordinated water molecule (one strong acceptor bond to the core H₃O⁺ and two donor bonds to the protein residues) performs a rotation about the O–O axis of the acceptor H-bond to the core H₃O⁺ after about 9 ps.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs25

183

This results in a concerted “flip” in all four H-bond angles describing the two donor H-bonds to Tyr57 and Tyr205.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs26

184

Note that the concerted H-bond cleavage event is triggered by large fluctuations evident in the OH⁽²⁾–O^(Tyr205) H-bond angle, which also exhibits a recrossing event after about 0.4 ps.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res17

185

Eventually, this process results in the transient rupture of both H-bonds and subsequent reformation in a different H-bond topology.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs27

186

Overall, the observed H-bond dynamics of the protonated water network to the various protein residues follows a complex pattern with the H-bonds of the cluster itself being surprisingly stable considering the strong fluctuations occurring in the next solvation shell.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res18

187

Thus, the Eigen H₃O⁺·(H₂O)₃ QM cluster embedded in the MM environment can be considered to be dynamically stable on the picosecond time scale accessible via Car–Parrinello molecular dynamics.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con6

188

Equally important, the protonated network is still free to perform non-trivial fluctuations that involve the QM to MM interface such as H-bond exchange flips.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con7

189

Concerning the solvated Zundel complex, H₅O₂⁺· (H₂O)₄, there is a crucial difference in the set-up as compared to the Eigen case.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj16

190

The Zundel-like configuration was formed by taking the equilibrated Eigen-like complex, picking two water molecules from the hydration layer and moving them into positions in the QM pocket that allow for proper H-bonding to the core H₅O₂.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod58

191

Note that this procedure preserves the overall number of particle in our simulation.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod59

192

This new configuration has then been equilibrated via geometry optimization followed by running short QM/MM trajectories in the NVT ensemble.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod60

193

As opposed to the Eigen-like complex, which formed spontaneously during the QM/MM simulation, see Section 2.3, this approach did not result in a Zundel-like complex with |δ| being close to zero.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs28

194

In order to nevertheless analyze Zundel-like H-bonding situations in BR, the coordination number of the core Zundel H₅O₂⁺ was enforced to be four in the hope of stabilizing a symmetric H-bond with |δ| ≈ 0 Å.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod61

195

This was achieved by introducing four distance constraints connecting each of the four donor hydrogen atoms of H₅O₂⁺ to the acceptor oxygen atom of the respective solvent water molecule thus forming a solvated Zundel complex H₅O₂⁺·(H₂O)₄; here the equilibrium value of the corresponding gas phase cluster, 1.68 Å, was used.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod62

196

This constrained protonated water network stabilized itself by forming H-bonds to the residues Arg82, Tyr79, Glu194, Glu204, and Tyr205 inside the hydrophilic pocket.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs29

197

As will be shown by a detailed analysis in Section 3.4, the central H-bond in this solvated Zundel complex is close to symmetric.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs30

198

These findings suggest that the strong anisotropic and inhomogeneous nature of the protein matrix seems to lowers the probability of forming a Zundel-like complex, which is quite different to the situation found in the gas phase and in bulk liquid water.^46,47

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con8

199

Note that similar behavior was observed in studies of Eigen-like and Zundel-like complexes confined in hydrophobic channel environments.⁴⁸

Type: Background | Advantage: None | Novelty: None | ConceptID: Bac24

200

Despite these suggestive findings, we stress here that no conclusions concerning the absolute stability of Zundel-type versus Eigen-type protonated water networks can be drawn based on this qualitative finding.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con9

Structure and dynamics of Eigen and Zundel complexes: BR versus gas phase

201

After all these preliminary checks, the dynamical behavior of the Eigen and solvated Zundel complexes in BR as compared to the gas phase can be analyzed.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj17

202

In Fig. 13 the time evolution of their δ coordinate for only a very short time interval of half a picosecond out of more than 10 ps equilibrated QM/MM trajectories is displayed.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs31

203

Note again that the latter complex in BR is solvated by four additional water molecules saturating the H-bonds of the two terminal water molecules, i.e. H₅O₂⁺·(H₂O)₄.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs32

204

First, the general qualitative behavior seems quite similar in both the gas phase and in BR.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res19

205

In particular, the Eigen complexes are strongly asymmetric, δ ≪ 0, and thus display well-defined H₃O⁺ cores, whereas the Zundel-type complexes have δ fluctuating close to zero.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res20

206

However, it can already be gleaned that the H-bond in the solvated Zundel complex in BR is slightly asymmetric, although clearly not as pronounced as in the Eigen case, see below for more thorough evidence.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res21

207

The fluctuation pattern of the two complexes is fairly different since the Eigen complexes show high-frequency fluctuations clearly superimposed as a modulation on lower frequency oscillations.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res22

208

This pattern is much more congested and thus complex in the Zundel-type cases in the sense that there appears to be no clear-cut separation of time scales, which might also give a hint that the H-bond is more fluxional in this case since the proton can be easily shifted between “donor” and “acceptor”.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res23

209

It is revealing to look at the time evolution of the symmetry of the H-bond, i.e. |δ|, together with its instantaneous length R_OO as done in Fig. 14 for the Zundel complex in the gas phase.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj18

210

In this way one can nicely see the dynamical interplay of a lengthening the H-bond by increasing R_OO with a more asymmetric proton location in view of the concurrently increasing |δ| value.

Type: Object | Advantage: Yes | Novelty: New | ConceptID: Obj18

211

Thus, in particular the symmetric H-bond displays a complex dynamics where many modes are strongly coupled.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs33

212

A more quantitative analysis of the H-bonds is achieved by averaging over fluctuations and by considering two-dimensional probability distribution functions in terms of δ and R_OO, see Fig. 15 and 16.

Type: Model | Advantage: None | Novelty: None | ConceptID: Mod63

213

For the Eigen complex the one H-bond out of the three that has the smallest R_OO distance is used for each configuration; this is the “most active H-bond” as introduced in .^{ref. 46}

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs34

214

In case of the Eigen cation, it can be seen that on average a shortening of the H-bond clearly correlates with a more symmetric H-bond, i.e. smaller R_OO values are preferentially found in conjunction with smaller |δ| values.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res24

215

In addition, the distribution in BR is astonishingly similar to that in the gas phase.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res25

216

Concerning the slightly more localized character of the Eigen complex in BR it could be speculated that the strongly inhomogeneous electric field inside the protein favors localized charges—similar to what is known from the liquid phase.^46,47

Type: Result | Advantage: None | Novelty: None | ConceptID: Res26

217

The distribution function of the Zundel complex in the gas phase is nicely symmetric, see Fig. 16(a) whereas this is not the case in BR, see panel (b).

Type: Result | Advantage: None | Novelty: None | ConceptID: Res27

218

In both cases, for any given donor–acceptor distance R_OO the distribution function is much broader along the δ coordinate than the Eigen cation.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs35

219

This implies that the proton in the former complexes is clearly much more mobile than the latter, i.e. the proton transfer potential is expected to be fairly “flat”, which is still true for the less symmetric solvated Zundel ion in BR.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res28

220

The characteristics of the most active H-bond, in particular concerning the asymmetry in the case of the solvated Zundel complex in BR, comes out most clearly after averaging over the fluctuations of the donor–acceptor distance, see Fig. 17.

Type: Object | Advantage: None | Novelty: New | ConceptID: Obj19

221

For both Zundel-type complexes the peak is much broader than in the Eigen cases.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs36

222

This clearly shows that the proton is well localized on the central water molecule in the Eigen complex, whereas it is essentially equally shared (or “delocalized”) between the two neighboring water molecules in the Zundel complex.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs37

223

Since the relative free energy, i.e. the free energy profile or the potential of mean force, is related to the minus of the logarithm of such distribution functions, it is evident that the energetically most stable position of the proton in BR at room temperature is the hydronium core in the case of Eigen-like solvation, whereas the “delocalized” position is clearly favored in Zundel-like situations.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res29

224

Note, however, that the solvated Zundel complex could only be stabilized in BR by applying suitable (though mild) constraints on the four solvating water molecules, whereas the Eigen complex turned out to be very stable from the outset.

Type: Observation | Advantage: None | Novelty: None | ConceptID: Obs38

225

This is not too surprising in view of the inhomogeneous electric field due to the protein environment that most probably introduces asymmetry and thus destabilizes symmetrical solvation complexes.

Type: Result | Advantage: None | Novelty: None | ConceptID: Res30

226

Overall, there seems to be an astonishing structural similarity of the two protonated water networks in BR as compared to gas phase reference data of free protonated clusters.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con10

Conclusions and outlook

227

A detailed account of building and equilibrating a realistic QM/MM model of bacteriorhodopsin in a solvated lipid bilayer membrane was given where the protonated water network within a hydrophilic pocket on the extracellular side of BR was treated as the QM part.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con11

228

In addition, we have presented first results from QM/MM Car–Parrinello molecular dynamics simulations of two qualitatively different protonated water networks.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con12

229

We find that an Eigen-type cluster, i.e. H₃O⁺·(H₂O)₃ connected via H-bonds to suitable BR residues in this pocket, is a stable entity on the time scale of several picoseconds.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con13

230

In the same pocket, it was difficult to stabilize a solvated Zundel-type complex, i.e. H₅O₂⁺·(H₂O)₄, which was achieved only after introducing distance constraints on the hydrogen bonds that connect the solvating water molecules to the H₅O₂⁺ core.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con14

231

In this case, however, an essentially symmetric H-bond in the Zundel-core [H₂O⋯H⋯OH₂]⁺ was found, where on average the proton resides close to the bond midpoint between the two bridging oxygen atoms.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con15

232

The dynamics of this particular H-bond appears to be much more complex than that involving the Eigen-core H₃O⁺, in particular featuring fairly slow large-amplitude motion.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con16

233

Interestingly, striking similarities in structural and dynamic properties of these protonated networks embedded in BR and the analogues protonated water clusters in the gas phase exist.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con17

234

Most importantly these findings are consistent with the idea that protonated water networks might play a major role in the extracellular proton release process—possibly being at the heart of the proton release group “XH”.

Type: Conclusion | Advantage: None | Novelty: None | ConceptID: Con18