INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Proton transport across cell membranes is an essential process in cell physiology. This process is facilitated and regulated by large transmembrane (TM) proteins, such as ATP synthases or bacteriorhodopsin. The structural complexity of these proteins, however, makes it extremely difficult to dissect the molecular mechanisms of protein-mediated proton transport through the nonpolar membrane interior. It is thus desirable to have a protein model that is small, has a well-known structural motif, yet operates with the efficiency and control of more complex proteins. This has led to the study of the influenza A M₂ proteina small, homotetrameric, voltage-gated ion channel that transports protons with high efficiency and selectivity (Pinto et al., 1992; Wang et al., 1993; Sakaguchi et al., 1997). Each monomer is built of 97 amino acids and contains a single TM domain 19 residues long. Not all residues, however, are essential for transport. Active channels have been reconstituted from a synthetic peptide containing only a subset of 25 residues, including the TM region, with no loss in specificity or efficiency (Duff and Ashley, 1992). The sequence of amino acids in the peptide is given in Fig. 1. The remarkable combination of simplicity and efficiency makes M₂ an attractive subject for studies on the mechanism of charge transfer across membranes and a potential target for reengineering to create a simple proton pump.

View larger version (37K): [in this window] [in a new window]   FIGURE 1   Cartoon of the M2 bundle with indicators to locations of key residues. The protein sequence is indicated below.

Although no high-resolution structural data for M₂ are available to date, recent NMR and CD studies of this protein in phospholipid bilayers strongly suggest that the TM region is -helical and has the quaternary structure of a four-helix bundle (Duff et al., 1992; Kovacs and Cross, 1997). The center of the bundle forms the transmembrane pore. This model, shown schematically in Fig. 1, is supported by a recent simulation of M₂ in a water-octane lamella (Zhong et al., 1998). If the external environment is brought to a pH of ~5.3, and a suitable pH gradient exists across the membrane, protons are transferred through the channel from the external environment, thus acidifying the interior of the cell. Although initial measurements showed that M₂ is also permeable to other small cations, such as Na⁺ (Shimbo et al., 1996), more recent experiments have pointed to possible errors in these measurements (Chizmakov et al., 1996) and have yielded very low efficiencies of transfer of these cations (Dieckmann and DeGrado, 1997; Chizhmakov et al., 1996). In fact, the permeability of M₂ to Na⁺ relative to H⁺ was estimated at 6 × 10⁷ (Chizhmakov et al., 1996). This, in turn, points to the presence of a gating mechanism within the channel, which allows protons to move through the pore but blocks all larger ions. The titratable group implicated in the gating is the intraluminal His³⁷ residue (Pinto et al., 1997; Wang et al., 1995). In the most widely accepted structural models of M₂, a ring of four histidine residues, one from each of the helices, occludes the pore (Sansom et al., 1997; Pinto et al., 1997; Zhong et al., 1998), thus providing an ideal control point inside the channel. The histidines are oriented approximately perpendicular to the water-membrane interface with the unprotonated -nitrogen atoms of the imidazole rings located on the extracellular side of the gate and the protonated -nitrogen atoms located on the intracellular side.

The main aim of this study is to gain insight into the molecular mechanism of proton transport and gating in M₂. To date, two mechanisms have been proposed (Sansom et al., 1997; Pinto et al., 1997). For both it has been assumed that protons are translocated along the network of water molecules that fills the pore of the channel. This is analogous to the mechanism of proton transfer through another channel, gramicidin A (Akeson and Deamer, 1990). The two mechanisms, however, differ substantially in the proposed way in which protons navigate the gate formed by the histidine residues. In one mechanism (further referred to as mechanism I), protons are transferred via a "proton shuttle" involving a histidine residue. This is shown schematically in Fig. 2 a. In the initial state, an -protonated histidine can accept a proton from the external environment on the unprotonated -nitrogen atom of the imidazole ring. This results in a biprotonated intermediate that is positively charged. This intermediate can relax by eliminating either the newly acquired -proton back into the environment or the -proton into the opposite side of the channel, thus shuttling the proton across the membrane. Then the system returns to the initial state through tautomerization of the imidazole ring.

View larger version (24K): [in this window] [in a new window]   FIGURE 2   Sketch of two proposed mechanisms of channel gating by M2. (a) A "proton shuttle" process (mechanism I in this text), whereby proton transfer is mediated by the intraluminal His37 residues and is followed by regeneration of the starting configuration through tautomerization of the imidazole ring. (b) A purely structural model (mechanism II) in which the open and closed states of the gate are mediated by steric forces resulting from protonation of one or more His37 residues.

In the alternative mechanism, all four histidines acquire an additional proton and, because of repulsion between their positive charges, move away from one another, thus opening the channel. In the open state, a chain of water molecules penetrates the gate to span the full length of the pore. Protons can be efficiently transferred along this chain as long as the gate remains open. This mechanism, shown schematically in Fig. 2 b, will be called mechanism II. It has been proposed on the basis of the results of constrained molecular dynamics simulations of a model of M₂, solvated at both ends by water molecules and placed in a dielectric medium parameterized to represent a membrane core (Sansom et al., 1997).

