Research Article

Studying the Effect of Amino Acid Substitutions in the M2 Ion Channel of the Influenza Virus on the Antiviral Activity of the Aminoadamantane Derivative In Vitro and In Silico

Introduction

The M2 proteins of Influenza A virus are transmembrane proteins that form proton-selective channels in the lipid membrane. M2-channel is a homotetramer structure assembled from four subunits of the M2 protein. Polypeptide chains with a length of 97 amino acid residues are partially helical in the region of the transmembrane domain. ¹ Proton conductivity of M2 channels is necessary for virus replication. Proton-selective M2 channels regulate pH inside the viral particle during virion penetration into the endosome of the cell by acidification the inner space of the viral particle, which results in dissociation of the matrix protein M1 complex with ribonucleotide and leads to the release of genetic material into the cytoplasm of the host cell. ^2,3 Another important role of M2 protein appears at a later stage during viral assembly, when newly synthesized viral proteins are transported to the cell membrane surface. At low pH values, the virus proteins, in particular, hemagglutinin (HA), can undergo premature and undesirable transformation; in this case the M2 proton channels pump hydrogen ions from the Golgi complex to maintain a sufficiently high pH value and preserve the transport form of the newly formed HA. ^4,5

The M2 viral protein channel is a target for antiviral therapeutic drugs, namely amantadine and rimantadine, which belong to the class of adamantane compounds. Aminoadamantanes have been used for the treatment and prevention of influenza A since the 1970s. Their economic and synthetic availability made them ideal for treating seasonal influenza epidemics around the world. However, because of the widespread use of these drugs, the influenza A virus, which has a high mutation rate, underwent genetic modification making it resistant to the action of adamantane-type drugs. In laboratory strains, scientists observed various mutations in the transmembrane domain of the M2 channel; however, the wild type is characterized by three main mutations, the amino acid substitutions of the serine residue with asparagine in 31 positions of the M2 protein being considered critical for the ability of aminoadamantanes to block the functions of the channel.

The Centers for Disease Control and Prevention (CDC), in association with WHO, recommend no use of rimantadine and amantadine () for the treatment of influenza because their resistance has exceeded 90%.

Close

Figure 1. Structure of aminoadamantans - amantadine (a), rimantadine(b) and amino acid derivative of rimantadine L-histidil-1-adamantayl ethylamine (HCl*H-His-Rim) (c).

The aim of this work is to study the mechanism of antiviral effect of HCl*H-His-Rim compound () on the proton-conducting channel M2 of influenza A virus through the identification of amino acid substitutions in the pore of M2 channel, affecting the antiviral activity of the compound and analyze the most significant mutations by molecular docking methods.

Materials and Methods

Compound. HCl*H-His-Rim

The formation of a peptide bond between the carbocycle of rimantadine and di-tert-butyloxycarbonyl-L-histidine was carried out in a single stage under the conditions of mixed anhydride method in the equimolar ratio according to the procedure described. ^6-8 The removal of protective tert-butyloxycarbonyl groups was carried out in a solution of 4 N HCl in ethyl acetate followed by vacuum drying to a dry foam. Yield: 92%. M.p. 210°С (decomposition); [α] ²⁰ _D = +6° (c = 1; C₂H₅OH), IR: υ(NH) 3248 cm^-1; υ(NH^im) 3139 cm^-1; υ(C=O) 1672 cm^-1. m/z found [M + H]⁺: 317,746; [M + Na]⁺: 339,761; calculated M (C₁₈H₂₈N₄O) 316,441 (Figure S1). ¹H NMR (D₂O, ppm): 8.81 (s, 1H, NH-CH=NH), 7.52 (s, 1H, NH-CH=C), 4.79 (H₂O (HDO) + NH-atoms), 4.31 (t, 1H, CH2-CH-NH₂, 8 Hz), 3.49 (q, 1H, NH-CH-Ad, 7Hz), 3.40 (d, 2H, His-CH2-CH, 8 Hz), 1.90 (s, 3H, Adamantyl β-H), 1.60 (dd, 6H, Adamantyl γ-H, 56 Hz, 11 Hz), 1.22 (dd, 6H, Adamantyl α-H, 41 Hz, 11 Hz), 1.04 (d, 3H, CH₃-CH-Ad, 7 Hz);¹³C NMR (D₂O, ppm): 166.9 (C=O), 134.1 (NH-CH=NH), 126.4 (NH-C=CH), 118.3 (NH-C=CH), 54.7 (CH3-C-Ad), 52.5 (CH₂-CH-NH₂), 37.7 (Adamantyl α-C), 36.3 (Adamantyl γ-C), 35.1 (CH₂-CH-NH₂), 27.8 (Adamantyl β-C), 26.1 (Adamantyl C-CH-NH₂), 13.1 (CH₃-C-Ad) (Figure S2).

Molecular docking, dynamics and free energy calculations

