INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS CLASSICAL MD SIMULATION OF... SIMULATIONS OF PROTON TRANSPORT... DISCUSSION REFERENCES

Aqueous proton transport plays an important role in a number of biological functions, such as bioenergetics and cell signaling. Proton transport across biological membranes is facilitated by integral transmembrane proteins, which form ion channels. The interior of channel-forming membrane proteins contains a column of water molecules through which protons and other small ions can diffuse across the membrane (Hille, 1992). Usually, ion channels are formed by large transmembrane proteins that are quite large (such as ATP synthases and bacteriorhodopsin), making studies of their selectivity and regulatory functions difficult at the atomic level.

The influenza A virus M2 protein is a relatively small (97 residue) protein, which plays an essential role in the life cycle of the influenza A virus. The transmembrane domain of the M2 protein consists of 19 residues, and it has been shown that a 25-residue synthetic peptide can form single proton selective channel in lipid bilayers (Duff and Ashley, 1992). The M2 channel is composed of a tetrameric bundle of transmembrane peptides, each having an -helical structure (Duff and Ashley, 1992; Kovacs and Cross, 1997) and can be activated by low pH (Wang et al., 1995). Spectroscopic studies (Kovacs and Cross, 1997; Kukol et al., 1999) indicated that the transmembrane -helical segments of the M2 protein are tilted with respect to the bilayer normal.

Molecular modeling of the M2 channel (Pinto et al., 1997) also suggested that the column of water molecules in the channel is interrupted by a ring of four His³⁷ residues, which was implicated in two possible gating mechanisms (Pinto et al., 1997; Sansom et al., 1997). In one mechanism the protons are transferred through the gate via a "proton-shuttle," which involves a histidine residue. In this scenario an -protonated histidine residue accepts a proton from the extracellular environment, becoming doubly protonated. This step can be followed by the release of one of the protons into the intracellular environment and subsequent relaxation of the histidine into the initial state. In the alternate mechanism, protonation of the histidine residues results in the pore opening, which enhances the proton diffusion across the membrane via shuttling through water molecules in the channel. Electrophysiological studies have also demonstrated that the influenza A virus M2 channel can be blocked by small drugs such as amantadine (Wang et al., 1993).

Due to its relative simplicity, the influenza A virus M2 channel has been extensively studied using computer simulation techniques. Earlier studies were performed in vacuo (Pinto et al., 1997; Sansom et al., 1997), while the later ones (Zhong et al., 1998) attempted to incorporate the effect of membrane environment using a mimetic bilayer (octane slab solvated by water molecules). However, due to the absence of the polar phospholipid headgroups, the mimetic bilayer may not be realistic. To address this deficiency, more recent simulations of the M2 channel were performed in explicit phospholipid (POPC or DMPC) membrane solvated by water molecules (Forrest et al., 2000; Schweighofer and Pohorille, 2000). While these simulations investigated the structural properties of the M2 channel in great detail, the actual proton transport through the water wire spanning the membrane could not be modeled directly using their classical molecular dynamics (MD) approach.

Another, more phenomenological, method to study the problem of proton transport in a channel environment was recently proposed by Schumaker et al. (2000, 2001), who used a combined molecular dynamics and diffusion model to predict the conductance of the gramicidin channel. This approach is based on a number of assumptions regarding the proton transfer mechanism through a single water wire which may not be valid, especially in channels of larger pore size.

Our group recently developed the multi-state empirical valence bond model (MS-EVB) for proton transport (Schmitt and Voth, 1998, 1999a,b), which we have incorporated into biomolecular MD simulations to study the proton transport across the M2 channel and in other biomolecular environments. The MS-EVB potential can successfully reproduce various important features of the excess proton in bulk water, such as the proton hopping rate, the density of vibrational states, and the relative stabilities of the solvated Eigen and Zundel complexes. More recently, the MS-EVB model was applied to study the proton wire dynamics in a model hydrophobic channel environment (Brewer et al., 2001) and to simulate the properties of the solvated proton near the surface of a phospholipid membrane (Smondyrev and Voth, 2002).