In principle, the correct mechanism can be identified through a series of accurate simulation studies of the energetics and dynamics of proton transfer through the M₂ channel embedded in a hydrated phospholipid bilayer. In practice, however, such an approach might not be feasible and certainly would be extremely demanding. Because both proposed mechanisms involve breaking and forming chemical bonds, a methodology that combines classical statistical mechanics and quantum mechanics would be desirable. Even if the system is treated in a classical approximation (Sagnella and Tuckerman, 1998), many complications arise, for example, the need to explicitly include polarization effects (Pomes and Roux, 1998). It is therefore appropriate to initially take a simpler approach, exploiting the fact that each of the proposed mechanisms involves specific but different predictions about the intermediate states of the protein during proton transfer. In this approach, we require that all intermediate states must be stable and that the orientation of water molecules inside the channel must be conducive to proton transfer. If a candidate mechanism does not fulfill these requirements, it is unlikely that it provides a correct explanation of channel action.

Each of the intermediate states for mechanisms I and II has been investigated by molecular dynamics computer simulations. For mechanism I, we studied the system corresponding to the initial state consisting of the four -protonated His³⁷ residues (system 1a). This state has already been simulated in a membrane-mimetic octane lamella (Zhong et al., 1998) and, recently, in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) membrane (Forrest and Sansom, 1998). Next, we considered a mixed state, consisting of three  and one biprotonated histidine (system 1b). This is an intermediate in which the proton being transferred is located at the gate (see also Fig. 2 a). Finally, we considered a three , one  configuration (system 1c). This state results from a successful proton transfer but precedes tautomerization of the histidine residue, which restores the initial state. To test mechanism II, we simulated an intermediate containing four biprotonated His³⁷ residues (system 2a). The same system was studied previously in octane lamella (Zhong et al., 1998). In addition, we also simulated a state with two  and two biprotonated species (system 2b) to investigate the possibility that partial protonation of the histidines may be sufficient to open the gate and establish a stable water wire spanning the channel. The simulations performed are summarized in Table 1.

View this table: [in this window] [in a new window]   TABLE 1   Summary of system configurationsprotonation states of His37

After describing the method of calculations, we present the results for system 1, which test mechanism I, followed by the results for system 2, which corresponds to mechanism II. Then we discuss implications of our results for the mechanisms of proton transport in M₂.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The protein was represented as an all-atom model terminally blocked by Ac- and -NHMe groups on the N- and C-terminal, respectively. The initial structure was a left-handed, coiled-coil, -helical bundle composed of the 19 TM residues (residues Leu²⁶ through Asp⁴⁴) of M₂. This structure was obtained by computer model building (Dieckmann and DeGrado, personal communication). The modeled protein fragment was subsequently extended to 25 residues per monomer (Ser²² through Leu⁴⁶), which coincides with the sequence used by Duff and Ashley and Kovacs and Cross in their reconstitution and structural studies (Duff et al., 1992; Kovacs and Cross, 1997). This modification was required because preliminary simulations of the channel in a water-hexane membrane mimetic showed that the original, shorter model was unstable (the monomers diffused away from each other). There were also functional reasons to extend the sequence. It contains three sets of charged residuesAsp²⁴, Asp⁴⁴, and Arg⁴⁵which form "charged rings" at the ends of the bundle. The rings are schematically shown in Fig. 1. All residues forming these rings were simulated in the ionized form. Ionization states of these residues, however, have not yet been determined experimentally and may be different from those assumed here. Recent theoretical calculations of the ionizable side chains in the channel-forming helices of the nicotinic acetylcholine receptor have led to the suggestion that the residues comprising the charged rings on the extracellular side of the channel are protonated at neutral pH, while the intracellular and intermediate rings are fully ionized (Adcock et al., 1998). It is possible that this is also the case for influenza M₂ protein. It should be kept in mind, however, that the viral channel lacks a significant extramembrane domain, which may influence ionization states of the nearby residues in the acetylcholine receptor.

The protein was embedded in a bilayer built of 90 dimyristoylphosphatidylcholine (DMPC) molecules and placed between two lamellae of water, each containing 2510 molecules. An appropriate number of sodium ions was added to maintain system neutrality. DMPC is a particularly suitable membrane model because the same phospholipids were used in experimental studies on M₂ (Kovacs and Cross, 1997). The dimensions of the simulation box were 65.0 × 65.0 Å² in the x,y directions, parallel to the interface, and 200.0 Å in the z direction, perpendicular to the interface. This yields an "open" geometry of the system in which water vapor exists in equilibrium with the water lamellae. To eliminate undesired edge effects, periodic boundary conditions were applied in the three spatial directions.

Water molecules were represented by the TIP3P model (Jorgensen et al., 1983); the DMPC membrane model was developed and tested by Berkowitz et al. (Essmann et al., 1995a). To describe the protein, the AMBER force field of Cornell et al. (1995) was used with a Lennard-Jones cutoff of 9.0 Å. The Lennard-Jones parameters between different components of the system were obtained from the standard Lorentz-Bertholot combination rules.

For each system listed in Table 1, molecular dynamics simulations were performed at the physiological temperature of 310 K in the NVT ensemble for ~1.3 ns after equilibration. This temperature is 13 K above the phase transition temperature of DMPC (Shen et al., 1997), which ensures that the bilayer is in the liquid crystal state. Equations of motion were integrated with a time step of 2 fs. The AMBER 4.1 simulation package (Pearlman et al., 1995), with long-range electrostatic interactions included via the particle mesh Ewald (PME) algorithm (Essmann et al., 1995b), was used. Including the long-range correction was necessary because the phospholipid headgroups carry large atomic charges that interact significantly over distances extending beyond the primary simulation cell. All simulations were preformed on a four-processor SGI Origin 200.

Before the production simulations were started the system was equilibrated. This was done in several steps. First, close interatomic contacts resulting in large van der Waals repulsive energies, found in certain parts of the initial model of M₂, were relieved by performing several short simulations (~20 ps each) of the protein in vacuo. To ensure that the overall geometry of the bundle did not significantly change during relaxation, only atoms with close contacts were allowed to move, and the atomic velocities were frequently rescaled.