The quantum-mechanical model of the ligand was generated in HyperChem 8.0.8 software product by Hypercube (https://hyper.com/).

The structure of the transmembrane domain M2 S31N from RCSB Protein Data Bank (structure code 2KIH) was taken as a target for docking to be performed. The docking of the ligand model into the protein channel was carried out using an ant colony optimization algorithm in PLANTS program. ⁹ For the ligand, flexible docking with the possibility of cis-to-trans transition via an amide bond was used. For the initial assessment and ranking of docking results, the ChemPLP scoring function was used. ¹⁰ To study the binding stability of HCl*H-His-Rim in conformations selected from molecular docking solutions with M2 S31N and M2 S31N_A30T, molecular dynamic (MD) modeling of this complex was performed in GROMACS software package (version 2018.8). ¹¹ Ligand models with arrangement of partial electric charges of atoms, parameters of Van der Waals interactions, stiffness constants for valence bonds, valence and dihedral corners were processed using the CGenFF WEB service. ¹² The layout and parameterization of MD-systems of protein-ligand complexes in lipid bilayer was carried out in the CHARMM-GUI: Membrane Builder WEB service ¹³ using a set of CHARMM ¹⁴ full atomic force fields of version C36m. ¹⁵ The model lipid bilayer consisted of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine molecules. The model was built in a hexagonal box with dimensions X = Y = 100 Å. The parameter Z was selected automatically based on the condition that the distance from lipid atoms to the upper and lower box should be at least 22.5 Å. These volumes were filled with water molecules in accordance with the transferable intermolecular potential with 3 point (TIP3P) model and with the addition of neutralizing Na⁺ and Cl- ions of 0.1M concentration. The same concentration of counter ions was introduced into the ligand-solvent system.

The simulation was carried out by the molecular dynamics method. To integrate Newton’s equations of motion with a step of 2 fs, a leap-frog algorithm was used. The path length was 20 ns. The long-range electrostatic interactions were calculated by the particle-mesh-Ewald method (PME) using the cubic interpolation function. To describe the van der Waals interactions, we used the force-switch smoothing function in the range from 10 to 12 Å. The pressure in the system was controlled by the Parrinello-Rahman pressure coupling ¹⁶ at the level of 1 bar. The temperature of the MD-system was kept constant by the velocity rescale temperature coupling ¹⁷ at 310 K. Molecular dynamics procedure was preceded optimizing atomic parameters by the method of steepest descent and simulation of the system with a position restraints of heavy protein and ligand atoms in canonical (NVT) ensemble and isothermal–isobaric (NPT) ensemble (the Berendsen’s pressure coupling ¹⁸ was used) according to Table S1.

The density map of the radial distribution of water averaged over the path of the MD simulation was constructed using the g_ri3Dc program. ¹⁹

The relative affinity of ligand binding was determined by calculating the linear interaction energy ^20,21 (LIE) using the gmx lie program of GROMACS software package. To make this, the average energies of van der Waals and electrostatic ligand-water interactions in the ligand-water system were previously calculated during MD simulation of 20 ns. Protein-ligand MD-simulation trajectories were re-processed using shift potential at 12 Å cutoff to calculate the potential energy of Coulomb interaction between the ligand and the environment. The empirical constants α and β, which determine the contribution of van der Waals and electrostatic interactions to the estimation of the binding energy, were selected in accordance with Hansson et al ²¹ and were equal for HCl*H-His-Rim: α = 0.18 and β = 0.5. The shift coefficient γ was taken equal to zero due to the impossibility to select it at this stage of the study.

The ligand-protein free energy of binding was calculated usingthe thermodynamic cycle described Aldeghi et al. ²² The first stage consisted of collecting the derivative values of the system Hamiltonian by the Kirkwood coupling parameter ²³ (λ-parameter) (for Thermodynamic Integration (TI)) and the difference in the potential energy of the system between the current and other states (for Free Energy Perturbation (FEP)-based methods). For protein-ligand and protein-water systems, MD simulations were performed in the NPT ensemble for 2.4 ns using a leap-frog integration algorithm of stochastic dynamics equations for each set of λ-parameters. In order to retain the ligand in the binding region during its decoupling, additional intermolecular restrains were introduced as described Aldeghi et al. ²² For the protein-ligand system, λ changed according to the sequence of λrestr (0.0, 0.01, 0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0) and λ_coul (0, 0.2, 0.4, 0.6, 0.8; 1); λ_VdW (0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0). In the case of the ligand-water system, coupling parameters changed only for the Coulomb and van der Waals ligand-water interactions: λ_coul (0, 0.25, 0.5, 0.75, 1); λ_VdW (0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0). The correction for the energy of the restrains introduced for the ligand-water (ΔG_{solv_rest}) was calculated analytically using the formula reported Boresh et al. ²⁴ Prior to the main stage of MD simulations, in addition to the general equilibration of the system described above, for each value of the coupling parameter λ, energy minimization and equilibration were carried out with restrictions imposed on the positions of heavy protein and ligand atoms in the NVT and NPT ensemble for 100 ps. Van der Waals interactions were calculated using the Lennard-Jones (LJ) “soft core” potential. The PME method was used to calculation both long-range Coulomb and van der Waals interactions.

Data of FEP MD simulations were processed in a program written in Jupyter Notebook using alchemlib and pymbar libraries. The first step was the parsing of xvg files with the extraction of reduced potentials using alchemlyb. ²⁵ The second step was to select uncorrelated values of reduced potentials (U_ij) with cutting off the non-equilibrated region of trajectories using the detect equilibration function of the pymbar library. ²⁶ For each window, the difference of reduced potentials between current and the neighboring window was used to calculate statistical inefficiency. In the third step, the statistically independent values of U_ij were processed using the Multistate Bennett Acceptance Ratio (MBAR) method. ²⁷ The correctness of the calculations was controlled by the overlap matrix between the λ-states, as well as by comparing the obtained ΔG values with the values obtained by other methods: Exponential averaging, ²⁸ Bennett Acceptance Ratio, ²⁹ TI. ³⁰

Since the total charge of ligand atoms in the MD simulation was +1 e, it was necessary to calculate the corrections for Ewald summation and periodic boundary conditions. They were calculated using analytical scheme as described Rocklin et al. ³¹ The protein charge was completely compensated by counterions and effective Q_P = 0. The box was converted to orthorhombic to simplify these calculations. The average linear box parameters t were 9.682 × 8.384 × 10.395 nm for the protein-ligand system, and 3.433 × 3.237 × 2.804 nm for the ligand-solvent system. The difference in linear parameters between the different complexes was less than 5%. Roughly the same fluctuations in the parameters were observed as a result of the pressure coupling algorithm work in the NPT simulation. 1/(4πε_o) was assumed to be 138.93545585 kJ nm e^-2 mol^-1, ε_s for TIP3P water was 82, ρ_310K = 993 kg/m³, and γ_S = 0.0764 e nm². Residual integrated potential (RIP) correction was calculated using APBS-mem ³² and APBS. ³³ In calculating RIP and discrete solvent effects (DSF) corrections, the following parameters of the model lipid bilayer were used: D_M = 4.51 nm, D_HEAD = 0.9 nm. ³⁴ The relative permittivity of protein and solvent were ε_s = 2, ε_H = 80.

Plasmid transfection

The source virus and mutants based on it were obtained by reverse genetics. A system of eight plasmids containing DNA copies of genome segments of the influenza virus strain A/PuertoRico/8/34 (H1N1) was used in this work. ³⁵ Nucleotide substitutions for mutant production were introduced using the commercial QuikChange XL Site-Directed Mutagenesis Kit (“Stratagene”, USA). Plasmid transfection was performed in mixed culture HEK-293T-MDCK. ³⁶

Viruses obtained by plasmid transfection were passivated to the allantois cavity in 10-day-old chicken embryos. Infected chicken embryos were incubated for 48 hours at 37°C and then cooled overnight at 4°C. Virus-containing allantois fluid was collected under sterile conditions and the virus content was evaluated as a titer obtained in the hemagglutination reaction ³⁷ (HGA) expressed in terms of the number of hemagglutinating units (GAUs). The preparations were stored at -80°C.

Polymerase chain reaction and sequencing

Viral RNA was isolated from the virus-containing allantois liquid using the RNeasy Mini Kit (Qiagen). Reverse transcription and PCR were performed with primers to the influenza A virus M gene. Amplification products were purified using the QIAquick PCR Purification Kit (Qiagen). Sequencing was performed using the DNA ABI Prism 3130 automatic DNA sequencer (Applied Biosystems) and the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems). Nucleotide sequences were analyzed using the DNASTAR Sequence Analysis Software Package (DNASTAR Inc.).