The aim of the present work is to simulate the proton transport in the influenza A virus M2 channel, embedded into explicit lipid bilayer (DMPC bilayer solvated by water molecules), using the next generation MS-EVB2 model (Cuma et al., 2001; Day et al., to be published). The paper is organized as follows: first, the details of constant pressure (NPT) MD simulations of the M2 channel in DMPC bilayer will be presented. In the subsequent sections the details of MD simulations using the MS-EVB2 model will be given. Structural and dynamical properties of solvated proton in the realistic M2 channel environment will be presented and discussed. Finally, the effect of the M2 channel conformational flexibility on the proton transport is studied. Conclusions that may be drawn from this work, a consideration of both the model's limitations and the relevance of our results to the gating mechanism, and several future directions will be given in the final section.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS CLASSICAL MD SIMULATION OF... SIMULATIONS OF PROTON TRANSPORT... DISCUSSION REFERENCES

Influenza A virus M2 channel model

Several models for the M2 channel have been proposed by various research groups. In the present work, we adopted a 22-residue model of the M2 protein terminally blocked by Ace and NH₂ groups at the N- and C-terminals, respectively. The sequence of the M2 protein corresponds to the Udorn Influenza A strain: Ace-SSDPLVVAASIIGILHLILWIL-NH₂, which runs from Ser²² to Leu⁴³ residues. Our selection of this particular strain is motivated by the finding that the Udorn strain M2 channels are blocked more efficiently than the channels from the Weybridge strain (Wang et al., 1993). In a lipid bilayer environment, the M2 protein is -helical. The protonation states of the ring of Asp²⁴ residues at the N-terminal and of the channel blocking His³⁷ residues were adopted from previous work (Forrest et al., 2000). Ionization state calculations reveal that the net charge over the four Asp residues is 1 (Forrest et al., 2000), indicating that one of the Asp residues is ionized, while the other three are protonated. It has also been suggested that the protonation of the His³⁷ residues may be responsible for channel opening. The channel is presumed to be in its closed state when all four His residues are deprotonated. Protonation of one or more His³⁷ residues therefore should result in the channel opening to some degree. Computer simulations have also suggested that tetrameric bundle stability may also depend on the protonation of His residues (Schweighofer and Pohorille, 2000). When all four His residues are protonated, their mutual repulsion destabilizes the channel structure.

In the present work we focused on the model when only one His residue is protonated and one Asp residue is ionized. Thus the histidine ring occluding the channel pore in our model contained three -His residues and one His⁺ residue. The neighboring helix to the one containing the protonated His residue contains the ionized Asp residue (i.e., the one neighboring in the clockwise direction when viewing the channel down its axis from the intracellular domain). The initial configuration of the M2 channel was also adopted from previous work (Forrest et al., 2000). This channel model was based on spectroscopic information and was generated using a combination of MD and simulated annealing (SA).

NPT MD simulations of the M2 channel in DMPC bilayer

The influenza A M2 channel was embedded into a lipid bilayer, which contained 114 DMPC molecules hydrated by 3552 water molecules. The transmembrane channel was represented as an all-atom model, and the AMBER force field (Cornell et al., 1995) was used to describe the protein interactions. An AMBER-like united atom model (Smondyrev and Berkowitz, 1999b) was used to model the DMPC molecules. All bonds were constrained using the SHAKE algorithm (Ryckaert et al., 1977) with a tolerance of 10⁴, allowing us to integrate the equations of motion using a 2.0 fs timestep. Water molecules were described using a TIP3P model (Jorgensen et al., 1983). Simulations were performed at constant pressure and temperature (T = 308 K), which is above the main transition temperature (T_m = 293 K) for fully hydrated DMPC bilayers. Temperature and pressure were controlled using the Nose-Hoover thermostat (Hoover, 1985) with relaxation times of 0.2 and 0.5 psec, respectively. Electrostatic interactions were calculated using the smooth particle mesh Ewald (SPME) summation technique (Essmann et al., 1995) with a tolerance of 10⁴. The real-space part of the Ewald sum and short-range van der Waals forces were truncated at 10 Å. Periodic boundary conditions were applied in all directions.