Once close contacts between atoms were removed, the protein was inserted into an equilibrated configuration of a water-DMPC membrane system. Water and membrane molecules that overlapped with the protein were removed and Na⁺ ions were added to keep the system neutral. Starting from this configuration, the system was simulated for 500 ps at constant pressure with independent scaling of the x, y, and z dimensions of the box. This allowed the water-membrane system to relax around the bundle, which was held in place to keep the monomers from diffusing away as the membrane equilibrated around them. Once the initial equilibration was achieved, the constraints on the bundle were released and the system was simulated for an additional 800 ps. Simulations of the same length were performed to equilibrate other systems (in different protonation states of His³⁷ residues).

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Stability of system 1

The stability of the whole system and, specifically, the protein can be characterized by following the time evolution of several aggregate structural indices. The packing arrangement of the bundle can be conveniently described by four types of parameters: the average number of residues per turn n_i, the average rise per residue d_i (i = 1, 2, 3, 4) of the individual helices, the interhelical separation R_ij^min, and the crossing angle _ij between those members of the bundle that are in contact (Chothia et al., 1981). In the initial model, the values of d_i and n_i for each helix were 1.51 Å and 3.86, respectively. By comparison, the corresponding values for an ideal -helical domain are 1.5 Å and 3.6 (Creighton, 1993). In real systems, d_i varies from 1.4 to 2.0 Å, and n_i varies from 3.2 to 3.9 (Chothia et al., 1981). The pairs of helices that are in contact in the model of M₂ are (0, 1), (0, 2), (1, 3), and (2, 3). In this model, all of them are separated by an interhelical distance R_ij^min of 9.5 Å. The crossing angle for all four pairs of helices is 17°. In other proteins these parameters can vary considerably. The values of parameters for the initial model are given here as a reference rather than as targets. Because the model is only an approximate structural solution, consistent with certain experimental findings rather than the precise description of the M₂ structure, deviations from the initial values are expected.

The time evolution of the four structural properties is shown in Fig. 3. The three systems exhibit very similar structural properties; therefore we present the data for system 1a only. The average values for n_i, d_i, R_ij^min, and _ij are ~3.9, ~1.5 Å, ~9.5 Å, and ~38°, respectively. The only quantity that markedly deviates from that in the model is the crossing angle. The flux of water in the channel, the degree of solvation of the protein ends, and the interaction with the membrane are all likely to influence this parameter. It is unclear, however, which is the primary factor affecting these changes. The main point emerging from the data shown in Fig. 3 is that the structural features of the bundle are stable over the 1-ns simulation. No systematic drift that would indicate disruption of the channel is observed. Substantial fluctuations occur only in the values for n₁ and n₃ for system 1a. Visual inspection of the protein in this system reveals that helices 1 and 3 are distorted slightly, which yields an overall bend along their long axes. This causes difficulties in determining structural parameters because they depend on these axes, which, in bent structures, can be defined only approximately.

View larger version (43K): [in this window] [in a new window]   FIGURE 3   Helix parameters as a function of simulation time for system 1a. For the top two plots, the residues per turn and the rise per residue are plotted for each helix in the bundle. For the lower two plots, the separation between helices and the crossing angles are shown for monomers that are in contact. Individual monomer helices are labeled 0, 1, 2, and 3.

The conclusion about stability of the three protein systems is reinforced by the data on the helicity of the monomers as a function of time (not shown). Calculations of this quantity with the DSSP program (Kabsch and Sander, 1995) reveal that each of the monomers remain -helical throughout the simulations. The ends of the monomers fray slightly, because of solvation by water, but the overall degree of helicity for each monomer in each of the three simulations is greater than 90%.

Another measure of the stability of the protein is the calculated root mean square deviation (RMSD) as a function of time. This is shown in Fig. 4. We chose to use the last configuration of each trajectory as the reference structure. Proceeding backward in time from the last time point, there is an initial adjustment of the structure for all systems, followed by a stable RMSD. No strong drift in RMSD, which would signal a potential structural instability, was observed.

View larger version (16K): [in this window] [in a new window]   FIGURE 4   Root mean square deviation of the enitire bundle as a function of simulation time for systems 1a, b, and c (legend shown as inset). The RMSD is calculated with respect to the last configuration of a given simulation; thus the curves converge to 0 as the ends of the simulations are approached.

Next we examine the orientation of the helical bundle in the lipid bilayer. In all of the simulations the bundle remains tilted, on average, by ~30° with respect to the bilayer normal. The tilt of individual helices as a function of time in system 1a is shown in Fig. 5. The results for other systems are quite similar. Each helix in the bundle is tilted to a different degree (helix 1 is tilted ~10°, while helix 3 is tilted 45°), which is expected if the bundle as a whole is tilted. In addition, the orientation of each helix in the bilayer is stable, i.e., it is not changing markedly during the simulations.

View larger version (20K): [in this window] [in a new window]   FIGURE 5   Time evolution of the tilt angle of each helix in system 1. The cosine of the angle is plotted, where the angle between helix axes is with respect to the interface normal (z axis).

