Evaluating how rimantadines control the proton gating of the influenza A M2-proton port via allosteric binding outside of the M2-channel: MD simulations

Abstract

In order to understand how rimantadine (RMT) inhibits the proton conductance in the influenza A M2 channel via the recently proposed “allosteric mechanism”, molecular dynamics simulations were applied to the M2-tetrameric protein with four RMTs bound outside the channel at the three protonation states: the 0H-closed, 1H-intermediate and 3H-open situations. In the 0H-closed state, a narrow channel with the RMT-Asp44-Trp41 H-bond network was formed, therefore the water penetration through the channel was completely blocked. The Trp41-Asp44 interaction was absent in the 1H-intermediate state, whilst the binding of RMT to Asp44 remained, which resulted in a weakened helix-helix packing, therefore the channel was partially prevented. In the 3H-open state it was found that the electrostatic repulsion from the three charged His37 residues allowed the Trp41 gate to open, permitting water to penetrate through the channel. This agreed well with the potential of the means force which is in the following order: 0H > 1H > 3H.

Keywords::

• H1N1 strain
• protonation state
• potential of means force
• proton transport

Introduction

The influenza A viruses, both the impending isolates such as H5N1 and the recent new strain (H1N1), are a vital and emergent problem of the global flu pandemic [1–2]. One of the important targets for disruption of the replication process is based on the fact that the integral M2 membrane protein pH gated channel is required in the early and late stages of infection [3–5]. Amantadine (AMT) and rimantadine (RMT) are the first effective drugs licensed for influenza treatment and function as antivirals via inhibition of the function of the M2 proton transporter. Besides pore blocking [6–7], which is a conventional mechanism, a new novel allosteric mechanism where four RMTs were found to bind outside the pore, was recently proposed based on an NMR study [8]. This highlights the need to understand how RMT controls the gating residue of the M2-channel. The present study aims to provide information at the molecular level for the development of more effective drugs.

The M2 protein of the influenza A virus functions as a proton selective channel that is activated at the low pH of the endosome after endocytosis of the virus. In the early stages of infection, the M2 proton transporter leads to acidification of the viral interior from the acidic endosomal compartment following cellular endocytosis and this is required for the unpacking of the viral genome via release of viral nucleoproteins prior to viral replication. Whilst in the later stages of infection, the M2 proton port is required to transport protons from the transgolgi membrane to the host cell cytoplasm to equalise the pH and so prevent premature conformational rearrangement of the newly synthesised haemagglutinin during transport to the host cell surface membrane [5]. The M2 protein is structurally a homotetramer which is stabilised in part by the disulphide bridges between the N terminal domain cysteine residues near the membrane, with each 97 amino acid monomer being composed of a 24-residue N-terminal extracellular domain, an α-helix single transmembrane (TM) domain of 19 residues and a 54-residue cytoplasmic domain [9]. The TM helix of M2 (M2-TM) spans the hydrophobic region of the membrane and includes a few hydrophilic residues at either end. The M2-TM forms tetrameric bundles and binds adamantane-based drugs, acting as the full M2 protein does, both in the micelles [10] and in the lipid bilayers [11–14].

Although the model structures of M2-TM are experimentally [12–17] and computationally available [18–29], high resolution structures have only recently been determined [7,8]. X-ray diffraction analysis has been used to solve the M2-TM structures under neutral and low pH conditions. The crystal structure revealed that a single AMT molecule in the pore of the channel was surrounded by Val27, Ala30, Ser31 and Gly34 [27,28]. In contrast, Schnell and Chou reported an NMR structure for the peptide spanning residues 18–60 in detergent micelles at high pH [8], where four RMT molecules were bound outside of the protein helix, facing the lipid bilayer and located in the membrane environment at the end of the helix towards the cytoplasmic face of the channel (Figure 1). Drug binding includes interactions with Leu40-Arg45 residues and, in particular the hydrogen bond formation between the –NH₂ group of RMT and the −COO⁻ group of Asp44 appears to be important [8]. This newly proposed channel inhibition model is known as the allosteric mechanism [29]. The novel external drug binding was suggested to stabilise the closed state, thus making the channel more rigid and difficult to rearrange the four helices to allow the channel to open. Although the NMR structure of the M2-TM region seems to be very well solved, the drug-binding site around Asp44 is not well defined and RMT adopts too many possible conformations. Even within a single NMR structure, each of the four RMT molecules displays a noticeable variability in the hydrogen bond and hydrophobic interactions with the M2 binding residues [30].