Determination of antiviral activity

The antiviral activity of the test compound (HCl*H-His-Rim) against influenza A/PuertoRico/8/34(H1N1) was determined using immunoassay (EIA) in MDCK cell culture, during the dilution of the virus 10^-4-10^-6 breeding. The compound was added in concentrations of 10.0 µg/mL (31.5 μM) at the same time as the cell monolayer was infected. The test compound and dilution of the virus were prepared in a medium supplemented with TPCK trypsin (2.5 µg/mL) in tests using influenza virus A/PuertoRico/8/34(H1N1) strains and mutants of this virus. The plates were incubated for 24 h at 37°C, after which the reaction was stopped by treating the cells with 80% acetone in a phosphate buffer. Cellular EIA was performed according to the previously described method. ^38,39 The percentage of viral activity inhibition of the compound was determined as the ratio of [experimental optical density of OD₄₅₀ – OD₄₅₀ cell control / OD₄₅₀ virus control – OD₄₅₀ cell control] × 100%.

Results and Discussion

The development of new therapeutic agents is an urgent task due to the fact that classical direct-acting drugs such as rimantadine and amantadine are not effective because of drug resistance formed. The function of the ion channels can be blocked by small inhibitor molecules (in particular, adamantane-series preparations), which leads to a significant decline in the reproduction of viral particles. The action mechanism of the aminoadamantane derivative is the subject of this study.

The mechanism of inhibition of amantadine and rimantadine has been intensively studied during the last two decades. Before the appearance of high-resolution crystallographic structures, the drug binding site was predicted from the position of mutations that led to the emergence of resistant strains of influenza A virus. Mutations that are responsible for resistance to amantadine or rimantadine are found in the following protein M2 positions: 26, 27, 30, 31, 34 and 38. ^40-42 The authors ⁴² generalized all hypotheses known to date regarding the binding sites of amantadine with the M2 protein. Most of these positions, including 27, 30, 31 and 34, are located on the inner surface of the M2 channel, which led to the hypothesis of drug penetration into the channel.