Tilting of individual -helices in lipid bilayers has previously been observed in molecular dynamics simulations of polyalanine embedded in a membrane (Shen et al., 1997). This was rationalized by pointing out that if the hydrophobic core of the membrane is shorter than the length of the TM region of the protein, the helix will tilt to maximize its contact with the core. Furthermore, our results are consistent with NMR data indicating that the individual helices of M₂ exhibit a 30° average tilt relative to the bilayer normal (Kovacs and Cross, 1997). It should be pointed out, however, that the tilt of individual helices does not necessarily imply that the whole bundle must be tilted. Because the M₂ geometry is somewhat funnel-shaped, the helices would be tilted even if the pore of the channel was aligned with the bilayer normal. In recent simulations of M₂ in a membrane mimetic, the bundle did not appear to exhibit significant tilt (Zhong et al., 1998). However, the time evolution of this quantity was not presented, but, instead, only a few snapshots of the system were shown, making it difficult to ascertain the arrangements of the protein during the simulation.

The discussion of stability, limited so far to the protein, can be further extended to include other components of the systems. Good global indices of system stability as a whole are density distributions (Z) of each component measured along the bilayer normal (the z axis) and centered in the middle of the bilayer. To test for a systematic drift in these profiles, the trajectories were broken into segments of several hundred picoseconds, and the density profiles were compared. No drift in these quantities was observed (data not shown), thus indicating the stability of the water-membrane-protein assembly. The distributions for system 1a, averaged over the full length of the trajectory, are shown in Fig. 6. The distributions for systems 1b and 1c are quite similar and, therefore, are not shown here. In each case, the membrane component exhibits well-defined headgroup and hydrocarbon tail regions. The observed structure of the water-membrane interface is consistent with results of other simulations of membranes in the absence of integral membrane proteins (Pohorille and Wilson, 1995; Essmann et al., 1995a; Perera et al., 1996; Marrink et al., 1996). The only difference is the apparent larger width of the interface in the presence of the protein. This is not a result of the deeper penetration of water into the membrane headgroup region, but reflects the presence of water inside the protein channel. Similarly, the asymmetry of the water density profiles on both sides of the membrane is associated with unequal distributions of water molecules along the pore of the protein.

View larger version (30K): [in this window] [in a new window]   FIGURE 6   Z dependent density profiles for components of system 1, with the center of the membrane at Z = 0. (Top) Water and membrane distributions and total protein distribution. (Bottom) Distributions of selected protein residues. Legends are shown as insets.

In each system, the protein is symmetrically distributed in the bilayer, with its extracellular N-terminus on the right (z > 0) and intracellular C-terminus on the left (z < 0). The overall length of the protein along z is ~45 Å, and the ends of the protein are well solvated by water. By comparing distributions in the two panels of Fig. 6, it can be found that charged residues of the protein (Asp²⁴, Asp⁴, Arg⁴⁵) are colocated with the zwitterionic membrane headgroups, which form a highly polar environment. The resulting strong electrostatic interactions between these residues and the membrane may help anchor the protein in the bilayer and limit the diffusion of the monomers to the plane of the membrane. In contrast to the charged residues, the gate-forming histidines are located in the hydrocarbon core of the membrane, slightly shifted from the center of the bilayer toward the intracellular side. The distributions of the charged residues and the histidines do not change with time, once again indicating that the system is stable during the simulations.

Evaluation of gating predicted by mechanism I

As already mentioned, it is required not only that all intermediate states must be stable, but also that the orientation of water molecules inside the channel must be conducive to proton transfer. Recall that in mechanism I, a proton reaching the gate becomes attached to a His³⁷ residue on the extracellular side in the  position, and, subsequently, the  proton is released on the intracellular side of the gate. The cycle is completed through tautomerization of the histidine residue. Thus this mechanism implies that the gate remains closed to water at all stages of the cycle. Whether this is the case for states 1a, 1b, and 1c can be examined by calculating the water density profiles in the channel. Because the channel is tilted in the bilayer, the channel axis  is a better coordinate to represent these profiles than the bilayer normal. The positions of water molecules are measured relative to the centroid of the His³⁷ residues. As can be seen from the distributions shown in Fig. 7, the channel is nearly completely occluded near the gate in all three states. In general, the average number of water molecules is small in most of the transmembrane portion of the channel but increases markedly toward the ends, as the channel widens.

View larger version (20K): [in this window] [in a new window]   FIGURE 7   Number of water molecules as a function of distance along the bundle axis (), with the histidine centroid as the center. Plotted separately are data for systems 1a (top), 1b (middle), and 1c (bottom).

Although the number of water molecules in the gate region is quite small, it is not precisely zero. Thus it may be suggested that occasionally a water wire forms across the gate that is sufficiently long-lived to transfer protons across the occlusion zone, providing a transport mechanism alternative to the "proton shuttle" mediated by a histidine residue. This hypothesis can be tested by calculating the persistence time of water clusters spanning the gate.

The calculations required several steps. For each time step, the water molecules, the oxygen atoms of which were within 15 Å of the channel axis, were selected. Among them were all water molecules inside the channel but also some that were located near its edges. Next, for each selected water molecule, its neighbors within an oxygen-oxygen cutoff distance of 3.5 Å were identified. This allowed for the construction of an adjacency table whose ith, jth entry was 1 if water molecules i and j were neighbors or 0 otherwise. Such a table is a representation of a directed graph (Kruse et al., 1997), and thus water molecules could be assigned to clusters by a process known as graph traversal. We implemented a depth-first algorithm to traverse the graph (Kruse et al., 1997). If all water molecules within a 6.0-Å cutoff from the His³⁷ centroid were found to belong to one cluster, then the gate was considered to be open. If the water molecules were split between two or more clusters, then the gate was considered to be closed. Once a water wire across the gate was identified, its persistence time was calculated from the number of consecutive configurations (time steps) in which it was observed.