Figure 1.  The initial structure of the tetrameric M2 protein complexed with four rimantadines (RMTs) bound outside the channel in the pre-equilibrated lipid bilayer-water pieces where L represents the distance along a channel axis starting from the extracellular site. The structure of the His37 tetrad has been coloured orange. Close view of the RMT-Asp44-Trp41 hydrogen bond network with definitions of the d1–d4 distances are also shown.

In the present study, molecular dynamics (MD) simulations were carried out for the tetrameric M2 protein channel complexed with four RMTs at the three protonation states of the His37 tetrad, corresponding to: the closed (0H), intermediate (1H) and open (3H) channel conformations. The main goal was to provide detailed information at the molecular level of how RMT inhibits the M2 proton conductance through the allosteric mechanism in terms of drug-target interactions, the conformation of the Trp41 channel gate, the water density and the potential of mean force (PMF) of water permeation along the channel axis.

Methodology

Preparation of the simulated systems

The NMR structures of the M2 tetrameric helix (residues 18–60) with the C-terminal base complexed with four RMTs, as taken from the Protein Data Bank [PDB entry code: 2RLF [8], were used as the initial coordinates for the MD simulations. The simulations were carried out for the three different protonation states of the His37 residues, where all four His37 residues were considered neutral entities for the closed channel, one was positively charged for the intermediate channel and three of them were positively charged for the open channels, as a function of pH [8,12,31,32]. For simplicity, the simulations were denoted hereafter as 0H, 1H and 3H for the non-, mono- and triple-protonated states, respectively. In all the complexes, the protonation state of the ionisable amino acid residues were assigned in the following manner: side chain of Asp(24,44) and Glu(56) were considered to be negatively charged while that of Arg(45,53) was treated as positively charged. Therefore, the interhelical salt bridge between the Asp44 and Arg45 was well formed [8]. Regarding drug interactions, the relatively high pK_a values of 10.4 [33] made the side chain of the RMT have a protonated form (R–NH₃⁺) [28] and consequently this –NH₃⁺ group was in contact with the –COO⁻ group of Asp44 at the allosteric binding region [8]. Each system was separately built according to the designed protonation state of the His37 tetrad with the RMT inhibitor bound and then inserted into a pre-equilibrated lipid bilayer, initially made up of 77 molecules of 1-palmitoyl-2-oleoyl-sn-glycerol-3-phosphatidylcholine (POPC) lipid [34] embedded in 3760 molecules of FLEXSPC water [35]. The simulated systems were neutralised by counterions and the solvated box dimensions were set to 70 Å × 70 Å × 80 Å.

MD simulations

All calculations were performed using the GROMACS 3.2.1 package [36] with the GROMACS force field [37]. The inhibitor’s topology file was taken from our previous studies [28, 38]. The whole structure was optimised by the steepest descent minimisations. In the next stage, the system was equilibrated for 0.5 ns with position restraints on the protein atoms to improve the packing of the lipid bilayer around the protein, followed by a free MD simulation of 16 ns during which time the structural coordinates were saved every 0.5 ps.

For all the simulations, the periodic boundary condition with fixed pressure P, temperature T, and number of atoms N (NPT ensemble) was employed. The LINCS algorithm [39] was applied to constrain the bond lengths and the angles involving hydrogen atoms and a 2 fs time step was used. The systems were coupled separately to a Berendsen temperature bath [40] at 310 K using a coupling constant τ_T = 0.1 ps. The pressure (1 bar) was kept constant by semi-isotropic coupling of the system to a Berendsen pressure bath [40]. The long-range interactions were restricted to within the twin-range cutoffs of 1.2 nm for both van der Waals interactions and electrostatic interactions, computed using the Particle Mesh Ewald (PME) algorithm [41]. The analysis phase was from 8 ns to 16 ns, in which the convergences of energies, temperature, pressure and global root mean square displacement (RMSD) were used to verify the equilibrium of the system.

Free energy barrier for water permeation

The energy barrier of the water transport throughout the M2 channel was calculated with an analogous treatment according to Raschke and Levitt [43]. The PMF of the water permeation along the channel axis can be derived from the water density plot according to:

1

where R is the gas constant, T is the temperature in Kelvin and w(r) is the probability density of the water distribution in a spherical region within a radius of r. Here, the probability density is given as:

2

where ρ_r and ρ_total are the local and total water densities, respectively. In the water density calculation, r is extended to 80 Å to ensure that water molecules in the aqueous zone (bulk water) are taken into account. ρ_total represents the bulk water, which presumably has a total probability equal or near to 1. By employing this assumption, ρ_total has the maximum value of ρ_r .

Results and discussion

System stability

The RMSD for the M2-RMT complexes in the 0H, 1H and 3H states were calculated and used for monitoring the stability of the systems. The calculations were carried out relative to the initial structure (heavy atoms only) of the M2-RMT complexes and the results are shown in Figure 2. The RMSD plots indicated that all systems reached equilibrium after the 8 ns simulations.

Figure 2.  The root mean square displacements (RMSD) for the M2-RMT complexes in the 0H, 1H and 3H states.

Protein-protein and drug-protein interactions

The hydrogen bond network of RMT-Asp44-Trp41 (Figure 1) has been proposed to be a vital determinant in the inhibitory mechanism of RMT when bound outside of the M2 channel [8]. To extract such information from the MD simulations, the percentage hydrogen bonding of the important central residue (Asp44) with the RMT inhibitor (RMT-Asp44), and the channel gating residue Trp41 (Asp44-Trp41) at the adjacent subunit, were evaluated based on the criteria of a proton donor-acceptor distance ≤ 3.5 Å and a donor-H-acceptor bond angle of ≥ 120° (Table 1). The RMT-Asp44-Trp41 hydrogen bond network was described in terms of the four distances, d1-d4 (Figure 1), which represent the RMT-Asp44 (d1-d2) and Asp44-Trp41 (d3-d4) pair interactions (Figure 3). The plots were separately determined for the four subunits (I–IV) of the 0H, 1H and 3H protonated states of the M2 protein.

Table 1.  Percentage of hydrogen bonds in the three (0H, 1H and 3H) protonation states between (i) the ammonium group of RMT and M2 residues and (ii) the Asp44 and Trp41 gating residues.

Interaction	% Hydrogen bond occupancy
(i) RMT-M2:
0H state
d1: (RMT-II)NH...OD1(Asp44-II)	31
d2: (RMT-II)NH...OD2(Asp44-II)	14
d1: (RMT-III)NH...OD1(Asp44-III)	10
d2: (RMT-III)NH...OD2(Asp44-III)	9
(RMT-III)NH...O(Trp41-IV)	11
(RMT-IV)NH...O(Phe48-I)	44
1H state
d1: (RMT-I)NH...OD1(Asp44-I)	11
(RMT-I)NH...O(Asp44-I)	25
(RMT-I)NH...O(Phe47-I)	49
d1: (RMT-II)NH...OD1(Asp44-II)	12
d2: (RMT-II)NH...OD2(Asp44-II)	16
(RMT-II)NH...O(Asp44-II)	47
(RMT-III)NH...O(His57-IV)	58
3H state
(RMT-III)NH... O(Arg45-IV)	26
(RMT-IV)NH... O(Asp44-IV)	5
(RMT-IV)NH... O(Leu46-IV)	10
(ii) Asp44-Trp41*:
0H state
d3: (Asp44-II)OD1…NE1(Trp41-III)	80

*No hydrogen bond was found for the 1H and 3H states

Figure 3.  The distributions of the d1–d4 distances, as defined in Figure 1, for the three simulated systems: 0H, 1H and 3H, at the subunits I–IV of the M2 channel.

The complete RMT-Asp44-Trp41 hydrogen bond network was only present in the 0H-closed state, where the RMT-Asp44 and Asp44-Trp41 interactions were strongly and clearly detected. The sharp and narrow peaks for subunit II at the d1 and d2 distances of ~3.5 Å (Figure 3B) indicated the two hydrogen bonds between the RMT and Asp44, with occupations of 31% and 14% (Table 1). At the same time, a strong and more secure hydrogen bond was formed between the Asp44 of subunit II and the Trp41 of subunit III at d3 of about 3.5 Å (Figure 3B), with an occupancy of 80% (Table 1). The presence of this hydrogen bond network among RMT, Asp44 and Trp41 indicated that close contact between the RMT ammonium group and the Asp44 carboxylate group could bring the bulky Trp41 indole ring into van der Waals contact to form a more stable Trp41 channel gate. In addition, the RMT-Asp44 interactions at subunit III (Figure 3C), indicated by the sharp peaks of the d1 and d2 at ~3.5 Å, also facilitate the formation of the closed conformation of the M2 channel.

In the case of the intermediate channel conformation seen with the mono-protonated (1H) state, the shortest hydrogen bond distances were found for the RMT-Asp44 (d1-d2) interactions at the three subunits, I, II and III (Figures 3E–3G), whilst the Asp44-Trp41 (d3-d4) distances were between ~7 and ~10 Å (Figures 3E–3H). Therefore, only the rather weak interactions between the RMT and Asp44 were found (<20% occupation, Table 1). The results indicated that Asp44 in the 1H state cannot play a role as the centre of the hydrogen bond network between the inhibitor and the gating residue Trp41. The situation is totally different for the higher protonated state, 3H, where almost no interaction amongst these three molecules was detected. This can be seen as the broad peak positions at longer distances of d1-d4 (Figures 3I–3L), without any hydrogen bond occupancy (Table 1). This observation implies that the strong electrostatic repulsion among the three protonated His37 imidazole rings in the 3H state destabilises the helix–helix packing, leading to the conformational rearrangement to break interactions between Trp41 and Asp44, and so allow the Trp41 gate to flip open (more details in the “Conformation of the Trp41 gating residue” section) and reduces the ability of RMT to bind to Asp44.

Water density across the M2 channel

To determine the inhibition ability of RMT at the allosteric site of the M2 channel, the water density profile as a function of the distance (L) along the pore-axis of the channel, starting from the N-terminal site, was evaluated In addition, the distribution patterns of the His37 proton selectivity and the Trp41 gate opening residues were analysed in order to illustrate the movement of these critical residues (Figure 4).

Figure 4.  The water densities and the distribution patterns of the His37 selectivity and Trp41 gating residue positions in the free M2 protein and M2-RMTs complex.

In the 0H-closed state of the M2 channel, a zero water density was observed in the channel pore at 35 Å < L < 45 Å for both the free and the RMT-bound forms. This indicated that water cannot penetrate through the channel when all four His37 residues are uncharged. With respect to the intermediate state (1H), although water entry was fully inhibited at ~30 Å in the free form (Figure 4), a degree of water density was observed in the area of the His37 and Trp41 residues. Interestingly, in contrast a zero water density was found across this area in the M2-RMT complex. The results detected for both 0H and 1H states lead us to conclude that the four RMTs bind to the external pore side (lipid facing pocket) in the systems with neutral (pH = 7.5) or one charged (pH = 7) histidine residues and this not only helps to stabilise the M2-channel in its closed conformation but also blocks water penetration.

In contrast to the 0H and 1H protonation states discussed above, a zero water density in the RMT bound M2 channels was not observed in the 3H state (Figure 4). Although the water density in the M2-RMT complex was significantly lower than that of the free channel, the four RMT drugs bound outside the pore did not inhibit the water transport. This observation is in good agreement with a previous electrophysiological study where the binding of RMT outside of the M2 channel was found to not be the primary site for pharmacological inhibition or proton transport [30]. In contrast, water penetration was completely prevented in the 3H state when RMT was bound inside the M2 pore [28]. A clear picture of water passing through the pH sensor His37 and proton gate Trp41 residues in the 3H state can be ascribed as follows. At pH ≤ 6, the His37 imidazole ring is protonated leading to destabilisation of the helix-helix packing due to the strong electrostatic repulsion. This conformational rearrangement leads to a breaking of the interaction between Trp41 and Asp44 residues (see Figure 3 and details discussed in the “Protein-protein and drug-protein interactions” section) and allows the gate to flip to an open form. Therefore, adding of RMT to the channel at pH ≤ 6, when the channel is in the open form, cannot facilitate the tetrad Asp44 to play a role as a locking residue for the M2 channel.

The above notion agrees very well with the protein-protein and drug-protein interactions, in terms of the RMT-Asp44-Trp41 hydrogen bond network, where the complete hydrogen bond network was only found in the 0H state, and only moderate RMT-Asp44 hydrogen bonds were detected in the 1H state, while those interactions were completely lost for the highest (3H) protonated state.

Conformation of the Trp41 gating residue

Figure 5 shows the CA-CB-CG-CD2 torsion angle of the indole ring of the gating residue Trp41 in the three different channels. In the 0H and 1H states, the orientations of the torsional angles of the Trp41 indole rings of the four M2 subunits are considerably similar. All rotations occur within the range of 30° to 120° with a maximum at 90°. The implication is that, at the low (0H and 1H) protonation states, the tryptophans mostly lie in a configuration that is almost perpendicular to the channel axis, i.e. the channel is closed by the four indole rings, explaining why water cannot pass through the 0H and 1H channels as discussed previously. In contrast, considerable variation in the torsion angles of the four tryptophan sidechains was found between the four subunits in the 3H state; namely at −110°, 110°, 130° and 170° for subunits I–IV, respectively. The dramatic rotation by ~180° in the Trp41 of M2 subunit-I and the moderate rotation in the Trp41 residues of the other three M2 subunits on changing from the 0H and 1H state to the 3H state indicates the transformation of the channel from the closed (0H and 1H) to the open (3H) conformation, supported by the water density detected along the pore (Figure 5).

Figure 5.  The rotational angle defined by CA, CB, CG and CD2 of the M2 Trp41 gate when M2 is complexed with the four RMTs outside the pore at different protonation states (0H, 1H and 3H). I–IV denote the four M2 subunits. Here, the vertical grey solid line (at 75°) represents the tortional angle of the closed conformation of the Trp41 indole ring obtained from NMR structure.

Affinity of water permeation

To investigate how difficult it may be for water molecules to move through the M2 channel in the free and RMT complexed forms of both four RMTs binding outside- (M2-RMT_out) and a RMT binding inside- (M2-RMT_in) mediate in which the data of M2-RMT_in was taken from our previous study [28], the PMFs of water permeation along the channel axis (L) starting from the N-terminal site were evaluated and compared in Figure 6.

Figure 6.  The potential of mean force (PMF) of the excess proton along the pore channel axis (L) starting from the N-terminal site for the three protonation states, 0H, 1H and 3H of the free M2 (dashed line) and RMT complexes where the four RMTs were located outside (black solid line, M2-RMT_out) and inside the ion channel (grey solid line, M2-RMT_in).

It is clearly seen that the energy barriers between the M2-RMT_out complex and the free M2 protein for the 0H and 3H states are not notably different. However, this is not the case for the 1H state, where the maximum energy barrier of ~7.5 kcal·mol⁻¹ for the complex at L ~36 Å and ~43 Å is significantly higher than that of the free form (~4 kcal·mol⁻¹ at L ~32 Å). In addition, the energy barriers for these free and M2-RMT_out complexed forms were found in the following order: 0H > 1H > 3H. In contrast, the PMFs for M2-RMT_in observed at L~40, 42, 40 Å of 0H, 1H and 3H states, respectively, are significantly higher than those of the free and M2-RMT_out complexes for all the protonation staes.

In the free form, the PMF barrier (taken from the maximum of each plot) of ~6 kcal mol⁻¹ for the 0H state was significantly reduced to ~3-4 kcal·mol⁻¹ in the 1H and 3H states. This supports the experimental evidence, where in the absence of RMT binding to the M2 channel water permeation was observed in the mono- (1H) and triple- (3H), but not the non- (0H) protonated states [22,28].

Regarding the M2 channel with RMTs bound outside, the 0H and 1H states display relatively high PMF energies of ~4-9 kcal·mol⁻¹, suggesting that water permeation is energetically unfavourable in these systems. This is consistent with the zero water density plots of both the 0H and 1H states (Figure 4). In contrast, the energy barrier for water transport is significantly lower at ~3 kcal·mol⁻¹ in the 3H state of the M2-RMT_out complex, consistent with the non-zero water density observed throughout the channel. Interestingly, the PMFs for M2-RMT_in complexes in all the protonation states are significantly higher than ~7.5 kcal·mol⁻¹ indicating no water transport was processed where a single RMT blocking inside the M2 channel. These results clearly demonstrated that the RMT binding inside the M2 pore has potentially inhibited the M2 machinery rather than the RMTs allosteric binding outside.

In summary, the information here leads us to conclude that the RMT-mediated inhibition via their binding to the outside of the pore in the closed and intermediate conformations at high pH conditions can prevent proton conductance by interrupting the formation of the water wire along the channel. In contrast, the transportation cannot be inhibited at a low pH condition, i.e. in the 3H-open conformation.

Conclusion

Based on the MD simulated results presented here, the inhibition efficiency of the four RMTs bound outside the M2 channel pore was considerably lower than that observed when bound inside the channel. The key parameter determining the drug outside binding efficiency is based on the hydrogen bond network centred on the Asp44 residue which interacts with RMT to further stabilise the gating channel by inter-subunit hydrogen bonding with the Trp41 gating residue. This network was eliminated when the positive charge was increased on the selective His37 residue, due to the decreased van der Waals contact of the Trp41 indole ring being affected by the strong electrostatic repulsion. The effect was manifested at pH ≤ 6 (triple protonation state on the His37 tetrad) where the RMT-Asp44-Trp41 hydrogen bond network was completely lost, leading to a rearrangement of the indole ring of the Trp41 gate to flip to an open conformation. The water permeation throughout the M2 channel that was consequently observed was ranked as 0H < 1H < 3H which is inversely consistent with the observed PMF. Taken together, the MD-based simulation results have clearly explained the role of RMT at the molecular level in the allosteric binding site of the M2 channel at high pH conditions but not at a low pH.

Acknowledgements

P.I. thanks the Post-Doctoral programme from the Commission on Higher Education. S.H. thanks the Centre of Excellence for Petroleum, Petrochemicals, and Advanced Materials, Chulalongkorn University. Research facilities, software packages and computing time were provided by the Computer Centre for Advanced Research and the Computational Chemistry Unit Cell, Faculty of Science, Chulalongkorn University. Molecular graphics images were produced using the UCSF Chimera package from the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco (supported by NIH P41 RR-01081).

Declaration of interest

This work was supported by The Thailand Research Fund.