The initial M2/DMPC/water configuration was obtained by putting the M2 channel into a lipid bilayer that contained 128 DMPC molecules (64 in each leaflet). This membrane was constructed such that the average area per DMPC molecule was 60.3 Å², as taken from previous constant pressure simulations of pure DMPC membrane (Smondyrev and Berkowitz, 1999a). A total of 14 DMPC molecules (7 in each leaflet) were removed to avoid bad contacts between the M2 channel and lipid bilayer. The system was equilibrated for 200 ps at constant volume and temperature. The system was then subjected to a 1.0 ns NPT equilibrium simulation using the methods described in the previous paragraph.

MS-EVB2 simulations of proton transport through M2 channel

Simulations of proton transport through the influenza A virus M2 channel were performed using a combination of the classical MD and MS-EVB approaches, which is described in greater detail elsewhere (Schmitt and Voth, 1998, 1999a,b). Here we present only a brief description of essential features of the MS-EVB approach. At each MD step, one constructs a number of possible EVB states. Each of those states consists of a hydronium cation and n water molecules. For example, in the gas phase, the Zundel cation H₅O is represented by two EVB states, while four EVB states are required to represent an Eigen cation H₉O (see Fig. 1). In solution, more states are required to allow for the proton hopping between multiple water molecules. An empirical valence bond potential can then be used to construct an (N × N) Hamiltonian matrix H, such that the ground state potential energy (E₀) is given by:
where c is the ground state eigenvector. The "pivot" hydronium for the EVB state selection process is determined based on the state with the highest probability c from the set of eigenvectors. The off-diagonal elements of the Hamiltonian matrix are defined to reproduce various important properties of the system, such as cluster energies, geometries, and proton transfer barriers. After determining the MS-EVB Hamiltonian matrix and ground state eigenvector, forces can be calculated using the Hellmann-Feynman theorem:
where h_mn(x) are the elements of Hamiltonian matrix and x represents all system coordinates. The interactions between hydronium and surrounding water molecules were modeled here using the MS-EVB2 parameter set (Cuma et al., 2001; Day et al., to be published). In the present work, dynamical proton transfer is not allowed between the hydronium and the His³⁷ residues in the M2 channel. The effects of this process will be studied in future research.

View larger version (9K): [in this window] [in a new window]   FIGURE 1   Schematic pictures of the Eigen cation H9O (A) and the solvated Zundel cation H5O (B).

Initial configurations for the proton transport simulations were prepared by replacing a water molecule inside the M2 channel with hydronium. The coordinates of atoms in the M2/bilayer/water system were taken from an equilibrium NPT trajectory, described in the previous section. An initial equilibration was done by constraining the hydronium oxygen atom and performing a short classical MD simulation on the order of several tens of picoseconds. This was done to adjust the water solution structure around the inserted hydronium molecule. Subsequent production runs were carried out using the MS-EVB approach within the MD. At this stage, proton transfer between water molecules was allowed. These MS-EVB simulations were done in the NVT ensemble, using a Nose-Hoover thermostat (Hoover, 1985) with a thermostat relaxation time of 0.2 ps. The system temperature was the same as in the equilibrium classical MD simulations. Electrostatic interactions were calculated using the Ewald summation technique with a tolerance of 10⁴. The real-space part of the Ewald sum and short-range van der Waals forces were truncated at 10 Å, and periodic boundary conditions were applied in all directions. The equations of motion were integrated using a timestep of 1.0 fs, which is consistent with the fastest vibration modes in the system. A constant electric field was applied along the z-axis (electrophysiological experiments are typically conducted under non-zero transmembrane potential) and a pH gradient applied across the membrane to mimic experimental conditions. The electric field applied in our simulations corresponds to a transmembrane potential of ~100-200 mV.

	CLASSICAL MD SIMULATION OF AN M2 CHANNEL IN A DMPC BILAYER

TOP ABSTRACT INTRODUCTION METHODS CLASSICAL MD SIMULATION OF... SIMULATIONS OF PROTON TRANSPORT... DISCUSSION REFERENCES

