Docking and Molecular Dynamics Simulations in Potential Drugs Discovery: An Application to Influenza Virus M2 Protein

Biological Faculty, Lomonosov Moscow State University (MSU), Leninskie gory 1, Moscow 119234, Russia Department of Physics, University of Osnabrück, Barbarastr 7, Osnabrück 49076, Germany Department of Health Sciences (DISSAL), University of Genoa, Via Antonio Pastore 1, Genoa 16132, Italy Laboratories of Biophysics and Nanotechnologies (LBN), Department of Experimental Medicine (DIMES), University of Genoa, Via Antonio Pastore 3, Genoa 16132, Italy Nanoworld Institute Fondazione EL.B.A.Nicolini (NWI-FEN), Largo Redaelli 7, Pradalunga, Bergamo 24100, Italy


Introduction
With 10-20% of the worldwide population catching Influenza-Like Illness (ILI) every year (Peasah et al., 2013), influenza is an important issue in Public Health, having a huge impact on healthcare systems and on society, both in terms of disease burden and costs (Gasparini et al., 2000;2002;2003;Lai et al., 2011;Gasparini et al., 2012;. Recently, Peasash and collaborators carried out a systematic review of the literature and found that the per capita cost of a case of influenza illness ranges from $30 to $64.22 (Peasah et al., 2013).
Altough influenza vaccines are an effective weapon against influenza and are cost-saving (Peasah et al., 2013), antiviral drugs could offer an opportunity to alleviate the burden of influenza, both for treating influenza symptoms and for post-exposure prophylaxis (Jackson et al., 2011;Tappenden et al., 2009).
In particular, Neuraminidase-Inhibitors (NI) such as oseltamivir and zanamivir proved to be cost-effective and clinically effective, whilst evidence for adamantanes M2-inhibitors such as amantadine and rimantadine is scarser (Jackson et al., 2011;Burch et al., 2009).
However, clinical resistance against anti-viral drugs is emerging and for this reason there is an urgent need to widen and diversify the armamentarium of antivirals (Ison, 2011).
Considering that in the field of medicinal chemistry the discovery and development of a completely New Molecular Entity (NME) or compound is particularly expensive, in term of time and costs, this could be carried out using two different approaches: To design/optimize new derivatives from an existing lead (such as the second-generation NI laninamivir and peramivir) and to repurpose/reposition old and already existing drugs for new potential clinical applications (Bastos and Coelho, 2014;Ekins and Williams, 2011).
The latter approach, also termed as drug retasking or reprofiling, has already given promising results. While drug retargeting was initially serendipitous, it has been later more systematically developed and exploited, also combining advanced biochemical, biophysical and bioinformatics/cheminformatics techniques.
Since Omeprazole Family Compounds (OFC) block the "proton pump" (hydrogen-potassium ATPase), we hypothesized that they could interfere with the mechanism of fusion of the virus envelope and endosomal membrane, thereby hindering the M2 proton pump mechanism of influenza viruses as well.
Recently, to test this hypothesis, we carried out a matched case-control study which showed that subjects treated with OFC displayed a lower risk of catching Influenza-Like Illness (ILI) (adjusted odd ratio = 0.29, 95% CI: 0.15-0.52), whilst this risk was about six times higher in unvaccinated non-OFC users (Gasparini et al., 2014).
This epidemiological finding seems to suggest that OFC could exert a protective effect against ILI, putatively blocking M2 ion channel, which is a homotetrameric type III integral membrane tetrameric pH-gated protein. It contains a small N-terminal ectodomain, a single transmembrane domain and Cterminal cytoplasmic tail (Pinto and Lamb, 1995). The transmembrane domain acts as both a signal sequence and a membrane anchor during protein synthesis. The HxxxW motif of the inner membrane spanning residues proves to be critical to the ion channel activity, which is highly selective for protons, depending on a histidine residue in the transmembrane domain.
However, the exact mechanism for transport of protons with high selectivity is not known.
Briefly speaking, the activity of the ion channel in the viral lipid envelope is essential for the life cycle of the virus. The low pH of an endosome activates the M2 channel prior to hemagglutinin-mediated fusion. Conductance of protons acidifies the viral interior and thereby facilitates dissociation of the matrix protein from the viral nucleoproteins-a required process for unpacking of the viral genome. In addition to its role in release of viral nucleoproteins, M2 in the Trans-Golgi Network (TGN) membrane prevents premature conformational rearrangement of newly synthesized hemagglutinin during transport to the cell surface by equilibrating the pH of the TGN with that of the host cell cytoplasm. H37 and W41 act putatively as a primary gate, while V27 acts as a secondary gate.
The antiviral drug amantadine inhibits the replication of the virus by putatively binding to the Transmembrane Domain (TMD) of the M2 proton channel. Other scholars suggest a potential role also of Ser 31 and Ala 30. Rimantadine binds to the pocket of Asp 44 and Arg 45 (forming also hydrophobic interactions with residues 40, 42, 43) (Gu et al., 2013;Kolocouris et al., 2008;Kozakov et al., 2010;Intharathep et al., 2008;Schnell and Chou, 2008).
Summarizing, two main different sites for drug interactions have been proposed. One is a lipid-facing pocket between two adjacent TM helices (around Asp-44), at which the drug binds and inhibits proton conductance allosterically. The other is inside the pore (around Ser-31), at which the drug directly blocks proton passage.
All three ligands were docked to experimentally known binding sites. Then, we used MD simulations to equilibrate the obtained structures in a model lipid membrane. Binding free energies to the M2 protein tetramer in water environment (with backbone atoms constrained to the original X-ray structure) and to the M2 tetramer pre-equilibrated in membrane were calculated. We calculated with the usage of the thermodynamic integration method the free energies of binding for amantadine, rimantadine, omeprazole compounds ( Fig. 1): A correlation with experimental data was shown. We demonstrate that the thermodynamic integration method predicts free energies of ligand binding better than molecular docking while embedding of M2 protein in a membrane further improves the calculated free energy values. Free energy calculations imply omeprazole as a potential anti-viral drug.

M2 Protein Preparation
M2 protein was retrieved and downloaded from Protein Data-Bank (PDB code 3C9J) in protonation state as it was get at pH 5,3. For systems 1-3 (Table 1) M2 protein was put in the cubic box with water and ions, for systems 4-6 M2 protein was placed inside the fully hydrated DPPC membrane.

Docking Procedure
The sites of binding for rimantadine and amantadine are well known from crystallographic data (Schnell and Chou, 2008;Cady et al., 2010), for omeprazole we chose the same site of binding as for amantadine. Docking procedure was provided with the usage of DOCK6 program (Brozell et al., 2012).

Ligand Preparation
Charges for all three ligands were calculated with GAMESS program (Alexeev et al., 2012) using Mulliken charges and Hartree-Fock method in 6-31G* basis with the pre-optimizing (with DFT method with B3LYP functional) molecule geometry. Protonation states were predicted at pH 5,3 as the protein.

MD Simulation Details
The MD simulation was carried out using the program Gromacs, version 4.5 (Pronk et al., 2013). The GROMOS force field was used to characterize the compounds (Van Gunsteren et al., 1998). The molecules were placed in the center of a dodecahedron box, which was subsequently filled with TIP3P water molecules. All ligands were neutral while M2 monomer was negatively charged. About 4 additional 4 Na + atoms were added in each box. After standard equilibration procedure a MD cycle was run. All simulations are presented in Table 1.
We used standard isothermal-isobaric conditions (NPT ensemble): Isotropic pressure coupling (Parrinello-Rahman barostat, time constant 2 ps) and constant temperature at 298 K (Nose-Hoover thermostat time constant 2 ps). Protein with ligands and water were coupled independently to the heat bath. Periodic boundary conditions were applied in all three dimensions. All bond lengths were kept constant using the Linear Constraint Solver (LINCS) algorithm. Time step was 2 fs. Long-range electrostatic interactions were treated with the Particle-Mesh Ewald (PME) algorithm (real space cutoff 1 nm, FFT grid spacing 0.18 nm). The Lennard-Jones potentials were computed by using a cutoff length of 1.2 nm.

Binding Constant Calculation: Thermodynamic Integration
Thermodynamic integration (De Ruiter and Oostenbrink, 2011) is a method used to compare the difference in free energy between two given states X and Y whose potential energies have different dependencies on the spatial coordinates. This method can be used to calculate binding constants between a ligand and a substrate.
In the thermodynamic integration method, the free energy A can be written as a function of the coupling parameter λ Equation (1): where, k B is the Boltzmann's constant, T is the temperature, Q the partition function of the canonical ensemble (canonical partition function).  By differentiating this expression and considering the free energy between two states X and Y we obtain Equation (2): It should be noted that, formally, the energy E in Equation (2) should be replaced by the Hamiltonian H of the system. H can be considered as the total energy of the system, that is the sum of the kinetic and potential energies. However, if only the potential energy of the system changes as a function of λ, the free energy difference will depend in practice only on the potential energy.
Thermodynamic integration method allows calculation of the binding energy ∆G b . For this purpose the thermodynamic cycle is applied and the binding free energy, G b , is given by Equation (3): ∆G 1 and ∆G 2 are the free energy changes associated with the decoupling of the ligand (and counter-ions) from the solvent (no protein present) and with the decoupling of the ligand (and counter-ions) from the solvent and the protein, respectively.
Decoupling in this context means that in the final state the ligand and eventual counter-ions do not "see" the rest of the system and are therefore effectively in a gas phase state. This procedure is referred to as the double decoupling method. The free energy difference, ∆G, was estimated using the thermodynamic integration formula in Equation (3). The integration was carried out numerically by means of the trapezoidal method.
To describe the decoupling of the system, 9 separate λ points were simulated (0, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 1). The Hamiltonian of the system can be written as the sum of the Hamiltonians of the ligand, H L (λ) and of the counter-ions, H I (λ): H(λ) = H L (λ)+H I (λ) and the free energy of binding becomes Equation (4):

Docking Results
Adamantanes and omeprazole are docked inside the channel-in hydrophobic "pockets" of the protein (Fig. 4, grey representation of protein-ligand complex). The free energies calculated from docking score are presented in Table 2. These free energy values significantly overestimate experimentally known binding constants.

MD Simulation Results
Equilibrium MD simulation in water shows that M2 protein changed its structure dramatically (Fig. 2). The Root-Mean-Square-Deviation (RMSD) parameter is increasing up to 0,46 nm and number of H-bonds is decreasing (at least 7 hydrogen bonds are cut). The protein in water becomes less stable and the secondary structure is disrupted.
Equilibrium MD simulation of the M2 protein in membrane demonstrates more stable structure (Fig. 3). According to our calculations inside the membrane crystal structure becomes more compact what is in a good agreement with NMR data (PDB code 2RLF).
According to the method of thermodynamic integration, binding energies are calculated (Table 2). Since water environment significantly affects channel structure as we have shown before, we fixed Ca atoms of the protein in order to prevent its disrupter during free energy calculations procedure (Fig. 4, blue representations of protein-ligand complexes). The results we obtained from the TI simulations in aqueous environment are in better agreement with experimental constants of binding (Table 2) comparing to those obtained from molecular docking but still differ from the experimental values in several magnitudes of the binding constant.
We have chosen fully hydrated dipalmitoylphosphatidylcholine (DPPC) membrane as a model membrane in our research ( Fig. 3a and b). TI into X-ray TI into DPPC Docking into TI into X-ray TI into DPPC Docking into TI into X-ray TI into DPPC Complex X-Ray structure structure in water membrane X-Ray structure structure in water membrane X-Ray structure structure in water membrane ∆G, kJ/mol -12. It was prequilibrated during the 100 ns. After this, M2 protein was inserted inside the bilayer with the usage of standard GROMACS utilities. The whole system was equilibrated again for 5 ns and then productive MD simulations were run. Analogously to the simulations in the water environment we calculated binding free energies for the membrane embedded M2 tetramer using TI approach. Calculated binding constants (Table 2) appear in much better agreement with the experimentally estimated constants. The number of hydrogen bonds represented in Table 2 demonstrates weak dependence on the ligand structure and should not be taken in account for the explaining of the mechanism of binding. The M2 tetramer inserted into lipid bilayer adopts more compact conformation as we have shown before. This conformational change leads to reshaping of the channel interior, which is known as the binding site for the M2 inhibitors studied in the present study. This "improvement" of the original X-ray structure allows more realistic positioning of the ligands in the interior of the channel. Also it is important to notice that ligands positioning are all quite different on the Figure 4. Especially for such a big ligand as omeprazole. Docking procedure places the ligand inside the pocket and docking free energy provides score of all interactions which are described inside the algorythm. Molecular dynamics simulations allow side chains and ligand to be more flexible. Thus, molecular dynamics simulations of the transmembrane region of M2 protein tetramer in its natural lipid environment seem to be essential for correct prediction of ligand affinities.
Drugs belonging to OFC group such as omeprazole, lansoprazole and pantoprazole selectively and irreversibly inhibit the part of the "proton pump" that performs the final step in the acid secretory process. In 2005, Sasaki and collaborators demonstrated an anti-Rhinovirus activity of lansoprazole, which was probably due to an endosomal anti-acidic mechanism (Sasaki et al., 2005). However, the link between OFC and ILI is not very clear. Some scholars speculate from laboratory and clinical evidence that the increase in gastric pH caused by OFC may be linked to increased bacterial colonization of the stomach and may predispose patients to an increased risk for respiratory infections (Sultan et al., 2008). However, Sultan and collaborators (Sultan et al., 2008) conducted a systematic review and meta-analysis of the Randomized Controlled clinical Trials (RCT) and found that there is no evidence of this putative link (OR 1.42, 95% CI 0.86 to 2.35; p-value = 0.17).
Our study seems to suggest a potential role of OFC in the armamentarium of drugs against influenza.

Conclusion
MD simulations can provide insights of ligandprotein interactions. We show that free energy of binding should be calculated for a protein in its natural environment especially for membrane proteins. Calculated binding constants and structural features shed light on atomistic properties and functionality of the proteins. At the same time they appear helpful for the drug repositioning: • All three ligands bind inside the channel pore and are located in hydrophobic "pockets" of the ion channel • MD simulations in membrane are essential in order to fix the X-ray structure of the M2 protein channel and predict binding free energies in the most precise way • Ligands form just few or no H-bonds with the M2 tetramer, the main contribution in the binding free energy is given by the hydrophobic interactions • TI method allows to predict binding constants for Rimantadine and Amantadine very close to the experimental values • Calculated free energy of binding for Omeprazole is the lowest among the three studied compounds, what, along with the validation of the method against the known binding constants for Rimantadine and Amantadine, implies high putative anti-viral activity of this compound The fact that the omeprazole binds two sites in the molecule of the M2 protein could make the mechanism of action of this compound less sensitive to variations of the M2 protein, which as is known, have led to the emergence of many strains of influenza virus resistant against amantadanes (Sheu et al., 2011;Pielak and Chou, 2010).