The number of events where a water wire formed across the gate and persisted for at least time t ps is shown in Fig. 8. For all systems, these values decay to almost zero in less than 4 ps, indicating that there were no water wires that lasted longer than this. In fact, the longest-lived water wire observed in the simulations lasted for slightly less than 7 ps. Both the scarcity and the short lifetime of water wires across the gate suggest that such wires do not provide an efficient mechanism for proton transport if only one histidine undergoes protonation.

View larger version (24K): [in this window] [in a new window]   FIGURE 8   Probabilities of finding water wires across the gate that persist for at least t ps. The details of how this was calculated are described in the text. For all systems, these probabilities decay to almost zero in less than 4 ps. The longest-lived water wire observed in the simulations lasted for less than 7 ps.

Water molecules near the gate can be further characterized by calculating the radial distribution functions, P(r_ij), between the nitrogen atoms in the His³⁷ imidazole ring and water molecules in the channel within a 10-Å cutoff. Four different distributions can be defined for each histidine residuetwo N-O^wat and two N-H^wat correlations (using either the  or  nitrogen atom of a His³⁷ as the origin). Note that system 1a contains four His³⁷ residue of one type only. They are in an  protonation state and will be further abbreviated as HIE. In system 1b, there are three such residues, and one double-protonated histidine (HIP) carrying a net positive charge. In system 1c, there are three HIE and one  protonated histidine (HID) residue. The protonation states for each system are listed in Table 1.

All four radial distribution functions for the three types of histidine residues are shown in Fig. 9. All P(r_NO) exhibit a peak at a typical N-O hydrogen bonding distance of 2.86 Å, which indicates that water molecules in contact with the His³⁷ residues form a well-defined solvation shell around the imidazole nitrogen atoms. In system 1a, the peak in P(r_NO) is small, indicating that the average number of water molecules residing near each -nitrogen atom of HIE is relatively small. In fact, the corresponding hydration number, obtained by integrating P(r_NO) from zero to the first minimum located at 3.5 Å, is equal to 1.0. The first peak in P(r_NO) is even smaller and the corresponding hydration number is equal to 0.55. For both types of nitrogen atoms, there is also a small peak at 4.8 Å. The 3-Å separation between the first and second peaks is close to a typical oxygen-oxygen distance between two hydrogen-bonded molecules of water.

View larger version (36K): [in this window] [in a new window]   FIGURE 9   Radial probability distributions between imidazole nitrogen atoms and surrounding water molecules for His37. The origin is placed on the nitrogen atoms. Distributions for the  N are on the left, and those for the  N are on the right.

The suggestion that the  nitrogen of HIE forms a hydrogen bond with the neighboring water molecules is confirmed by examining P(r_NH). This distribution exhibits peaks at 1.96 and 3.26 Å from the nitrogen atom, consistent with an arrangement in which one hydrogen atom of the water molecule is approximately aligned with the N-O axis and the other forms an angle of ~105°. In contrast, P(r_NH) exhibits only one peak at 3.5 Å, which indicates that both hydrogen atoms of the water molecule point away from the N-H group of HIE. The distribution functions for HIE in systems 1b and 1c are similar and, therefore, are not shown here. The only notable difference is that the first peak in P(r_NO) in system 1b is somewhat smaller than in the other two systems.

The radial distribution functions around HID in system 1c are also similar. The main difference is that the patterns of hydration of  and  nitrogen atoms are reversed, compared to HIE, which is consistent with the fact that, in HID, the  nitrogen atom is protonated while the  nitrogen atom is not. Furthermore, the second peaks in the distribution functions are somewhat larger, indicating a better long-range ordering of water molecules.

The distribution functions for HIP in system 1b clearly differ from those for the other two types of histidine residues. The distributions around  and  nitrogen atoms are nearly identical and suggest the existence of N-H···O hydrogen bonds, with both hydrogen atoms of the water molecule pointing away from the histidine. This agrees well with the fact that both nitrogen atoms in HIP are protonated. Furthermore, the first peaks in the distribution functions are much larger than the peaks in distributions around nitrogen atoms in the other two types of histidine residue (HIE and HID). The corresponding hydration numbers are nearly identical (2.0 and 1.9, respectively), because both  and  positions are approximately equivalent. In addition, the second peaks are very well defined. This shows that the presence of a charge on the histidine exerts an attracting influence on water molecules and induces an order in water molecules in the channel that extends for at least two hydration shells.

To confirm our deductions regarding water orientation around the His³⁷ residues we calculated distributions of the angle N-O^wat-H^wat for each type of imidazole nitrogen. These distributions are shown in Fig. 10. For all unprotonated nitrogen atoms, the orientational distributions exhibit two clear peaks located at angles of 0° and 104°. This means that one hydrogen atom in the water molecule points directly at the unprotonated site along the N-O axis. It further indicates that such an orientation is fairly rigid. The orientational distributions of water molecules near the protonated nitrogen atoms are much broader, suggesting that these water molecules undergo a wide range of motions. The common characteristic of these motions is that hydrogen atoms preferentially point away from the histidine. It appears that the most common arrangement of water molecules around HIP is such that its C_2v vector lies along the N-O axis and points away from the nitrogen atom (the N-O^wat-H^wat angle is approximately equal to 130°). For uncharged HIE and HID residues, there is also a significant fraction of water molecules oriented such that one hydrogen atom points directly away from the nitrogen atom, forming the N-O^wat-H^wat angle of 180°. All of these results are consistent with the previously discussed radial distribution functions. The orientations of water molecules near different types of histidine residues are schematically summarized in Fig. 11.