Because the initial configurations for the simulations of the proton transfer through the M2 channel were taken from the classical MD trajectory, it was important to check the stability of the membrane and particularly, the stability of the four helix bundle. In Fig. 2 the time evolution of the membrane area over the 1-ns NPT simulation is shown. The average membrane area is 3959 ± 39 Å² and the relative fluctuations of the bilayer area are ~1-2%, which are comparable with surface area fluctuations of pure phospholipid membranes. For example, the average area per headgroup in pure DMPC membrane under the same simulation conditions is 60.3 ± 0.9 Å², as determined from NPT molecular dynamics simulations using the same force field (Smondyrev and Berkowitz, 1999a).

View larger version (19K): [in this window] [in a new window]   FIGURE 2   Time evolution of the membrane area in a classical NPT molecular dynamics simulation of the M2 channel in DMPC bilayer.

The stability of the four-helix bundle in lipid bilayers can also be characterized using various structural properties. The root-mean-square deviations (RMSD) of C backbone alpha carbons are shown in Fig. 3. The reference structure was taken from the start of the constant pressure MD simulation. The RMSD levels rapidly after the initial increase and remains stable over most of the 1.0-ns simulation. No strong drift was observed, which might indicate that the bundle was disintegrating. The steady-state RMSD is in the range of 0.15 nm, which is comparable with the values (within the range of 0.2-0.25 nm) reported in previous M2 simulations of the similar complex (Forrest et al., 2000). All four peptides forming the M2 channel remained -helical during the 1.0-ns simulation, with the exception of at most two residues at the peptide terminals. Greater RMSD fluctuations at the helix termini has been previously reported (Forrest et al., 2000), which is consistent with our observation of the uncoiling in these regions.

View larger version (13K): [in this window] [in a new window]   FIGURE 3   C carbon RMSD versus time for the classical NPT MD simulation of the M2 channel.

Another characteristic of the bundle stability is the time-dependent tilt angle of individual helices. The tilt of the -helical segment can be defined by the angle between the bilayer normal and the vector connecting the C backbone atoms in the Asp²⁴ and Ile³⁹ residues. This definition assumes that the average number of residues per helical turn is 3.6; hence the 15-residue spacing generates a vector that is a reasonable approximation to the direction of the helix axis (Schweighofer and Pohorille, 2000). The time evolution of the tilt angles is shown in Fig. 4 for each of the four helices. The average tilt angle is 38.7 ± 2.9°, in excellent agreement with the recently reported solid-state NMR data 38 ± 3° (Wang et al., 2001). Previously reported experimental values were somewhat lower: 33 ± 3° for both solid-state NMR (Kovacs and Cross, 1997) and 32 ± 6° for FTIR (Kukol et al., 1999). The average tilt angles obtained in simulations of M2 channel in POPC membrane (Forrest et al., 2000) were somewhat smaller (~20°) for the 22-mer B. However, it should be mentioned that larger tilt angles at the end of their simulations were observed. Simulations of the M2 channel in DMPC membrane (Schweighofer and Pohorille, 2000) produced different tilt angles for individual helices ranging from ~10° to ~45°, suggesting that the whole bundle is tilted with respect to the bilayer normal. Our simulations did not show significant spread between the values of the individual helix tilt angles, which may indicate that the bundle does not have a notable overall tilt.

View larger version (32K): [in this window] [in a new window]   FIGURE 4   Time evolution of individual helix tilt angles with respect to the bilayer normal in the classical NPT MD simulation of the M2 channel.

Fig. 5 shows a snapshot of the backbone structure of the tetrameric M2 bundle, taken after a 1.0-ns NPT simulation of the M2/bilayer/water system. It shows a striking similarity to the refined structure obtained from solid-state NMR (Wang et al., 2001). One can see that the tetramer obtained in our simulations has a high degree of symmetry, which was lacking in some previous simulations. The reason for such good agreement with experiment may not be immediately clear, but a number of factors such as simulation conditions, treatment of electrostatic interactions, and the potential model used for the phospholipid membrane may have had a significant effect on the resulting structure.

View larger version (36K): [in this window] [in a new window]   FIGURE 5   Snapshot of a representative structure of the M2 channel as the excess proton (yellow ball), forming a Zundel cation with two water molecules, is shuttling through the ring of histidine residues (shown as ball and sticks). Note that two histidine residues (including the biprotonated histidine) are flipped toward the extracellular part of the membrane, opening, at least partially, the channel pore. Also shown are water molecules 237 and 1305, which participate in the "water wire" through with the proton shuttles.

	SIMULATIONS OF PROTON TRANSPORT THROUGH THE M2 CHANNEL