Studies ⁴³ showed that the formation of a classical hydrogen bond is not a prerequisite for the inhibitory effect on proton conductivity. There is an evidence of a mechanism in which the amantadine amino group does not interact with any of the M2 channel pore polar groups. For example, the S31A mutant was sensitive to rimantadine, ⁴⁴ which indicates that S31 does not play a decisive role in binding to the drug.

The inner surface of the M2 canal is composed mainly of non-polar amino acid residues, and the formation of reliable hydrophobic interaction with the ligand is sufficient for antiviral activity to occur. Duff et al ⁴⁵ showed that amantadine amino group can interact with the imidazole imidazole pairing of His37 residues; therefore the presence of an additional imidazole ring in the channel pore may lead to disruptions in the protonation process due to competitive reactions with 2HCl*H-His-Rim imidazole ring.

Thus, the main supposed mechanism of antiviral action of synthetic compound 2HCl*His-Rim is probably similar to that of amantadine on the M2 ion channel. ⁴¹

Based on the predictions of primary molecular docking by reverse genetics methods, ⁴⁶ mutants of influenza virus A/PuertoRico/8/34(H1N1) containing amino acid substitutions in position V27, I28 and A30 were created.

The ability of the ligand to inhibit mutant virus replication in comparison with native A/PuertoRico/8/34(H1N1) influenza virus has been studied by computer simulation of ligand (2HCl*His-Rim) interaction with target protein (M2 channel) in the presence of amino acid substitutes V27A and A30T. Moreover, the mutation of A30T seemed most likely to achieve a decrease in the antiviral properties of the ligand compound due to the maximum difficulties in binding to the target protein. This is due to the large steric restrictions because of the increase in the volume of the side group. The V27A mutant binds optimally closer to the N-terminus of the channel, but in general the total binding energy was the same as that of the reference structure.

The situation with the double replacement of V27A + I28V was less favorable to achieve a result than A30T. Replacing the aliphatic group 2-butyl with the aliphatic group isopropyl could not significantly affect the stereometry inside the channel at position 28. The residue I28 was not turned into the channel, therefore it had little effect on the steric difficulties for the ligand, but could lead to a decrease in the strength of the hydrophobic interaction of the binding site with the adamantane core of the inhibitor molecule.

A reverse genetics system was used to prove the effects. This system is successfully used for scientific and applied research. ⁴⁵

The obtained mutants of influenza virus A/PuertoRico/8/34(H1N1) were tested in vitro for sensitivity to HCl*H-His-Rim. Influenza virus A/PuertoRico/8/34(H1N1) was sensitive to HCl*H-His-Rim; the compound successfully inhibited the replication of this strain at a concentration of 10.0 μg/mL (31.5 μM), the percentage of inhibition in the breeding of 10^-5 virus was about 70% (Table 1). In turn, the three proposed mutants were less sensitive to the action of the compound than native A/PuertoRico/8/34(H1N1) influenza virus.

Table 1 shows the average of the twelve parallel experiments. The compound acquired the greatest loss of activity was mutant A30T. The mutants of influenza virus A/PuertoRico/8/34(H1N1) V27A and V27A + I28V were more sensitive to the compound, which showed comparable antiviral activity data.

For an in-depth study of the antiviral mechanism of the HCl*H-His-Rim compound, variants with native virus (reference) and mutant virus (A30T) were studied by molecular docking methods.

Molecular docking, energy minimization and equilibration results

Based on the value of the ChemPLP scoring function, 10 docking solutions were selected for the reference structure and for the mutant channel structure A30T (Table 2). The solutions form distinct clusters according to the position of the ligand in the channel, taking into account the fact that the channel structure is symmetrical in the direction perpendicular to the bilayer plane (4th order axis). From each cluster of solutions, one solution was selected that was the most identical in position and conformation in the reference and A30T. Further the structures of the complexes were placed in bilayer, solvated and procedures of energy minimization and equilibrating were performed.

The first variant of possible ligand binding by a channel (surface) is characterized by the surface arrangement of the ligand (). The adamantane skeleton is slightly submerged in the channel and is at V27 level. The histidine part of the ligand is adjacent to the side chain D24 and interacts with it through hydrogen bonds and electrostatic interactions. This part of the ligand also forms hydrogen bonds with water surrounding the entrance to the channel. It should be noted that after equilibration, the water molecules were found below the ligand H37.

Close

Figure 2. Structure of the outer (upper) part of the proton pump in combination with histidine rimantadine. Ligand conformations were detemined by molecular docking followed by the incorporation of complexes in the lipid bilayer, energy minimization and equilibration in NVT and NPT ensembles.

The second variant (middle) is characterized by the adamantane core being between 27 and 30 amino acid residues, and the amino acid part interacts with the side group N31 and the main chain of 30 amino acid residues. At this stage of research, water molecules were not detected in the channel between V27 and H37.

The third option is the deep location of the ligand in the channel (deep). The imidazole ring of the histidine part of the ligand is located almost across the channel and is at the level of 30-31 aa. At the same time, the adamantane part of the ligand is located 0.35-0.4 nm higher than the plane of the H37 imidazole rings, which are essential for the functioning of the proton pump. ³ Water molecules are located only above 30 amino acid residues.