View larger version (27K): [in this window] [in a new window]   FIGURE 10   Orientational distributions for water in the first solvation shell of the imidazole nitrogen atoms (as defined by the first minima in the radial distribution functions shown in Fig. 9). The solid line corresponds to the angle defined by N-Owat-Hwat. The light dashed line shows the angle defined by using the  N in place of the  N.

View larger version (21K): [in this window] [in a new window]   FIGURE 11   Sketch showing the deduced orientation of water about the His37 imidazole ring, for the different protonation states.

Proton transfer capabilities depend not only on the correct orientation of water molecules in the channel but also on the lifetimes of these orientations. One way to gauge the orientational dynamics of a water molecule is by calculating the time correlation functions for its dipole and H-H vector, defined as

(1)

where P_l refers to the lth-order Legendre polynomial, and is a vector in the molecular frame (either the water dipole or H-H vector). Here C₁ for the dipole vector and C₂ for the H-H vector were calculated by averaging over all time origins  and over all molecules in a specified region. Three types of water molecules were consideredbulk water, water at the water-membrane interface, and solvation shell water (within 3.5 Å of any His³⁷ imidadzole nitrogen atom). The definition of the His³⁷ solvation shell is based on the first minimum in the N-O radial distance distributions shown in Fig. 9.

The results for system 1a are shown in Fig. 12. The correlations for other systems are very similar and thus are not shown. This similarity indicates that the protonation state of the His³⁷ residue does not markedly affect the water dynamics. The results further show that there is a substantial difference in the water dynamics between the three regions. In particular, water molecules in the solvation shell of the His³⁷ residues reorient much more slowly than those in the bulk region and, to a somewhat lesser extent, at the interface. This slower motion of water near the His³⁷ residues could aid proton transfer by keeping the geometry of the donor and acceptor molecules fixed.

View larger version (22K): [in this window] [in a new window]   FIGURE 12   Time correlation functions for reorientation of the water dipole and H-H bond for bulk water, water at the membrane surfaces, and in the solvation shell of the His37 residues.

Evaluation of gating in mechanism II

In this set of simulations we test the effects of multiple protonations of the His³⁷ residues. System 2a consists of the fully protonated state, whereas system 2b has only two of the His³⁷ residues protonated. It is expected that the confinement of the charged His³⁷ residues to a small region in the lumen of the channel will be electrostatically unfavorable in the absence of some neutralizing entities, such as counterions. This may cause the histidine residues forming the gate to "swing open," as shown schematically in Fig. 2 b, leaving sufficient space for water molecules to penetrate the gate.

For a fully protonated state, containing four HIP residues (system 2a), two molecular dynamics trajectories, started from different initial conditions, were obtained. In one trajectory, the gate was initially in the closed state. In the other trajectory, the protonated His³⁷ residues were first moved to the "open" position characterized in the simulations of Sansom et al. (1997). In both cases the geometry of the helical bundle was initially constrained to allow for adjustments of the gate structure without destroying the integrity of the channel. Once the constraints were released, the fate of the system was the same in the two simulationsthe channel disintegrated within 280-450 ps, with one monomer separating from the remaining three. The separation was preceded by increasing deformation of the bundle and is reflected in the RMSD, as shown in Fig. 13. A representative structure of the disrupted bundle is shown in the top panel of Fig. 14, whereas a typical configuration for the intact channel (with waters) is shown in the bottom panel.

View larger version (14K): [in this window] [in a new window]   FIGURE 13   Root mean square deviation of the enitire bundle as a function of simulation time for systems 2a and b. The rest of the caption is the same as for Fig. 4.

View larger version (31K): [in this window] [in a new window]   FIGURE 14   Snapshot of the M-2 bundle for the fully protonated system (system 2a) as it disintegrates (top) and an intact bundle with channel water (bottom). The carbon backbones are represented as ribbons, and the His37 residues are shown in CPK.

The disruptive effects of electrostatic repulsion should be markedly reduced if the gate is only partially protonated. This led us to examine a system in which only two His³⁷ residues are biprotonated (system 2b). The helix parameters for this system, shown in Fig. 15, are fairly stable, with the exception of n₂ (shown in the top panel). This helix appears to have a bend near the histidine residue, but the regions on both sides of the bend remain -helical. Because the helix is not straight, the helical parameters are less accurate. During the latter part of the simulation, the helix becomes less kinked, but simultaneously the histidine residues move apart slightly, creating a small gap in the gate region. Overall, the bundle flattens, as though the cylinder had been compressed, and loses its conelike geometry. These changes are reflected in the RMSD, shown in Fig. 13. Despite rather significant shape changes, the monomers remain associated and provide an adequate channel through the membrane. The positioning of the bundle in the membrane is not strongly affected by the changes in geometry.

View larger version (43K): [in this window] [in a new window]   FIGURE 15   Helix parameters as a function of simulation time for system 2b. The rest of the caption is the same as for Fig. 3.

It is expected that multiple protonation of His³⁷ residues will increase water penetration through the gate region. As can be seen in Fig. 16 for system 2b, the gate is only slightly more open relative to that seen in system 1. Furthermore, the lifetimes of water wires spanning the gate remain quite short (see Fig. 8). Both radial and orientational distributions of water molecules around HIP and HIE residues are similar to the distributions around the same residues in system 1b. In general, system 2b undergoes large changes in bundle shape, with little accompanying change in features important to the conduction of protons across the channel. Thus it appears that the channel is a robust structure capable of accommodating significant changes in bundle geometry. Simultaneous protonation of two His³⁷ residues does not appear to result in excessive amounts of water entering the gate. However, because the simulations do not extend beyond 1 ns, it is not possible to ascertain whether the system remains intact over longer time periods.