TOP ABSTRACT INTRODUCTION METHODS CLASSICAL MD SIMULATION OF... SIMULATIONS OF PROTON TRANSPORT... DISCUSSION REFERENCES

Initial conditions

We performed a total of seven simulations using the MS-EVB MD methodology, which differed from each other in initial simulation configuration only. In each simulation, a configuration from the equilibrium classical NPT MD trajectory of the M2 channel in DMPC bilayer was selected and an excess proton was placed near the N-terminus of the M2 bundle inside the channel. The details of the seven initial conditions are shown in Table 1. The question of how the proton enters the channel was not investigated in this work, as the primary goal of the present simulations was to study the proton transfer inside the channel environment. Channel entrance will be a topic of future research.

In simulations 1, 5, and 7 the excess proton passed through the M2 channel and exited into the bulk water on the intracellular side of the bilayer. In simulations 2 and 3, the excess proton diffused out of the channel into the bulk water on the extracellular side, while in simulation 4, the excess proton diffused toward the helix bundle "wall" of the M2 channel. The same initial configuration (before proton insertion) was used as a starting point in simulations 2, 3, and 4, but the excess proton was placed at three different points to ensure that the observed behavior is not an artifact of the initial conditions. In simulation 6, the excess proton drifted toward the channel "wall" and became immobilized without moving through the M2 channel. Our conclusions regarding the excess proton behavior in simulations 4 and 6 are based on the observation that the proton moved away from the channel axis and remained in the same position over the last several tens of picoseconds of simulation.

In the following discussions, only the properties of the systems that exhibited pronounced mobility of the proton along the channel axis (namely simulations 1, 5, and 7) will be analyzed. It must be stressed that the MS-EVB2 simulations were of limited duration because of their computational cost (e.g., <1 ns). Simulations 4 and 6, where the proton became "trapped," might well have exhibited eventual transport through the channel had they been longer. Furthermore, simulations 2 and 3 might have eventually seen the excess proton re-enter the channel due to the favorable transmembrane potential.

Proton mobility in the M2 channel

Due to the delocalized nature of the excess proton, we define a suitable coordinate to describe the position of the "excess proton" at each step. We can adopt the "center of excess charge" (CEC) (Cuma et al., 2001) as a coordinate, which is given by:

(1)

where c_i is the ith EVB state component of the ground-state MS-EVB eigenvector and r_i(t) is the center of charge position vector of the hydronium in the ith EVB state. This center of excess charge coordinate can be used to determine proton diffusion rates through the M2 channel. Because the M2 channel axis is parallel to the z-axis, one can define the diffusion constant, D_z, using the Einstein relation:

(2)

The diffusion constants obtained from the trajectories spanning the length of the M2 channel in simulations 1 and 5 were (1.88 ± 0.04) × 10⁹ m²/s and (1.12 ± 0.02) × 10⁹ m²/s, respectively, compared to a value of (3.5 ± 0.9) × 10⁹ m²/s for the MS-EVB2 model in bulk water. These values are only two to three times lower than the proton diffusion rates in bulk water defined from the MD simulations. The value of the diffusion constant D_z obtained from the trajectory in simulation 6 was (0.62 ± 0.01) × 10⁹ m²/s, which reflects the fact that the proton finally became immobilized near the channel wall. The value for the diffusion constant obtained in simulation 7 was (3.77 ± 0.03) × 10⁹ m²/s, which is comparable to the diffusion constant for the proton in bulk water. It should be noted that in simulation 7 the proton was initially placed deeper in the channel compared to the other simulations. It is further possible that a favorable channel conformation may have also facilitated the faster proton transfer through the channel.

It was found that the diffusion rates vary notably as the proton progressed through the channel. In Fig. 6 we show the time evolution of the center of excess charge position along the channel axis for both simulations 1 (top plot) and 5 (bottom plot). Each of these trajectories was split into 150-ps segments and diffusion constants were then determined separately for each one. The resulting values are listed in Table 2. These results indicate that the diffusion rates decrease as the proton moves toward the middle of the channel and then increase again as the proton is being released into the bulk water. Interestingly, the axial diffusion coefficients obtained in the present work are similar to the diffusion constants found in simulations of water in model channel environments of comparable radii (Brewer et al., 2001).