Free MD-simulation. The passage of water in the channel

shows the standard deviation plots of the non-hydrocarbon ligand atomic coordinates gen with the respect to the initial position as a function of the MD simulation time. As can be seen from the plot, the most stable conformation in both cases is deep. The ligand in the surface binding position is the most labile. In this conformation, there is a strong oscillation of the amino acid part of the ligand in the solvent, which surrounds the entrance to the channel, while the position of the adamantane component atoms exhibits less changes. In a medium binding situation, the position of the ligand atoms in the reference protein is more stable compared to the A30T mutant.

Close

Figure 3. Standart deviation plots coordinates of the non-hydrogen ligand atoms with the respect to the initial position as a function of the MD simulation time simulation

In order to study the distribution of water in the pore channel, a map of the dependence of water density on the distance to the center of the channel and the coordinates of the perpendicular plane of the membrane was created (). In the case of surface fixation of HCl*H-His-Rim, both in the reference and in A30T, the density of water molecules in the channel is high both above and below H37. Apparently, the ligand in this position poorly blocks the water flow and, as a result, the activity of the proton pump. A significantly lower density of water molecules is observed when the ligand is in intermediate (mid) conformation. However, attention should be paid to the water-filled cavity between N31 and H37. At the previous stage of equilibration, no water molecules were observed in this volume. In case of long-time MD-simulation, water molecules penetrate from the inner side of the channel and pass through H37 and solvate the histidine part of the ligand, establishing a network of hydrogen bonds between the ligand and the channel wall. Significantly less diffusion of water into this region occurs from the upper part of the channel, since this diffusion is prevented by the adamantane hydrophobic fragment of the ligand. Thus, the water flow in the channel, in the case of the middle binding variant, is either insignificant or absent.

Close

Figure 4. Map of the radial density distribution of water molecules averaged by DM-trajectory.

In the case of low ligand binding, water does not penetrate into the channel below 30 amino acid residue. From below, water flow is also not observed throughout the entire MD simulation trajectory. Thus, the ligand being in this conformation completely blocks the water current in the channel of both the reference structure and A30T mutant.

Evaluation of the binding affinity of HCl*H-His-Rim conformations by the LIE Method. The role of water in binding the imidazole fragment of the ligand

To estimate the affinity by the LIE method, the terms of the interaction energy of the ligand with the environment in the ligand-water system were first calculated: < U_Coul> = -591.0 ± 3.0 kJ/mol, < U_lj> = -77.8 ± 1.2 kJ/mol. These values were used to calculate the dependence of ΔG_LIE on the MD-stimulation time (). The statistical parameters of the distribution of ΔG_LIE values are given in (Table S2). For surface binding, there is a larger scatter of ΔGLIE than for deep binding. It is explained by the higher ligand mobility in the surface conformation. In the middle conformation, a sharp change in the bind energy on the graph cannot be explained by the change of the ligand position in the binding site assumed. The short-range part of the Coulomb and Van der Waals energies of the interaction between ligand and protein and between ligand and solvent with counterions was analyzed with simultaneous viewing of the molecular plot. It was revealed () that in MD-simulation the energy of ligand-protein interaction undergoes no significant changes. At the same time, the interaction energy of the ligand with the solvent varies strongly and stepwise, both in the case of the average conformation of ligand-to-reference protein binding and with A30T. In the visual analysis of MD-simulation, it was found that these changes are associated with the penetration of water molecules into the cavity between N31 and H37. Water molecules pass mainly from the inner part of the membrane and then 3-4 water molecules constantly hydrate the histidine part of the ligand, establishing a network of hydrogen bonds between the ligand with water and amino acid residues of the protein. The structures of complexes with solvated ligand in the middle conformation were also chosen to calculate the free binding energy.

Close

Figure 5. Dependence of the relative free binding energy of HCl*H-His-Rim on the time of MD-simulation

Close

Figure 6. Short-range interaction energy for Coulumb (black) and Lennard-Jones term (red) for ligand-protein complex in middle conformation.

A sharp jump in binding energy is due to the transition of water molecules within the cut-off radius of the ligand atoms. The jump in the Coulomb term of the ligand-environment interaction energy in the A30T structure at 16 ns is due to a random character of approaching of the Cl- ion to ligand within the cut-off radius and is fast process.

The LIE data should be considered skeptically in terms of comparing possible ligand binding sites with each other, due to the significant difference in the hydrophobic properties of the ligand environment affecting the real value of the shift parameter γ. ⁴⁷ It should be noted that the Coulomb shift potential takes into account only interactions between particles at a distance not exceeding the cutoff radius. In addition, unlike the reaction field method, the dielectric permittivity value of the medium is not used here. The reaction field method is impractical to use, since the system is heterogeneous in terms of the dielectric permittivity value.

Calculations of binding free energy

Table 3 shows the free energy of binding for ligand-solvent and ligand-receptor system obtained from MBAR method with electrostatics corrections and total ΔG_bind. Total ΔG_bind was calculated by the formula: ΔG_bind = ΔG^prot- ΔG^solv, where