View larger version (10K): [in this window] [in a new window]   FIGURE 16   Number of water molecules as a function of distance along the bundle axis () for system 2b.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

One of the main objectives of this paper was to examine the structural stability of different protonation states of the M₂ channel involved in mechanisms I and II. All three protonation states required for mechanism I were found to be stable, while the channel containing four biprotonated His³⁷ residues, crucial to mechanism II, disintegrated during the simulations. These results clearly lend support to the "proton shuttle" mechanism (mechanism I). This conclusion, however, should be drawn with some caution. Each molecular dynamics trajectory extended for only ~1 ns (excluding equilibration). This is a short time compared to time scales of various molecular motions that may influence structures of membrane-embedded proteins. Structures stable for a nanosecond may substantially reorganize or even disintegrate over longer periods of time. The reverse, however, is rarely true. If a protein assembly rapidly falls apart, it is highly unlikely that it will eventually adopt a structure similar to the initial one. Disintegration, of course, may be caused by choosing initial conditions for the simulations that are very far from equilibriumif there are large repulsive interactions in the starting configuration any system may rapidly undergo large deformations. Precisely for this reason we made considerable efforts to relieve strains in system 2a before starting unconstrained production trajectories.

Our conclusions about the stability of the M₂ channel are in agreement with the results of Zhong et al. (1998), who studied the same system in a membrane-mimetic lamella of octane. They used a different procedure to construct a structure similar to system 2a and found that the helix bundle was unstable. Furthermore, they also observed that the truncated peptide, composed of only 19 residues, did not form a proper channel. Both results are at variance with simulations reported by Sansom et al. (Sansom et al., 1997; Forrest and Sansom, 1998). They found that a system containing four biprotonated His³⁷ residues was stable and water molecules were able to penetrate the gate. In their simulations, however, the membrane was not represented explicitly but was modeled as a continuum dielectric medium, and several artificial restraints were applied to the protein bundle (Sansom et al., 1997). Recently, they reported a 1-ns simulation of the unprotonated, truncated form of the channel in the explicit membrane (Forrest and Sansom, 1998) and found this structure to be stable. Because their report is very brief, it is not possible at present to determine why their results may differ from the results obtained in this study and in the previous work of Zhong et al.

We further note that only a few discrete states of the system were considered in this work. The actual proton transport also involves a large ensemble of states in which a charge resides on a water molecule in the channel. For mechanism I, it is unlikely that these states are less stable than system 1b, because in this system the charge is located in the most crowded part of the channel, the gate. The situation is different in mechanism II. As the proton transported along the water wire reaches the gate it brings an additional, fifth positively charged species to this region, and thus the system should become even less stable.

Another objective of the simulations was to investigate the arrangement of water molecules inside the channel and their capacity to form hydrogen bonds that assist proton translocation. Transport of protons across membranes along networks of hydrogen bonds involving water molecules and/or protein side chains is quite common in biological systems. It has been postulated to play a role in the functioning of several proteins essential for cellular bioenergetics, such as F₁F₀ ATPases, the photosynthetic reaction center, cytochrome f, and bacteriorhodopsin (Akeson and Deamer, 1990; Baciou and Michel, 1995; Martinez et al., 1996; Kimura et al., 1997; Morgan and Wikstrom, 1994; Pomes et al., 1998). The same mechanism has also been implicated in proton transport through a simple transmembrane protein channel, gramicidin A (Akeson and Deamer, 1990), and involves a single file of water molecules filling the pore. The concept of hydrogen-bonded wires is very attractive because it simply explains how protons can translocate through membranes with high efficiency. The process only requires small displacements of several protons between consecutive hydrogen-bonded donors and acceptors along the chain, instead of slow diffusion of charged species. This idea is essentially the same as in the Grotthuss mechanism of proton conductivity in ice (Agmon, 1995). An important characteristic of proton transport along wires is that the barrier to the transfer of protons from the donor to the acceptor of a hydrogen bond is quite low. In fact, simulations of a model water wire hydrated at both ends showed that the free energy of proton transfer is almost constant everywhere along the wire, which makes proton transfer nearly barrierless and rapid (on a picosecond time scale) (Pomes and Roux, 1998).

The involvement of histidine residues in proton transport via the shuttle mechanism, proposed in mechanism I, is also well known in biology. It has been observed in carbonic anhydrase (Ren et al., 1995; Nair and Christianson, 1991) and postulated for respiratory heme-copper oxidases (Morgan and Wikstrom, 1994). The ability of histidine to tautomerize or change its protonation state in response to the molecular environment is well documented for a variety of enzymes (Ren et al., 1995; Nair and Christianson, 1991; Ash et al., 1997; Frey, 1995). As in proton wire mechanism, the barrier to proton transfer via the shuttle mechanism appears to be low. For human carbonic anhydrase III it was estimated at 1.3 kcal/mol (Ren et al., 1995). Quantum and statistical mechanical studies on model systems, such as NH₄⁺-imidazole-NH₃ (Li et al., 1998b) and HCOO-imidazole-H₂O (Li et al., 1998a) provide a detailed picture of double proton transfer mediated by an imidazole ring. It has been shown that a concerted mechanism, in which a proton is accepted at the N position and another proton is simultaneously donated by the N atom, is unlikely. Instead, proton shuttle proceeds through a double-protonated state of the imidazole ring. If both the proton donor and proton acceptor are properly aligned with the ring, transport is very fast; in a model system the characteristic time for the proton transfer reaction was found to be on the order of 1-2 ps (Li et al., 1998b).