View larger version (27K): [in this window] [in a new window]   FIGURE 6   Time evolution of the projection of the center of excess charge (Eq. 1) on the channel axis in simulations 1 (plot A) and 5 (plot B). Vertical lines separate the trajectory into segments that were used to calculate diffusion coefficients.

View this table: [in this window] [in a new window]   TABLE 2   Diffusion coefficients (×109 m2/s) along the channel axis in simulations 1 and 5, calculated for different 150-ps parts of the MD trajectory (see Fig. 6 for definitions of different regions)

MS-EVB amplitude analysis

The probability distribution of the largest EVB amplitude, c, averaged over the complete trajectories for the hydronium in simulations 1 and 5 are shown in Fig. 7. This provides a convenient way to characterize the structure of the solvated hydronium (Schmitt and Voth, 1999a). In bulk water, a value of c ~ 0.7 corresponds to an Eigen complex, whereas in the Zundel complex, two EVB states share the excess proton almost equally, resulting in an amplitude c ~ 0.5. The probability of Zundel cation formation is slightly higher in these two M2 simulations compared to the bulk water results (Schmitt and Voth, 1999a). The same trend was predicted from simulations of the excess proton in model hydrophobic channels of various radii (Brewer et al., 2001).

View larger version (23K): [in this window] [in a new window]   FIGURE 7   Probability of the largest EVB amplitude c for the excess proton in the M2 channel (in simulations 1, 5, and 6) and in bulk water (shown for a reference).

The EVB coefficient distribution obtained from simulation 6 is rather different. As the proton becomes immobilized near the channel wall, the probability of forming an Eigen cation decreases significantly, which may be due to both energetic reasons and steric interactions (Fig. 7). At the same time, the number of possible EVB states in the first solvation shell decreases and the probability of finding a complex with an amplitude ~0.8 increases. Similar results were obtained in simulations of an excess proton near the surface of a phospholipid membrane as the excess proton diffused toward the bilayer interface (Smondyrev and Voth, submitted for publication).

The amplitude distributions were also calculated for smaller segments of the total trajectory (as in the previous subsection for diffusion constants). No substantial difference was observed, despite varying average pore radii for each segment. This result is somewhat different from previous findings in which the amplitude distributions were more sensitive to the pore radius (Brewer et al., 2001). It is possible that this discrepancy highlights the effect of the polar channel environment and surrounding membrane on the structure of protonated water inside the M2 channel. We also calculated the free energy profiles associated with the proton transfer "reaction coordinate" (Cuma et al., 2001), which can be conveniently defined as

(3)

where i and j are the indices of the EVB states with the two largest amplitudes. The reaction coordinate is zero in the case of a Zundel cation and becomes q_react ~ 0.5 for an Eigen complex (in bulk water). Fig. 8 shows the free energy profile for the proton transfer in simulations 1, 5, and 6. The free energy profile for proton transfer in bulk water is presented on the same figure as a reference. One can see that the barrier to formation of a Zundel complex is lower by ~0.2 kcal/mol compared to bulk water. A similar trend was also observed in simulations of water wires in model hydrophobic channel environments as the radius of the channel was decreased (Brewer et al., 2001).

View larger version (25K): [in this window] [in a new window]   FIGURE 8   Free energy profiles as a function of reaction coordinate qreact (Eq. 3) for the hydronium in the M2 channel and in bulk water.

Pore radius analysis and proton passage through the HIS gate

We now analyze the relationship between channel conductance and the geometric properties of the pore of the M2 channel. Pore radius profiles were calculated using the HOLE program (Smart et al., 1997). First, we analyzed the pore profiles for the initial configurations used in simulations 1-6. Interestingly, we did not observe any significant differences in the regions corresponding to the initial positions of the excess proton. Thus, the diffusion of the excess proton back into the extracellular bulk water in simulations 2-4 cannot be explained using a simple geometric argument.