ΔG^x = ΔG^coul + ΔG^vdw+ΔG^rest+ΔΔG_DSF + (ΔΔG_NET+ΔΔG_USV) + ΔΔG_RIP,

ΔG^coul is the free energy Coulomb interaction,

ΔG^vdw is the free energy of Van der Waals interaction, defined as Lenard-Jones potential,

ΔG^rest is the free energy coupling restraints (in the case of the ligand-solvent system, it is calculated analytically according to the formula described by Boresch S, et.al. ²⁴ and includes a correction for the concentration of 1 N)

ΔΔG_DSF is the associated finite size discrete solvent correction term ³¹ ;

ΔΔG_NET is the correction for periodicity-induced net-charge interactions;

ΔΔG_USV is the correction for periodicity-induced undersolvation;

ΔΔG_RIP is the correction for residual integrated potential effects.

An infinite-system discrete solvent correction term (DSI) is not included in Table 3 because its value (-73.8 kJ/mol) is the same for protein-ligand and solvent ligand systems and is reset to zero when calculating ΔG_bind.

In the surface conformation of the ligand, binding is not energetically favorable since in both cases (reference and mutant) the value of the free binding energy is slightly higher than zero. In this conformation, the affinity to HCl*H-His-Rim is the same for both the mutant channel and the reference structure, which is explained by the remoteness of the position of the point mutation A30T from the ligand binding site.

In the case of middle binding without taking into account the solvation of the amino acid part of the ligand, ΔG_bind was found to be significantly positive, which suggests that this case is not energetically profitable. At the same time, when the amino acid part of HCl*H-His-Rim is solvated, the free energy of ligand binding to both the reference and mutant proteins decreases, which qualitatively corresponds to the LIE calculation results. In this case, ΔG_bind also have close values for the two channel structures, but this conformation, even taking into account ligand solvation, appears to be less energetically favorable than the surface one. According to the free energy calculations, the deep conformation is the most favorable in energy. In this case, according to the data of free dynamics, the water flow in the channel is completely blocked. In addition, in case of deep conformation, the binding of the ligand to the reference is much more advantageous than with A30T, which is consistent with the result of the experiment in vitro. In our opinion, this difference may be due to an increase in the volume of the side chain 30T and steric limitations when binding HCl*H-His-Rim to the mutant channel A30T. In addition, the potential barrier which is required to overcome the ligand upon binding in deep conformation appears to increase in the mutant form.

Conclusion

According to calculations of molecular dynamics, the conformation and position of the ligand in which the imidazole ring of the histidine part of HCl*H-His-Rim is at the level of 30-31 amino acid residues and the adamantane part of the ligand is 0.35-0.4 nm above the plane of imidazole rings H37 is the most energetically favorable and the most stable. In this ligand conformation, the water flow in the channel is completely blocked and binding the ligand to the reference is much more advantageous than with A30T, which is consistent with the result of the in vitro experiment. In our opinion, this difference may be associated with an increase in the volume of the side radical 30T and steric limitations upon binding HCl*H-His-Rim to the mutant channel A30T. In addition, the potential barrier which is necessary to overcome the ligand upon binding in deep conformation is apparently increasing in the mutant form.

Ethical Issues

Not applicable. This article does not describe any studies with human or animal subjects.

Conflict of Interest

The authors declare that they have no conflict of interest.

Acknowledgements

The publication has been prepared with the support of the “RUDN university program 5-100”.

This work was supported by the Ministry of Science and Higher Education within the State Assignment Federal Scientific Research Center Crystallography and Photonics of the Russian Academy of Sciences in part of molecular docking and modeling. The calculations were performed by Hybrid high-performance computing cluster of Federal Research Center Computer Science and Control of Russian Academy of Sciences. Available at: http://hhpcc.frccsc.ru).

Supplementary Materials

Supplementary file 1 contains Tables S1-S2 and Figures S1-S2.

References