These general considerations of proton transport across membranes form the basis for discussing implications of our simulations for the mechanism of M₂ action. In both mechanisms considered here, it is assumed that proton translocation proceeds, at least in part, through a water wire. Yet the simulations demonstrate (data not shown) that there is no obvious rigid hydrogen-bonded chain of water molecules spanning the channel. This does not contradict the water wire hypothesis, especially because no charged water species were present in the simulated systems, but it points to the highly dynamic nature of the wire. In addition, in contrast to gramicidin, water molecules are not arranged in a single file throughout the channel. Because the channel has a funnel-like shape, many water molecules fill the ends of the pore. In these regions, the orientations of water molecules are influenced by the electrostatic field created by highly polar membrane headgroups and the protein residues surrounding the pore. As a result, the dipole moments of water molecules point, on average, toward the middle of the bilayer. This orientation is not ideal for a water wire. Considering these arguments, it is possible that hydronium ions diffuse through the mouth of the M₂ channel and protons are transported through a proton wire/proton shuttle mechanism only through the central part of the pore.

We further note that water molecules near unprotonated nitrogen atoms of His³⁷ are oriented such that they form linear O-H···N hydrogen bonds with the imidazole ring. This orientation is fairly rigid and relatively long-lived, which is ideal for proton transfer. Similarly, orientations of water molecules around the N-H groups of the ring are also conducive to proton transfer. These arrangements may markedly increase the efficiency of proton translocation, compared with other biological systems in which thermodynamic work needs to be expended to align all of the groups involved in the proton shuttle. For example, proper alignment of donor and acceptor groups in human carbonic anhydrase III requires 11 kcal/mol (Ren et al., 1995).

Once a proton is transferred to a -nitrogen atom and an approximately symmetrical biprotonated HIP state is formed, a question arises: What prevents the backward reaction from occurring? Considering that the characteristic time for proton shuttle in model systems appears to be quite short, it is possible that a proton does move several times between water molecules on the opposite sides of a histidine residue before it is successfully transferred to the intracellular side. Because such movements are rapid, they will not have a substantial effect on the total efficiency of proton transfer (this is not the rate-determining step). Alternatively, there might be some changes in the position of the histidine residue that decrease the probability of the backward reaction but were not captured or identified in the simulations. In fact, rotation of His⁶⁴ involved in proton shuttle in carbonic anhydrase II was observed with the change of pH (Nair and Christianson, 1991), but it is not known whether it plays a role in mediating the reaction.

It is natural to assume that the rate-determining step in mechanism I is tautomerization of His³⁷, which is required to bring the system back to its initial state. Unfortunately, there is no estimate of the rate of this reaction in M₂, and no molecular mechanism of tautomerization has been proposed. Unquestionably, these issues are central to evaluating the likelihood of mechanism I. We only note that a new cycle of proton transport may start even before the previous cycle has been completed. This possibility is implied by the apparent stability of system 2b, containing two HIP residues, and by the fact that if a water molecule is located near a HIE residue its orientation is unaffected by the presence of the N-H group in another histidine residue on the same side of the gate. Clearly, the ability to process two cycles in parallel or nearly in parallel, rather than sequentially, would increase the rate of proton transport.

The nature of the rate-determining step clearly distinguishes mechanism I from mechanism II. In the water wire mechanism, this step is most likely the reorientation of water molecules in the wire, once a proton has been transferred. For a chain of water molecules spanning a membrane, it has been estimated that such a reorientation requires ~500 ps (Marrink et al., 1996). Note that this step is not of concern in mechanism I; after tautomerization, water molecules will spontaneously reorganize to their initial orientations.

Considering the structural simplicity of M₂, compared to most other systems transporting protons across membranes, understanding the mechanism of its action may be of considerable significance for constructing functionalized vesicles. These structures may be very useful in medicine and pharmacology, particularly for controlled drug delivery, and may provide important models in studies of the origin of cellular life. In particular, it may be possible to incorporate into the channel a chromophore capable of ejecting protons in response to illumination (Deamer, 1997) and reengineer the protein such that it will act as a simple proton pump, an essential component of a bioenergetics system.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

The simulations presented here are consistent with the "proton shuttle" mechanism of proton transfer across the closed gate in M₂, but not with the mechanism that requires the opening of the gate and the formation of a continuous water wire. The initial state and two distinct intermediates in the "proton shuttle" mechanism that involve different protonation states of one histidine residue were all found to be stable on the time scales of the simulations. Furthermore, the orientation of water molecules near the gate with respect to the nitrogen atoms in the neighboring histidine residues is conducive to proton transfer across the gate. The possibility that a doubly protonated state of the gate can also participate in proton transfer is not excluded. In contrast, protonation of all four His³⁷ residues, which was proposed to lead to gate opening, appears to cause structural disintegration of the channel due to large, electrostatic repulsion in the gate region. If this repulsion is reduced by allowing only two His³⁷ to be protonated, the gate does not seem to open sufficiently, at least on the time scales of the simulations, to yield enough water penetration for efficient proton transport through a continuous water wire mechanism.

Although these results point to a likely candidate mechanism of proton transfer in M₂, many questions regarding this mechanism remain unanswered. What is the free energy and kinetics of proton transfer between water molecules and the histidines in the gate? What are pK_a values of the  and  nitrogen atoms of the histidines at different stages of the transport process? What is the likelihood of forward versus backward transfer of a proton once it reaches the gate? What is the time needed for histidine tautomerization? What is the role of other residues in the protein in promoting proton transport? If we are to fully understand the mechanism of transport, all of these questions must be answered in subsequent experimental and theoretical studies.