In Fig. 9 we plot the pore radius profiles for four different positions of the excess proton along the M2 channel, shown by vertical arrows. The top plot shows the pore radius profile evolving as the excess proton passes through the M2 channel in simulation 1, while the bottom plot shows the same characteristics for simulation 5. The numbers with vertical arrows in each case correspond to the instantaneous positions of the excess proton in the channel. Similar trends have emerged for both cases. The excess proton passes through the ring of histidine residues (which may be responsible for channel gating) in ~100 ps or less. The approximate location of the ring of four histidine residues is Z = 7-8 Å in our simulations. The start of the proton diffusion through the ring of histidine residues (from position 3 to position 4, see Fig. 9) is preceded by the opening of the channel pore in the vicinity of the histidine ring (see plot 3 on both graphs, which is shown as a bold solid line). The channel pore then closes once the proton passes through the histidine gate (see plot 4, shown as thin solid lines). These results suggest that concerted changes in channel conformation may be required for the successful passage of an excess proton through the M2 influenza A channel.

View larger version (17K): [in this window] [in a new window]   FIGURE 9   Pore radius profiles calculated at different instances along the simulation trajectory in runs 1 (plot A) and 5 (plot B). Instantaneous positions of the center of excess charge corresponding to each of the profiles are shown with arrows. The His37 residues are located at ~z = 7-8 Å from the bilayer center (z = 0 Å).

It should also be mentioned that the passage of an excess proton through the ring of histidine residues is not accomplished by the simple diffusion of the hydronium ion. From visual inspection of two trajectories in which the proton passes through the channel (simulations 1 and 5), we concluded that the following mechanism facilitates the proton passage through the gate: the proton moves through the histidine ring by hopping between two water molecules that approach the ring of histidine residues from the intra and extracellular regions of the channel, respectively. When the proton hopping across the gate occurs, the two water molecules involved and the excess proton form a Zundel complex that spans the histidine ring region (see Fig. 5). It must be stated here, however, that as more detailed experimental structural data on the channel become available, this proposed mechanism may well need to be revisited.

We now examine how the changes in channel conformation affect the diffusion of the excess proton through the M2 channel. Fig. 5 shows a snapshot of the M2 channel taken from an MD trajectory (simulation 1) in which the excess proton passed through the ring of four histidine residues. It shows that the opening of the channel pore may involve the flipping of the histidine residues, which would otherwise block the channel. Similar behavior was observed in simulation 5. These results also suggest that efficient proton transfer through the M2 channel may require concerted changes in channel conformation. To test this hypothesis, we conducted another simulation of proton diffusion along the M2 channel in which the backbone atoms of the channel were constrained in space, but the amino acid residue side chains were free to move. The starting configuration was taken from 90 ps into simulation 1 (for details, see Table 1) such that the center of the excess charge was initially at Z_CEC ~ 7.0 Å. In Fig. 10, the CEC position along the z-axis is shown as a function of time for simulations of the proton transfer in both the unconstrained and the "frozen" M2 channels. In the simulation in which the M2 channel was frozen, the excess proton remained near its initial position for 400 ps, while in the simulation in which the M2 channel was unfrozen, the excess proton diffused the length of the channel in approximately the same amount of time. Although these results are only for single trajectories, their qualitative difference appears to be very important in illustrating the dynamical effect of the channel backbone fluctuations.

View larger version (24K): [in this window] [in a new window]   FIGURE 10   Time evolution of the projection of the center of excess charge on the channel axis in simulation 1 and a simulation with a frozen channel. The latter simulation started from the configuration corresponding to 90 ps of simulation 1.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS CLASSICAL MD SIMULATION OF... SIMULATIONS OF PROTON TRANSPORT... DISCUSSION REFERENCES

