Modeling protonated water networks in bacteriorhodopsin

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.

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

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

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

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.

It is found that hydrogen-bonded networks around both H3O+ and H5O2+ cores can be stabilized in the protein matrix, leading to so-called Eigen and solvated Zundel complexes, H3O+· (H2O)3 and H5O2+·(H2O)4, respectively.

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

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.

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.

Motivation: Vectorial proton transport through membranes

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

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

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

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.

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.

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.

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

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.

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

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”.

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.

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.

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.

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.

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”.

Such ab initio, and in particular Car–Parrinello molecular dynamics approaches10 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.11

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

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.

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.

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.

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

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.

To this end, a recently developed QM/MM methodology12 based on a Car–Parrinello type approach10,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 software13–15 is used; see ref. 16 for an overview of this novel methodology including several ongoing applications.

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.

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

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

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.

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.

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

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.

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.20

Most interestingly, FTIR experiments21 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.

In addition, theoretical studies22 of pKa values suggest that the proton is stored in terms of a protonated water cluster, in particular H5O2+, rather than at the close by glutamic acid residues Glu194 and Glu204.

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.

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.

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.

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

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.

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).

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).

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

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.

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.

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.

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 H2O5+ model complex as suggested earlier22.

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 –NH2 group of Lysine.

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.

This necessitated the parameterization of all-atom Charmm26,27 as well as Gromos13 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.gref. 28. for a parameterization of the usually needed protonated Schiff base.

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 functional29,30 with the standard 6-31G** basis set using the Gaussian9831 quantum chemistry package.

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.

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

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.

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 fitting33 or Bader analysis.34

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.

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.

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

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

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

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.

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.

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.

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.

For the classical relaxation and equilibration runs the H5O2+ 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.

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

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.

The atomic charges on the Zundel ion (O: −0.82, Hterminal: 0.53, and Hbridge: 0.52) and the O–Hbridge stretching and O–Hbridge–O bending force constants (147 kcal mol−1 Å−2 and 32 kcal mol−1 rad−2, respectively), have been obtained in the same fashion as for the retinal/lysine complex Lyr from DFT calculations of isolated H5O2+.

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 H5O2+.

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.

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 package35,36.

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.

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

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

Phase B involves the use of the Gromos force field and software13–15 for preparing equilibrated starting configurations being used in the QM/MM production runs.12

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

(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.

(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.

(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.

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.

(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.

(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.

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

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.

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

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.

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.

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 reported23 water molecule Wat405, being also located in the hydrophilic pocket close to Glu204, was fixed.

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.

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

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.

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 H3O+·(H2O)3 with several H-bonds to close by residues from the pocket, see Section 3.3 for a detailed analysis.

QM/MM interfacing

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.12

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

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

The QM/MM interface12 combines this code to the Gromos software14,15 that represents the MM part in terms of the Gromos force field.

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

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).

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.

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.

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.

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

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.

For further details we refer to the original literature.12

In particular, the QM subsystem is treated here using Troullier–Martins pseudopotentials41 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.11

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.

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.

All dynamical DFT simulations were performed in the NVT ensemble employing Nosé–Hoover chain thermostats11 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.

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,11i.e. ∼10 ps, were generated within the QM/MM framework.

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.

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.

Simulation of protonated water networks in BR

Protonated water clusters H+·(H2O)n are quite stable molecular ions in the gas phase42 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”.

In the Zundel complex, H5O2+, the most probable position of the proton is midway between the two solvating water molecules, schematically denoted as [H2O⋯H⋯OH2]+.

This situation is very different for the Eigen complex, where the proton and one water molecule form a hydronium core H3O+, which itself is microsolvated by the three more molecules, i.e. H3O+·(H2O)3.

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

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 H5O2+·(H2O)4.

Standard structural parameters to investigate H-bonds in general, and in particular in protonated complexes, are the so-called asymmetric stretch coordinate δ = ROHRHO, 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 ROO.

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

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

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

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.

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.

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

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⋯[OH2]MM where ⋯ denotes the H-bond.

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.

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.

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.

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.

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

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).

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.

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

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.

As a further check on the quality of the interface the microsolvated H3O+ ion, i.e. H3O+· (H2O)3, was treated in an all-QM approach and using a QM/MM approximation where the QM subsystem was the H3O+ core only.

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 ROO the smaller |δ|, is nicely captured.

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

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

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).

Evaluating the stability of the biomatrix

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.

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.

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.

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.

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

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.

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.

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−1 reported in ref. 45 (Table 2) for bulk SPC water at 300 K and 1 bar.

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.

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

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.

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.

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

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.

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.

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

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.

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

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.

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.

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.

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.

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.

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.

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).

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.

Embedding Eigen and Zundel complexes in BR

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.

In the present case, the three external water molecules of the Eigen complex H3O+·(H2O)3 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.

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

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.

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.

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 C2v axis, thereby trading in one H-bond for another one.

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).

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

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.

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).

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

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

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

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

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.

Thus, the Eigen H3O+·(H2O)3 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.

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.

Concerning the solvated Zundel complex, H5O2+· (H2O)4, there is a crucial difference in the set-up as compared to the Eigen case.

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 H5O2.

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

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

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.

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

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

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

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.

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

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

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.

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

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.

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.

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. H5O2+·(H2O)4.

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

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

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.

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.

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”.

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

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

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

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 ROO, see Fig. 15 and 16.

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

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 ROO values are preferentially found in conjunction with smaller |δ| values.

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

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

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).

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

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.

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.

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

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.

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.

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.

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.

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.

Conclusions and outlook

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.

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

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

In the same pocket, it was difficult to stabilize a solvated Zundel-type complex, i.e. H5O2+·(H2O)4, which was achieved only after introducing distance constraints on the hydrogen bonds that connect the solvating water molecules to the H5O2+ core.

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

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

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.

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”.