1. Nieto-Torres JL, Verdiá-Báguena C, Castaño-Rodriguez C, Aguilella VM, Enjuanes L. Relevance of viroporin ion channel activity on viral replication and pathogenesis. Viruses 2015; 7(7):3552-73. doi: 10.3390/v7072786 [Crossref]
2. Liang R, Swanson JMJ, Madsen JJ, Hong M, DeGrado WF, Voth GA. Acid activation mechanism of the influenza A M2 proton channel. Proc Natl Acad Sci U S A 2016; 113(45):E6955-E64. doi: 10.1073/pnas.1615471113 [Crossref]
3. Duong-Ly KC, Nanda V, Degrado WF, Howard KP. The conformation of the pore region of the M2 proton channel depends on lipid bilayer environment. Protein Sci 2005; 14(4):856-61. doi: 10.1110/ps.041185805 [Crossref]
4. Henkel JR, Gibson GA, Poland PA, Ellis MA, Hughey RP, Weisz OA. Influenza M2 proton channel activity selectively inhibits trans-Golgi network release of apical membrane and secreted proteins in polarized Madin-Darby canine kidney cells. J Cell Biol 2000; 148(3):495-504. doi: 10.1083/jcb.148.3.495 [Crossref]
5. Hay AJ, Wolstenholme AJ, Skehel JJ, Smith MH. The molecular basis of the specific anti-influenza action of amantadine. EMBO J 1985; 4(11):3021-4.
6. Shibnev VA, Garaev TM, Finogenova MP, Shevchenko ES, Burtseva EI. Some pathways to overcoming drug resistance of influenza A virus to adamantane derivatives. Pharm Chem J 2012; 46(1):1-5. doi: 10.1007/s11094-012-0723-2 [Crossref]
7. Shibnev VA, Garaev TM, Finogenova MP, Shevchenko ES, Burtseva EI. Derivatives of 1- (1-adamantyl) ethylamine and their antiviral activity. Patent 2011 RU 2461544 C1 (Russian).
8. Deryabin PG, Garaev TM, Finogenova MP, Odnovorov AI. [Assessment of the antiviral activity of 2HCl*H-His-Rim compound compared to the anti-influenza drug Arbidol for influenza caused by A/duck/Novosibirsk/56/05 (H5N1) (Influenza A virus, Alphainfluenzavirus, Orthomyxoviridae)]. Vopr Virusol 2019; 64(6):268-73. doi: 10.36233/0507-4088-2019-64-6-268-273 [Crossref]
9. Korb O, Stützle T, Exner TE. An ant colony optimization approach to flexible protein–ligand docking. Swarm Intell 2007; 1(2):115-34. doi: 10.1007/s11721-007-0006-9 [Crossref]
10. Korb O, Stützle T, Exner TE. Empirical scoring functions for advanced protein-ligand docking with PLANTS. J Chem Inf Model 2009; 49(1):84-96. doi: 10.1021/ci800298z [Crossref]
11. Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJ. GROMACS: fast, flexible, and free. J Comput Chem 2005; 26(16):1701-18. doi: 10.1002/jcc.20291 [Crossref]
12. Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem 2010; 31(4):671-90. doi: 10.1002/jcc.21367 [Crossref]
13. Lee J, Patel DS, Ståhle J, Park SJ, Kern NR, Kim S. CHARMM-GUI membrane builder for complex biological membrane simulations with glycolipids and lipoglycans. J Chem Theory Comput 2019; 15(1):775-86. doi: 10.1021/acs.jctc.8b01066 [Crossref]
14. MacKerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ. All-atom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B 1998; 102(18):3586-616. doi: 10.1021/jp973084f [Crossref]
15. Huang J, Rauscher S, Nawrocki G, Ran T, Feig M, de Groot BL. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat Methods 2017; 14(1):71-3. doi: 10.1038/nmeth.4067 [Crossref]
16. Parrinello M, Rahman A. Polymorphic transitions in single crystals: a new molecular dynamics method. J Appl Phys 1981; 52(12):7182-90. doi: 10.1063/1.328693 [Crossref]
17. Bussi G, Parrinello M. Stochastic thermostats: comparison of local and global schemes. Comput Phys Commun 2008; 179(1-3):26-9. doi: 10.1016/j.cpc.2008.01.006 [Crossref]
18. Berendsen HJ, Postma JP, van Gunsteren WF, DiNola AR, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys 1984; 81(8):3684-90. doi: 10.1063/1.448118 [Crossref]
19. Beckstein O, Sansom MS. Liquid-vapor oscillations of water in hydrophobic nanopores. Proc Natl Acad Sci U S A 2003; 100(12):7063-8. doi: 10.1073/pnas.1136844100 [Crossref]
20. Aqvist J, Marelius J. The linear interaction energy method for predicting ligand binding free energies. Comb Chem High Throughput Screen 2001; 4(8):613-26. doi: 10.2174/1386207013330661 [Crossref]
21. Hansson T, Marelius J, Aqvist J. Ligand binding affinity prediction by linear interaction energy methods. J Comput Aided Mol Des 1998; 12(1):27-35. doi: 10.1023/a:1007930623000 [Crossref]
22. Aldeghi M, Heifetz A, Bodkin MJ, Knapp S, Biggin PC. Accurate calculation of the absolute free energy of binding for drug molecules. Chem Sci 2016; 7(1):207-18. doi: 10.1039/c5sc02678d [Crossref]
23. Kirkwood JG. Statistical mechanics of fluid mixtures. J Chem Phys 1935; 3(5):300-13. doi: 10.1063/1.1749657 [Crossref]
24. Boresch S, Tettinger F, Leitgeb M, Karplus M. Absolute binding free energies: a quantitative approach for their calculation. J Phys Chem B 2003; 107(35):9535-51. doi: 10.1021/jp0217839 [Crossref]
25. Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015; 1-2:19-25. doi: 10.1016/j.softx.2015.06.001 [Crossref]
26. Chodera JD. A simple method for automated equilibration detection in molecular simulations. J Chem Theory Comput 2016; 12(4):1799-805. doi: 10.1021/acs.jctc.5b00784 [Crossref]
27. Shirts MR, Chodera JD. Statistically optimal analysis of samples from multiple equilibrium states. J Chem Phys 2008; 129(12):124105. doi: 10.1063/1.2978177 [Crossref]
28. Zwanzig RW. High‐temperature equation of state by a perturbation method I Nonpolar gases. J Chem Phys 1954; 22(8):1420-6. doi: 10.1063/1.1740409 [Crossref]
29. Bennett CH. Efficient estimation of free energy differences from Monte Carlo data. J Comput Phys 1976; 22(2):245-68. doi: 10.1016/0021-9991(76)90078-4 [Crossref]
30. Klimovich PV, Shirts MR, Mobley DL. Guidelines for the analysis of free energy calculations. J Comput Aided Mol Des 2015; 29(5):397-411. doi: 10.1007/s10822-015-9840-9 [Crossref]
31. Rocklin GJ, Mobley DL, Dill KA, Hünenberger PH. Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects. J Chem Phys 2013; 139(18):184103. doi: 10.1063/1.4826261 [Crossref]
32. Callenberg KM, Choudhary OP, de Forest GL, Gohara DW, Baker NA, Grabe M. APBSmem: a graphical interface for electrostatic calculations at the membrane. PLoS One 2010; 5(9). doi: 10.1371/journal.pone.0012722 [Crossref]
33. Baker NA, Sept D, Joseph S, Holst MJ, McCammon JA. Electrostatics of nanosystems: application to microtubules and the ribosome. Proc Natl Acad Sci U S A 2001; 98(18):10037-41. doi: 10.1073/pnas.181342398 [Crossref]
34. Kucerka N, Tristram-Nagle S, Nagle JF. Structure of fully hydrated fluid phase lipid bilayers with monounsaturated chains. J Membr Biol 2005; 208(3):193-202. doi: 10.1007/s00232-005-7006-8 [Crossref]
35. Timofeeva TA, Sadykova GK, Rudneva IA, Boravleva EY, Gambaryan AS, Lomakina NF. Changes in the phenotypic properties of highly pathogenic influenza A virus of H5N1 subtype induced by N186I and N186T point mutations in hemagglutinin. Mol Biol 2016; 50(5):755-61. doi: 10.1134/s0026893316050174 [Crossref]
36. Hoffmann E, Neumann G, Kawaoka Y, Hobom G, Webster RG. A DNA transfection system for generation of influenza A virus from eight plasmids. Proc Natl Acad Sci U S A 2000; 97(11):6108-13. doi: 10.1073/pnas.100133697 [Crossref]
37. Killian ML. Hemagglutination assay for influenza virus. Methods Mol Biol 2014; 1161:3-9. doi: 10.1007/978-1-4939-0758-8_1 [Crossref]
38. Zotova SA, Korneeva TM, Shvedov VI, Fadeeva NI, Leneva IA, Fedyakina IT. Synthesis and antiviral activity of indole and benzofuran sulfides. Pharm Chem J 1995; 29(1):57-9. doi: 10.1007/bf02219467 [Crossref]
39. Burtseva EI, Shevchenko ES, Leneva IA, Merkulova LN, Oskerko TA, Shliapnikova OV. [Rimantadine and arbidol sensitivity of influenza viruses that caused epidemic morbidity rise in Russia in the 2004-2005 season]. Vopr Virusol 2007; 52(2):24-9.
40. Dong G, Peng C, Luo J, Wang C, Han L, Wu B. Adamantane-resistant influenza A viruses in the world (1902-2013): frequency and distribution of M2 gene mutations. PLoS One 2015; 10(3):e0119115. doi: 10.1371/journal.pone.0119115 [Crossref]
41. Wang C, Takeuchi K, Pinto LH, Lamb RA. Ion channel activity of influenza A virus M2 protein: characterization of the amantadine block. J Virol 1993; 67(9):5585-94. doi: 10.1128/jvi.67.9.5585-5594.1993 [Crossref]
42. Ma C, Zhang J, Wang J. Pharmacological characterization of the spectrum of antiviral activity and genetic barrier to drug resistance of M2-S31N channel blockers. Mol Pharmacol 2016; 90(3):188-98. doi: 10.1124/mol.116.105346 [Crossref]
43. Pielak RM, Chou JJ. Flu channel drug resistance: a tale of two sites. Protein Cell 2010; 1(3):246-58. doi: 10.1007/s13238-010-0025-y [Crossref]
44. Pielak RM, Schnell JR, Chou JJ. Mechanism of drug inhibition and drug resistance of influenza A M2 channel. Proc Natl Acad Sci U S A 2009; 106(18):7379-84. doi: 10.1073/pnas.0902548106 [Crossref]
45. Duff KC, Gilchrist PJ, Saxena AM, Bradshaw JP. Neutron diffraction reveals the site of amantadine blockade in the influenza A M2 ion channel. Virology 1994; 202(1):287-93. doi: 10.1006/viro.1994.1345 [Crossref]
46. Zeberezhnyĭ AD, Grebennikova TV, Vorkunova GK, Yuzhakov AG, Kostina LV, Norkina SN. [Engineering by reverse genetics and characterization of the new reassortant influenza virus strain H5N1]. Vopr Virusol 2014; 59(6):23-7.
47. Almlöf M, Brandsdal BO, Aqvist J. Binding affinity prediction with different force fields: examination of the linear interaction energy method. J Comput Chem 2004; 25(10):1242-54. doi: 10.1002/jcc.20047 [Crossref]