In this paper we have presented the results of MD simulations of excess proton diffusion in the influenza A virus M2 channel. While a number of classical molecular dynamics simulation studies of the M2 channel have been presented in recent years, the present work includes the possibility of a proton hopping between water molecules using the MS-EVB methodology. Recent simulations of proton transport in a smooth, cylindrical, hydrophobic channel (Brewer et al., 2001) using the same approach suggested that changes in the structure of protonated complexes may be explained simply on geometric grounds. For example, the Zundel complex is stabilized with respect to the Eigen as the channel radius decreases. While a similar trend was observed in the present work, the effect was not as pronounced as in the simplified channel environment. This may be attributed to the effect of the polar environment inside the M2 channel. At the same time, the diffusion coefficients were found to be very different inside the channel compared to regions close to the channel termini. It was also found that the proton (even when initially placed inside the channel) may diffuse back into the bulk water layer even in the presence of a transmembrane electric field. It is possible that local electric fields may be responsible for this observed behavior and that the inclusion of the explicit membrane environment may also play a role in the process. Because the main goal of the present work was to study the structural and dynamical properties of the excess proton inside the channel, we have not yet investigated how the proton actually enters the channel in detail. Some important questions such as the protonation of Asp residues at the C-terminus of the M2 channel and its dependence on the channel conformation may arise in this context.

We have focused on a system consisting of three -protonated and one biprotonated histidine residues. The dynamical protonation/deprotonation of the histidine residues was not allowed. While this did not allow us to probe the postulated mechanism whereby the proton is transferred via a "proton shuttle" involving a histidine residue, it was found that the proton can still diffuse through the histidine ring via a water wire. It was also discovered that the changes in the channel conformation may have a profound effect on the proton transport, and that the proton is capable of transferring through a confined space (such as the ring of histidine residues in the M2 channel), which is impenetrable to water molecules. These results suggest that, in the case of proton channels, a geometric criterion that would otherwise be reasonable for predicting ion conductance properties may overestimate the energy barrier for proton conductance.

It was found that the passage of the excess proton through the channel may be accompanied by concerted changes of the channel conformation. As the proton approaches the ring of histidine residues, the pore of the channel opens, possibly resulting in the flipping of some of the histidine residues. In two events in which the proton passed through the channel, the biprotonated histidine was involved in the flipping motion. As the proton passed through the histidine gate, the pore of the channel closed behind it. This observed behavior may be attributed to the repulsion between the positively charged hydronium and the positively charged biprotonated histidine residue. This presents an interesting possibility for a second mechanism of the proton shuttle involving the histidine residue, involving two protons. In this mechanism, a proton is first bound to one of the histidines (making it biprotonated). As a second proton approaches the gate, the histidine residue flipping motion opens the channel and releases the first proton into the intracellular environment. The second proton may then protonate one of the histidine residues in turn, allowing the process to repeat. While our results suggest that such mechanism may not be required to explain the proton conductance, it could become important in a situation in which the extracellular environment is acidified and there is a higher concentration of free protons. These possibilities will be explored in future research.

Admittedly, there are a number of limitations in the simulation methods used in the present paper. While the use of the MS-EVB approach has allowed us to perform simulations on a time scale sufficiently long to observe proton diffusion through the M2 channel, the computational cost of such simulations is still relatively high. Thus, our data analysis could utilize only a limited number of trajectories. From this perspective, the present results should be viewed qualitatively and be used primarily to define the range of values associated with a particular process. At the same time, we observed a number of features of the proton transfer through a realistic channel environment that allows us to make certain preliminary conclusions regarding the possible proton transport mechanisms involved. In this sense, our results are the first of their kind.

The simulations presented in this work may also be sensitive to a number of factors, for example the protonation of key residues, such as Asp²⁴ and His³⁷. The proximity of the histidine residues forming the four-residue ring further complicates the treatment of ionization states of these residues because such calculations are sensitive to channel conformation.

Future simulations from our group will address some of these issues. For example, more detailed studies of the gating mechanism and the role of histidine residue protonation/deprotonation in the proton transport are already within reach. Future simulations will also focus on the role of histidine ionization states in the gate closing and in the proton conductance. The MS-EVB2 model will certainly provide a computationally efficient way to study a number of mechanisms that are not directly accessible using conventional MD techniques. Results from simulations using the MS-EVB methodology may also be used as input for other computational modeling methods. For example, calculated diffusion constants, combined with free energy profiles for the proton conductance (which can be obtained using umbrella sampling techniques or continuum electrostatics methodologies), may be used to calculate the conductance of proton channels. Work is currently underway in this direction as well.