Identification of novel inhibitor against endonuclease subunit of Influenza pH1N1 polymerase: A combined molecular docking, molecular dynamics, MMPBSA, QMMM and ADME studies to combat influenza A viruses

Accepted Manuscript Title: Identification of novel inhibitor against endonuclease subunit of Influenza pH1N1 Polymerase: A combined molecular docking, molecular dynamics, MMPBSA and ADME studies to combat influenza A viruses Authors: Sajad Mohseni, Fariborz Nasri, Kambiz Davari, Sako Mirzaie, Atousa Moradzadegan, Frahad Farzaneh PII: DOI: Reference:

S1476-9271(18)30032-X https://doi.org/10.1016/j.compbiolchem.2018.08.005 CBAC 6906

To appear in:

Computational Biology and Chemistry

Received date: Revised date: Accepted date:

14-1-2018 6-8-2018 10-8-2018

Please cite this article as: Mohseni S, Nasri F, Davari K, Mirzaie S, Moradzadegan A, Farzaneh F, Identification of novel inhibitor against endonuclease subunit of Influenza pH1N1 Polymerase: A combined molecular docking, molecular dynamics, MMPBSA and ADME studies to combat influenza A viruses, Computational Biology and Chemistry (2018), https://doi.org/10.1016/j.compbiolchem.2018.08.005 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

SC RI PT

Identification of novel inhibitor against endonuclease subunit of Influenza pH1N1 Polymerase: a combined molecular docking, molecular dynamics, MMPBSA and ADME studies to combat influenza A viruses Sajad Mohseni1, Fariborz Nasri2, Kambiz Davari3, Sako Mirzaie1*, Atousa Moradzadegan4*, Frahad Farzaneh1

Department of Biochemistry, Sanandaj Branch, Islamic Azad University, Sanandaj, Iran

2

Department of Chemistry, Sanandaj Branch, Islamic Azad University, Sanandaj, Iran

3

Department of Microbiology, Sanandaj Branch, Islamic Azad University, Sanandaj, Iran

4

Department of Experimental Sciences, Dezful Branch, Islamic Azad University, Dezful, Iran

M

A

N

U

1

A

CC

EP

TE

D

Corresponding Authors: Atousa Moradzadegan, Email: [email protected] Sako Mirzaie, Email: [email protected] , [email protected]

GRAPHICAL ABSTRACT

1

SC RI PT U N A M D

HIGHLIGHTS

TE

1. The new potent inhibitor against endonuclease subunit of Influenza pH1N1 Polymerase (PA) is suggested by virtual screening, QMMM and MD studies.

EP

2. The inhibition mechanism and physico-chemical properties of three systems including free PA, PA: R05-2 complex and PA: ZINC15340668 complex were evaluated by MD simulations.

A

.

CC

3. ADMET analysis data showed that compound ZINC15340668 could be a potential inhibitor of therapeutic targets of PA and further analysis can be carried out through various experimental studies.

Abstract The influenza H1N1 virus is the causative agent of the flu pandemic in the world. Due to the shortage of effective means of control, it is remained the serious threats to public and avian health. To battle the surge of viral outbreaks, new treatments are crucially needed. The viral RNA polymerase, which is responsible for transcription and replication of the RNA genome, is comprised of subunits PA, PB1 and PB2. PA has endonuclease activity and is a well known 2

SC RI PT

target for inhibitor and drug design. In the current study, we employed molecular docking, molecular dynamics (MD), MMPBSA and ADME studies to find and propose an inhibitor among 11873 structures against PA. Our molecular docking, MD and MMPBSA studies showed that ZINC15340668 has ideal characteristics as a potent PA inhibitor and can be used in experimental phase and further development. Also, ADME prediction demonstrated that all physic-chemical parameters are within the acceptable range defined for human use. Molecular mechanism based study revealed that upon inhibitor binding; the flexibility of PA backbone is increased. This observation demonstrates the plasticity of PA active site, and it should be noticed in drug design against PA Influenza A viruses. Keywords: Influenza H1N1, Endonuclease, Molecular docking, Molecular dynamics, ADME prediction

U

1- Introduction

A

CC

EP

TE

D

M

A

N

The common and ordinary causes of seasonal epidemics are influenza viruses. This virus can cause wide symptoms from simple cough, sore throat, nasal discharge, fever, headache, and muscle pain, to more sever and complicated conditions like bronchitis and pneumonia [1]. Just in the 20th century, four influenza pandemics were occurred and recorded. In pandemic “Spanish flu” in 1918, 500 million people were infected by flu, which led to over 30 million deaths [2]. In 2009, a pandemic “swine flu” affected human [3]. Additionally, “avian flu” shows an ongoing danger that may result in tragedy consequences if not controlled and limited [1]. A nucleus of infected cells is a site of virus genome replication. A viral RNA-dependent RNA polymerase, which is responsible for the replication and transcription, is comprised of subunits PA, PB1 and PB2 [4]. The mechanism of viral mRNA transcription is known as “cap-snatching” mechanism [5]. Up to now, this unusual mechanism has been reported just in negative strand, segmented RNA viruses, including orthomyxoviruses (notably influenza), bunyaviruses and arenaviruses [4]. In the case of influenza, cap-snatching initiates with the binding of host cell pre-mRNAs due to their 5ʹ cap structure to the PB2 subunit of the polymerase. Then, the first 10–13 nucleotides of the pre-mRNAs from 5ʹ end cleaved by the endonuclease activity of PA subunit (PA) reside in viral polymerase. Now, the resulted capped oligomer is used as primers for transcription of the viral mRNAs by the PB1 subunit of the polymerase [4]. Finally, the viral transcript is polyadenylated at a conserved U-rich region of the template RNA [6]. The resulted viral mRNA has both the 5ʹ and 3ʹ signals to translate in the host cell cytosol. There is some efforts in producing the crystal structure of the two functional domains involved in cap-snatching [7]. The PB2 subunit, which contains a cap-binding domain, binds to the 7-methyl guanosine ligand by means of an aromatic sandwich [8]. PA includes endonuclease domain at its N-terminus end and has a core fold like other two-metal dependent nucleases of the PD…D/E…K superfamily [9, 10]. Based on mutation studies, it has been shown that the endonuclease domain has high affinity to bind to divalent cation, manganese ions. This enzyme also has lower affinity to magnesium ions 3

M

A

N

U

SC RI PT

[11]. Transcription by cap-snatching is vital and necessary for virus replication and so, the inhibition of endonuclease domain is a good target for drug design against influenza viruses [12, 13]. With regard to the on-going circulation of highly pathogenic avian H5N1 strains, which could potentially change for human-to-human transmission [14], the unexpected appearance of the 2009 H1N1 pandemic strain [15], and the emergence of resistance to current and common treatments, the need for new treatment targeting influenza virus is now widely acknowledged [4]. Since, the inhibition of the cap-snatching endonuclease of influenza virus polymerase stall viral transcription and replication, this domain has been assigned as a target for drug design for approximately 15 years [4]. A number of 4-substituted 2,4-dioxobutanoic acid compounds [1618], a substituted 2,6-diketopiperazine natural compound (Flutimide) from the fungus Delitschia confertaspora [19] along with derivatives, and N-hydroxamic acid and N-hydroxyimide compounds are examples of designed and synthesized inhibitors against viral influenza virus endonuclease [20]. MD simulation methods have become a favorable approach to design new drugs. Nowadays large macromolecular systems with many ingredients can be simulating in the nano and microsecond scale [21]. Drug discovery through computational methods speeds up designing, optimizing and discovering the new drugs. By MD simulation, we can study the microscopic behavior of biological systems during the time. Classical MD, in relative to Newton’s physics, is a method for investigating the position exchange of atoms and molecules [22]. In the current study, we employed virtual screening approach to rational design the inhibitors from a large library of small molecules, against the endonuclease domain of influenza pH1N1 polymerase.

TE

D

2- Methods

2-1- Enzyme selection and preparation

A

CC

EP

The crystallographic structure of PA from influenza pH1N1 polymerase (PDB ID: 4AVG) was retrieved from the Protein Data Bank (PDB). Before docking, enzyme structure was prepared by removing water molecules and then, bond orders were assigned, and hydrogen atoms were added to the crystal structure. Gasteiger–Marsili method [23] was assigned as a charge method for the molecular docking study. For accurate molecular docking and to prevent clashes, the structure was minimized by the steepest descent method implemented in GROMACS 5.0.2 [24]. A cocrystallized inhibitor, 4-[4-[(4-chlorophenyl) methyl]-1-(cyclohexylmethyl)-4-piperidinyl]- 2hydroxy-4-oxo-2-butenoic acid (R05-2) [4], was chosen as a positive control and reference for the virtual screening (Figure 1). This molecule is the inhibitor of PA. 2-2- Library selection and preparation As mentioned in the previous section, R05-2 was selected as a positive control and compounds with 40 percent resemblance to this molecule, were obtained from ZINC database (totally 11873 structures) [25]. Gasteiger-Marsili method was used for calculating the charge of molecules. 4

Lastly, for molecular docking studies, the 3D structures were saved in pdb format. For the internal validation phase, the structure of R05-2 was extracted from the PDB file using program ViewerLite 5.0. After assigning bond orders, missing hydrogen atoms were added. Then, in the AutoDockTools package, the partial atomic charges were calculated using the Gasteiger–Marsili method.

SC RI PT

2-3- Molecular docking and virtual screening studies

U

AutoDock Vina [26] was used for molecular docking of the compounds of ZINC library (11873 molecules) and the active site of PA. This software which is free for academic researchers performs virtual screening and molecular docking by the experimental scoring functions. AutoDock Vina automatically determines the grid maps [26]. Default parameters of software were appropriated for the docking scope. For all docking runs, a grid map with 30×30×30 points and grid-point spacing of 1 Å was set up at the geometrical center of the R05-2 in the PDB structure. AutoDock Vina produces several poses with different binding energy for each compound and the pose with lowest energy was selected for MD simulations.

N

2-4- Molecular dynamics simulations

A

CC

EP

TE

D

M

A

All MD simulations were performed using GROMACS 5.0.2 program and AMBER 99SB force field [24]. We used MD for performing protein and ligand complexes as described [27-29].The pKa of PA residues was determined by PROPKA 2.0 server to select the appropriate ionization states of PA ionizable residues [30]. The topology files and partial charges of inhibitors were created by ACPYPE, a tool predicated on ANTECHAMBER [31]. [31]. Each system solvated in a cubic regular box with TIP3P water model [32] and enough sodium and/or chlorine ions were included to neutralize each system. Energy minimizations were achieved by the steepest descent integrator as well as the conjugate gradient algorithm to reach a maximum force of less than 1000 kJ mol-1 nm-1 on each atom. A twin-range cutoff scheme was designated to estimate shortrange, non-bonded interactions, with van der Waals interactions shortened at 1.4 nm and electrostatic interactions shortened at 0.9 nm. To deal with the Long-range electrostatic interactions, particle mesh Ewald (PME) technique was implemented [33, 34]. The temperature was fixed at 300 K (ref_t) by velocity rescaling associated with a stochastic term and coupling time constant (tau_t) of 0.1 ps [35]. Such thermostat can be compared with Berendsen coupling, with the same scaling utilizing tau_t, however the stochastic term certifies that a suitable canonical joint is constructed [24]. The pressure was set at 1.0 atm by using a Parrinello-Rahman barostat with a coupling constant of 2 ps [36]. Simulations were accomplished within a time period of 2 fs, and all bonds engaging hydrogen atoms were confined using a Linear Constraint Solver (LINCS). The Number of iterations to adjust for rotational lengthening in LINCS (lincs_iter) and highest order in the extension of the constraint coupling matrix (lincs_order) was allocated 1 and 4, in that order. In order to create "normal" MD simulations, an order of 4 is normally sufficient. All systems were equilibrated under a constant pressure (NPT) ensemble (100 ps) and a constant volume (NVT) ensemble (100 ps). Each MD simulation was carried out 5

Bi = 8π2Ui2

Eq. (1)

Where, Ui2 is the mean square displacement of atom i [41]. 2-5- Calculation of Inhibitor binding free energy

SC RI PT

for 30 ns. The analyses of trajectories were performed utilizing VMD software and the standard tools applied in the GROMACS program [37]. Anomalous scattering confirmed the presence of two bound manganese ions in the active site of PA [4]. The charge and non-bonded parameters for manganese are not defined in any force filed implemented in GROMACS. Based on Mulliken charges calculated using GAUSSIAN94 [38] and a 6-31G* basis set, a full +2 charge assigned for manganese ion. In addition, for the non-bonded parameters of Mn+2, a well depth of -0.014 kcal mol-1 [39] and van der Waals radii of 1.69 Å were employed [40]. Experimental RMSF values about the mean position of atom i were estimated from the crystallographic temperature factor Bi with Eq. (1).

M

Gbinding = Gcomplex – (Gprotein + Gligand)

A

N

U

Modified MM-PBSA or Molecular Mechanics–Poisson Boltzmann Surface Area is open source software, utilized for obtaining the binding free energy among receptor and inhibitors. As a scoring function, the MM-PBSA has been exerted for the drug design by computational methods [42]. MM-PBSA in this study was used to achieve the binding free energy of designed inhibitor: PA and R05-2: PA, separately. Following equation describes the free energy of binding:

TE

2-6- ADME prediction

D

Free energy of the protein–inhibitor complex showed by Gcomplex, Gprotein represents the free energies of the protein in solvent and Gligand is free energies of the inhibitor in solvent [42].

A

CC

EP

The absorption, distribution, metabolism and excretion (ADME) parameters of all selected compounds from NSC were computed using QikProp module. QikProp is a quick, accurate, easy-to-use absorption, distribution, metabolism, and excretion (ADME) prediction program designed by Professor William L. Jorgensen [43]. This software, which is implemented in Schrödinger’s Maestro Molecular modeling suit, predicts physically significant descriptors and pharmaceutically relevant properties of organic molecules, either individually or in batches.

3- Results and discussion 3-1- Internal validation phase and docking studies Before run the virtual screening, validation phase was done to evaluate the precise of molecular docking. In the internal validation phase, co-crystallized R05-2 was docked onto the PA (PDB: 6

N

U

SC RI PT

4AVG) and then, the highest racking pose was selected. Superimposing the experimental and predicted conformation expressed as root mean square deviation (RMSD). The RMSD for R05-2 was 0.98, which showed that the employed in silico method is robust and appropriate for estimating the interactions of such ligand with PA. A molecular surface view of the active pocket of PA with its co-crystallized inhibitor (R05-2) is demonstrated in Figure 2 and shows the similarity between experimental structure and docked pose [44]. Afterwards, the identification of lead compounds using AutoDock with ZINC library was done. Compound R05-2, which is employed as positive control of PA, is a relatively potent inhibitor with the IC50 of 0.33 µM [18]. With regard to the docking energy, the ZINC codes of the best ten structures, together with positive control were summarized in Table 1. Also, their 2D structures are depicted in Figure 3. Based on the docking energy (Table 1), just two structures, with ZINC codes of ZINC15340668 and ZINC17323185, have lower docking energy than R05-2. The binding poses of ZINC15340668 and ZINC17323185 in the active site of are shown in Figure 4. ZINC15340668 has the lowest docking energy and interacts with Ala 20, Met 21 and Tyr 24, Ile 138, Phe 105 via hydrophobic interaction. This molecule also coordinated with one of the divalent Mn located in the PA (Figure 4). 3-2- MD simulations of PA and validation

A

CC

EP

TE

D

M

A

After the molecular docking studies and virtual screening, the screened molecule with the lowest docking energy in complex with PA, along with R05-2: enzyme complex and free enzyme were introduced to MD simulations (30 ns). By performing the MD simulations of bio-molecules, the conformational changes, formed and disrupted interactions and their behaviors during a defined time could be analyzed in details. For all three systems mentioned above, structural root mean square deviations (RMSD) of backbone atoms were calculated during 30 ns to evaluate the stability of the MD simulations (Figure 5). The RMSD is calculated based on superimposed structures fitting to the first frame of the 30 ns MD simulation. As can be seen in Figure 5, the average RMSD value for free PA is 1.9 Å. Upon ZINC15340668 binding, the RMSD of enzyme was increased to 2.5 Å at 3610 ps. A sharp rise in RMSD value was happened at 26300 ps to 27640 ps (RMSD ~ 3.7 Å). The increasing of RMSD value in the case of R05-2 was occurred at 21820 ps from 1.5 to 3.2 Å. At the last 5000 ps, both ZINC15340668 and R05-2 bound PA have the same fashion in their RMSD profile (Figure 5). With regard to Figure 5, the fluctuation of backbone atoms of PA was increased upon ZINC15340668 binding. . In crystal structure of proteins, there is an important factor which is known as B-factor, temperature factor or DebyeWaller factor. It measures the ambiguity/mobility of an atom in dynamic protein 3D structures, namely, the displacement of the atomic positions from its mean position [45]. B factor increases as Ui2 increases. A low B factor demonstrates that the atom is in the well-ordered parts of the structure, while a large B factor generally implies a very high flexibility of this atom. We computed the experimental RMSF (root mean-square fluctuation) of PA crystal by employing Eq. (1) to validate our MD simulations (Figure 6). The last two residues were excluded from our validation phase because of the high flexibility of these two residues possibly result from the 7

A

CC

EP

TE

D

M

A

N

U

SC RI PT

poor ordered of residues 196-198 [4]. As can be seen in Eq. (1), the effects of static disorder in crystal are ignored [46]. Based on Figure 9, the overall pattern of fluctuations captured from the MD simulation is close to that observed from the crystallographic B-factors (Figure 6). However, there are some distinct differences between theoretical and experimental RMSF. Some residues in the PA crystal have lower RMSF than the values obtained by MD simulation. This effect may be due to the result from the crystal packing arrangement, which restricts the flexibility of certain regions. Figure 7 exhibits the RMSF values of all residues of the free PA, R05-2 and ZINC15340668 bound PA. A region with the highest RMSF value is located between residues Ile 62 and Ala 70. In the free PA, the RMSF value of Ile 62 is 3.2 Å, while for R05-2: PA and ZINC15340668: PA complexes, these values are 5 Å and 4.4 Å, respectively. The region mentioned above, is a part of a mobile loop (residues 53–73) of unknown function that is solvent exposed and usually disordered [4]. Glu 59, in the crystal form, interacts with a divalent cation belonging to the active site of a neighboring PA [9]. It has been shown that the replacement of residues 51–72 with a Gly-Gly-Ser does not significantly affect enzymatic activity [47]. Figure 7 demonstrates that upon R05-2 or ZINC15340668 binding to the PA active site, the RMSF values of residues 62-70 are sharply increased. The RMSF value of Lys 34 in the free PA is 0.95 Å. Upon binding of ZINC15340668 and R05-2 into the PA, this value for Lys 34 is elevated to 1.3 and 1 Å, respectively (Figure 7). Upon inhibitor binding, increasing in the RMSF values of residues Tyr 24, Ile 38 and Arg 84 were also observed. Residue Arg 84 is highly conserved among all influenza A strains [11, 48]. It has been shown that upon ligand binding into the pocket of PA, the orientations of side chains of residues Arg84, Lys34, Glu26 and Phe105 are changed. It shows a plasticity of the active site and an induced fit mode of ligand binding [4]. This elucidation has a good agreement with our RMSF data in Figure 7. For clarification of the increasing of the RMSF value of Lys 34 in three systems, Figure 8 is provided from the last frame of the MD simulation. In the free PA, Lys 34 is fixed by 5 hydrogen bonds (Figure 8A). The backbone NH of Lys 34 is interacts with the backbone carbonyl moiety of Glu 31 via hydrogen bond. Also, the side chain carboxyl moiety of Glu 26 interacts with the side chain NH3 of Lys 34 due to two hydrogen bonds. In addition, the backbone carbonyl group of Lys 34 makes two hydrogen bonds with the NH of Ala 37 and Ile 38 (Figure 8A). However, Lys 34 in both inhibitor:PA complexes, is fixed by four hydrogen bonds (Figure 8B and 8C). This fewer one hydrogen bond could be the reason of the higher flexibility of Lys 34 in the inhibitor bound PA. In the R05-2 bound PA, the backbone NH of Ile 38 makes a hydrogen bond with the carbonyl group of Lys 34. The side chain NH3 moiety of Lys 34 interacts with the side chain carboxyl moieties of Glu 26 and Glu 31 by means of two hydrogen bonds. The last hydrogen bond is formed between the backbone NH of Lys 34 and the backbone carbonyl of Glu 31 (Figure 8B). In the ZINC15340668 bound PA, the side chain NH3 of Lys 34 interacts with the side chain carboxyl group of Glu 26. In the same manner, the backbone NH of Lys 34 makes a hydrogen bond with the carbonyl moiety of Glu 31. In addition, the backbone carbonyl of Lys 34 interacts with the backbone NH of Ala 37 and Ile 38 via two hydrogen bonds (Figure 8C). Prior studies showed that Tyr 24, Glu 26, Lys 34, Ala37, Ile38, Arg84, Phe105, Tyr130 and Lys134 are highly 8

N

U

SC RI PT

conserved amongst all influenza A strains [4]. So, the inhibitor, which is proposed here, could be employed against all influenza A strains. The DSSP secondary structure analysis was carried out on the trajectory file in the term of simulation time as presented in Figure 9. In the free PA, the first turn (residue 25-28) is appeared after starting simulation time and keeps it till 28500 ps (Figure 9A). Then it is converted to the α-helix. In the ZINC15340668 bound PA (Figure 9B) and R05-2 bound PA (Figure 9C), this region shifts between turn and α-helix structures. Residues Tyr 24 and Glu 24 are located in this flexible turn and bridging the second and third αhelix structures. Upon Figure 8A, Lys 134, as an important catalytic residue, is located in a fifth α-helix. This α-helix is appeared after beginning simulation time and keeps it for a short time. Then it shifts between α-helix and turn, and continues to 22000 ps. Finally it shifts between turn and 310-helix. In the ZINC15340668 bound PA, after a while, the region containing Lys 134 fluctuates around turn, 310-helix and bend structures. After 22000 ps, it just shifts between turn and bend structure (Figure 9B). It seems this variation and alteration in the secondary structure, is responsible for the high RMSF value of Lys 134 (2Å) in comparison with free PA (RMSF of 1.3 Å) and R05-2 bound PA (RMSF of 0.7 Å) (Figure 7). In the R05-2 bound PA, the helix containing Lys 134 shifts between α-helix and turn structures. 3-3- Calculation of Binding free energy by MMPBSA

CC

EP

TE

D

M

A

The MM/PBSA methodology was integrated with MD simulations to estimate the binding-free energy of docked complexes, in order to increase the precision of prediction and explore detailed energy composition [29]. The snapshots were extracted from the last 10 ns of MD trajectories for the analysis of the binding-free energy. The binding-free energies (ΔGbind) of R05-2 and ZINC15340668 with YRS were calculated at -9.8 and -35.4 kJ mol-1, respectively (Table 2). As seen in Table 2, both inhibitors have approximately similar solvent-accessible surface area energy. The Van der Waals energy of ZINC15340668: PA system (-165.6 kJ mol-1) is lower than R05-2: PA complex (-130 kJ mol-1). The electrostatic energy of R05-2 is higher than ZINC15340668. However, our data summarized in Table 2, demonstrates that the optimizations of electrostatic interactions between the ZINC15340668 and PA may lead to the potent inhibitor. The challenging matter is that with the decreasing the electrostatic energy, the salvation energy could be increased and the total binding energy may be remain unchanged. 3-4- Theoretical ADME estimations

A

The selected compounds (two high ranked compounds in Table 1) and our positive control (R052) were further analyzed for their drug-like behavior by analysis of pharmacokinetic parameters required for ADMET by employing QikProp. For these compounds, Predicted octanol/water partition coefficient, Predicted blockage of hERG K + channel [49], Predicted apparent Caco-2 cell permeability in nm/s [50], Percentage of human oral absorption and Predicted central nervous system activity [51], were predicted and included in Table 3. With regard to the ADMET prediction results, the pharmacokinetic properties of the hits are within the acceptable range (Table 3). Based on pharmacokinetic parameters, ZINC15340668 has better human oral 9

absorption than R05-2. Also, R05-2 has deviation in Jorgensen Rule (Table 3). Our computational pharmacokinetic studies showed that ZINC15340668 is suitable for additional and further drug development process.

4- Conclusion

A

CC

EP

TE

D

M

A

N

U

SC RI PT

In this study, we occupied molecular docking, MD, MM/PBSA and ADME studies to identify potent hits against PA from Influenza A. Compound R05-2 was selected as a positive control and the structures with the similarity of 40 percent with the control were obtained from ZINC database. With regard to docking energy, the structure with the lowest docking energy (ZINC15340668) was chosen as the best hit. To elucidate the inhibition mechanism of PA in detail, three systems including free PA, R05-2 bound PA and ZINC15340668 bound PA were introduced to MD simulations. The RMSD analysis on the backbone of PA showed that upon inhibitor binding (R05-2 or ZINC15340668), the flexibility of enzyme is increased. RMSF analysis demonstrated that in the inhibitor bound PA, the flexibilities of Lys 34, Tyr 24, Ile 38 and Arg 84 are raised. The elevated flexibility of these residues is in correlated with the plasticity of the active site and an induced fit mode of ligand binding. To validate our MD simulations, we converted the B-factor of Cα atoms of PA crystal, to the RMSF values. Then, we compare them with our computational RMSF. The comparison showed a good correlation between the experimental values and theoretical RMSF. Based on MMPBSA computation, the free energy of binding of ZINC15340668 to the PA active site, is lower than R05-2. This analysis also revealed that the optimizations of electrostatic interactions between the ZINC15340668 and PA would lead to the potent inhibitor. Upon theoretical ADME data, ZINC15340668 is appropriate for additional and further drug development process.

10

References

A

CC

EP

TE

D

M

A

N

U

SC RI PT

[1] Shim JM, Kim J, Tenson T, Min J-Y, Kainov DE. Influenza Virus Infection, Interferon Response, Viral Counter-Response, and Apoptosis. Viruses 2017;9:223. [2] Taubenberger JK, Morens DM. 1918 Influenza: the mother of all pandemics. Rev Biomed 2006;17:69-79. [3] Fineberg HV. Pandemic preparedness and response—lessons from the H1N1 influenza of 2009. New England Journal of Medicine 2014;370:1335-42. [4] Kowalinski E, Zubieta C, Wolkerstorfer A, Szolar OH, Ruigrok RW, Cusack S. Structural analysis of specific metal chelating inhibitor binding to the endonuclease domain of influenza pH1N1 (2009) polymerase. PLoS pathogens 2012;8:e1002831. [5] Plotch SJ, Bouloy M, Ulmanen I, Krug RM. A unique cap (m7GpppXm)-dependent influenza virion endonuclease cleaves capped RNAs to generate the primers that initiate viral RNA transcription. Cell 1981;23:847-58. [6] Poon LL, Pritlove DC, Fodor E, Brownlee GG. Direct evidence that the poly (A) tail of influenza A virus mRNA is synthesized by reiterative copying of a U track in the virion RNA template. Journal of virology 1999;73:3473-6. [7] Ruigrok RW, Crépin T, Hart DJ, Cusack S. Towards an atomic resolution understanding of the influenza virus replication machinery. Current opinion in structural biology 2010;20:104-13. [8] Guilligay D, Tarendeau F, Resa-Infante P, Coloma R, Crepin T, Sehr P, et al. The structural basis for cap binding by influenza virus polymerase subunit PB2. Nature structural & molecular biology 2008;15:500-6. [9] Dias A, Bouvier D, Crépin T, McCarthy AA, Hart DJ, Baudin F, et al. The cap-snatching endonuclease of influenza virus polymerase resides in the PA subunit. Nature 2009;458:914. [10] Yuan P, Bartlam M, Lou Z, Chen S, Zhou J, He X, et al. Crystal structure of an avian influenza polymerase PA(N) reveals an endonuclease active site. Nature 2009;458:909-13. [11] Crépin T, Dias A, Palencia A, Swale C, Cusack S, Ruigrok RW. Mutational and metal binding analysis of the endonuclease domain of the influenza virus polymerase PA subunit. Journal of virology 2010;84:9096-104. [12] Das K, Aramini JM, Ma L-C, Krug RM, Arnold E. Structures of influenza A proteins and insights into antiviral drug targets. Nature structural & molecular biology 2010;17:530-8. [13] Krug RM, Aramini JM. Emerging antiviral targets for influenza A virus. Trends in pharmacological sciences 2009;30:269-77. [14] Palese P, Wang TT. H5N1 influenza viruses: facts, not fear. Proceedings of the National Academy of Sciences 2012;109:2211-3. [15] Itoh Y, Shinya K, Kiso M, Watanabe T, Sakoda Y, Hatta M, et al. In vitro and in vivo characterization of new swine-origin H1N1 influenza viruses. Nature 2009;460:1021. [16] Hastings J, Selnick H, Wolanski B, Tomassini J. Anti-influenza virus activities of 4substituted 2, 4-dioxobutanoic acid inhibitors. Antimicrobial agents and chemotherapy 1996;40:1304-7. [17] Nakazawa M, Kadowaki S-e, Watanabe I, Kadowaki Y, Takei M, Fukuda H. PA subunit of RNA polymerase as a promising target for anti-influenza virus agents. Antiviral research 2008;78:194-201. 11

A

CC

EP

TE

D

M

A

N

U

SC RI PT

[18] Tomassini J, Selnick H, Davies M, Armstrong M, Baldwin J, Bourgeois M, et al. Inhibition of cap (m7GpppXm)-dependent endonuclease of influenza virus by 4-substituted 2, 4dioxobutanoic acid compounds. Antimicrobial agents and chemotherapy 1994;38:2827-37. [19] Tomassini J, Davies M, Hastings J, Lingham R, Mojena M, Raghoobar S, et al. A novel antiviral agent which inhibits the endonuclease of influenza viruses. Antimicrobial agents and chemotherapy 1996;40:1189-93. [20] Cianci C, Chung T, Meanwell N, Putz H, Hagen M, Colonno R, et al. Identification of Nhydroxamic acid and N-hydroxyimide compounds that inhibit the influenza virus polymerase. Antiviral Chemistry and Chemotherapy 1996;7:353-60. [21] Mortier J, Rakers C, Bermudez M, Murgueitio MS, Riniker S, Wolber G. The impact of molecular dynamics on drug design: applications for the characterization of ligand– macromolecule complexes. Drug Discovery Today 2015;20:686-702. [22] De Vivo M, Masetti M, Bottegoni G, Cavalli A. Role of molecular dynamics and related methods in drug discovery. Journal of medicinal chemistry 2016;59:4035-61. [23] Gasteiger J, Marsili M. Iterative partial equalization of orbital electronegativity—a rapid access to atomic charges. Tetrahedron 1980;36:3219-28. [24] Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJ. GROMACS: fast, flexible, and free. Journal of computational chemistry 2005;26:1701-18. [25] Irwin JJ, Shoichet BK. ZINC− a free database of commercially available compounds for virtual screening. Journal of chemical information and modeling 2005;45:177-82. [26] Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. Journal of computational chemistry 2010;31:455-61. [27] Davari K, Nowroozi J, Hosseini F, Sepahy AA, Mirzaie S. Structure-based virtual screening to identify the beta-lactamase CTX-M-9 inhibitors: An in silico effort to overcome antibiotic resistance in E. coli. Computational Biology and Chemistry 2017;67:174-81. [28] Hosseini Y, Mollica A, Mirzaie S. Structure-based virtual screening efforts against HIV-1 reverse transcriptase to introduce the new potent non-nucleoside reverse transcriptase inhibitor. Journal of Molecular Structure 2016;1125:592-600. [29] Karami M, Jalali C, Mirzaie S. Combined virtual screening, MMPBSA, molecular docking and dynamics studies against deadly anthrax: An in silico effort to inhibit Bacillus anthracis nucleoside hydrolase. Journal of Theoretical Biology 2017;420:180-9. [30] Bas DC, Rogers DM, Jensen JH. Very fast prediction and rationalization of pKa values for protein–ligand complexes. Proteins: Structure, Function, and Bioinformatics 2008;73:765-83. [31] Wang J, Wang W, Kollman PA, Case DA. Automatic atom type and bond type perception in molecular mechanical calculations. Journal of molecular graphics and modelling 2006;25:24760. [32] Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. Comparison of simple potential functions for simulating liquid water. The Journal of chemical physics 1983;79:926-35. [33] Darden T, York D, Pedersen L. Particle mesh Ewald: An N⋅ log (N) method for Ewald sums in large systems. The Journal of chemical physics 1993;98:10089-92. [34] Essmann U, Perera L, Berkowitz ML, Darden T, Lee H, Pedersen LG. A smooth particle mesh Ewald method. The Journal of chemical physics 1995;103:8577-93. [35] Bussi G, Donadio D, Parrinello M. Canonical sampling through velocity rescaling. The Journal of chemical physics 2007;126:014101.

12

A

CC

EP

TE

D

M

A

N

U

SC RI PT

[36] Parrinello M, Rahman A. Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied physics 1981;52:7182-90. [37] Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics. Journal of molecular graphics 1996;14:33-8. [38] Frisch M, Trucks G, Schlegel H, Gill P, Johnson B, Robb M, et al. Gaussian 94, revision B. 3. Gaussian Inc, Pittsburgh, PA 1995;9. [39] Vedani A, Huhta DW. A new force field for modeling metalloproteins. Journal of the American Chemical Society 1990;112:4759-67. [40] Bradbrook GM, Gleichmann T, Harrop SJ, Habash J, Raftery J, Kalb J, et al. X-Ray and molecular dynamics studies of concanavalin-A glucoside and mannoside complexes Relating structure to thermodynamics of binding. Journal of the Chemical Society, Faraday Transactions 1998;94:1603-11. [41] Kuzmanic A, Zagrovic B. Determination of ensemble-average pairwise root mean-square deviation from experimental B-factors. Biophysical journal 2010;98:861-71. [42] Kumari R, Kumar R, Consortium OSDD, Lynn A. g_mmpbsa_ A GROMACS tool for highthroughput MM-PBSA calculations. Journal of chemical information and modeling 2014;54:1951-62. [43] Jorgensen WL, Maxwell DS, Tirado-Rives J. Development and testing of the OPLS allatom force field on conformational energetics and properties of organic liquids. J Am Chem Soc 1996;118:11225-36. [44] Mirzaie S, Rafii F, Yasunaga K, Yoshunaga K, Sepehrizadeh Z, Kanno S, et al. Prediction of the mode of interaction between monoterpenes and the nitroreductase from Enterobacter cloacae by docking simulation. Computers in biology and medicine 2012;42:414-21. [45] Liu Q, Li Z, Li J. Use B-factor related features for accurate classification between protein binding interfaces and crystal packing contacts. BMC bioinformatics 2014;15:S3. [46] Kapustina M, Carter CW. Computational studies of tryptophanyl-tRNA synthetase: activation of ATP by induced-fit. Journal of molecular biology 2006;362:1159-80. [47] DuBois RM, Slavish PJ, Baughman BM, Yun M-K, Bao J, Webby RJ, et al. Structural and biochemical basis for development of influenza virus inhibitors targeting the PA endonuclease. PLoS pathogens 2012;8:e1002830. [48] Yuan S, Chu H, Singh K, Zhao H, Zhang K, Kao RY, et al. A novel small-molecule inhibitor of influenza A virus acts by suppressing PA endonuclease activity of the viral polymerase. Scientific reports 2016;6:22880. [49] Cavalli A, Poluzzi E, De Ponti F, Recanatini M. Toward a pharmacophore for drugs inducing the long QT syndrome: insights from a CoMFA study of HERG K+ channel blockers. Journal of medicinal chemistry 2002;45:3844-53. [50] Yazdanian M, Glynn SL, Wright JL, Hawi A. Correlating partitioning and Caco-2 cell permeability of structurally diverse small molecular weight compounds. Pharmaceutical research 1998;15:1490-4. [51] Ajay, Bemis GW, Murcko MA. Designing libraries with CNS activity. Journal of medicinal chemistry 1999;42:4942-51.

13

Table 1. Zinc code, IUPAC name, molecular weight and docking energies of the best ten poses among 11873 structures, achieved by virtual screening. For comparison, the R05-2 has been included as positive control.

ZINC17323185

3

ZINC15741842

4

ZINC12977416

5

ZINC15619754

6

ZINC16645015

7

ZINC08086163

8

ZINC15619720

9

ZINC15619719

10

ZINC77342870

Positive control

R05-2

U

2

N

ZINC15340668

A

1

Molecular weight (g/mol) N-[(4-chlorophenyl)methyl]-4-oxo-4343.854 (4-propylphenyl)butanamide 4-(4-chlorophenyl)-N-[(1S)-1-(3350.245 chlorophenyl)ethyl]-4-oxo-butanamide 4-(4-chlorophenyl)-4-oxo-N-[(1R)-1-(p343.854 tolyl)propyl]butanamide 4-(4-chlorophenyl)-N-[(3336.218 chlorophenyl)methyl]-4-oxo-butanamide 4-(4-chlorophenyl)-N-[(1S)-1-(4350.245 chlorophenyl)ethyl]-4-oxo-butanamide 4-(4-chlorophenyl)-4-oxo-N-(p315.8 tolylmethyl)butanamide N-[1-(4-chlorophenyl)ethyl]-4-oxo-4315.8 phenyl-butanamide N-[(1S)-1-(4-chlorophenyl)ethyl]-4-oxo329.827 4-(p-tolyl)butanamide N-[(1R)-1-(4-chlorophenyl)ethyl]-4-oxo329.827 4-(p-tolyl)butanamide 4-(4-chlorophenyl)-N-[(1S)-2-hydroxy-1317.816 phenyl-ethyl]butanamide

Docking energy (kcal/mol) -16.2

SC RI PT

IUPAC name

M

ZINC code

D

Number

-14.2 -10.3 -9.6 -9.5 -9.1 -8.4 -8.1 -8.1 -7.9

A

CC

EP

TE

-11.8

14

Table 2. Molecular energy terms for the PA: R05-2 and PA: ZINC15340668 complexes.

vdw, van der Waals; elect, electrostatic; solv, polar solvation; SASA, solvent-accessible surface area.

PA: ZINC15340668

-130.086 90.637 49.245 -19.653 -9.857

-165.624 51.963 97.205 -18.972 -35.427

SC RI PT

PA: R05-2

A

CC

EP

TE

D

M

A

N

U

Energy (kJ/mol) ΔEvdw ΔEelect ΔEsolv ΔESASA ΔGbinding

15

Table 3. Physico-chemical properties of the hits obtained from Qikprop of Schrödinger.

QPlogPo/wa

hERGb

QPPCacoc

QPlogBBd

ZINC15340668 ZINC17323185 R05-2

4.170 3.597 2.448

-4.533 -4.613 -3.811

1755 2330 20

-0.345 -0.167 -0.917

Percent Human Oral Absorptione 100 100 65

Predicted octanol/water partition coefficient (reasonable value from -2.0 to 6.5).

b

Predicted blockage of hERG K + channel (concern below -6).

c

Predicted apparent Caco-2 cell permeability in nm/s (<25 poor, >500 great).

d

Predicted brain/blood partition coefficient (reasonable value from -3.0 to 1.2).

e

Percentage of human oral absorption (<25% is poor).

f

Predicted central nervous system activity on a –2 (inactive) to +2 (active) scale.

A

CC

EP

TE

D

M

A

N

U

a

16

Lipinski Rule of 5 Violations

Jorgensen Rule of 3 Violations

CNSf

0 0 0

0 0 1

+/+/-

SC RI PT

Molecule ID

Figure legends Figure 1. The structure of R05-2, which was selected as a template and positive control in virtual screening.

SC RI PT

Figure 2. Superposition of crystal (red stick representation) and docked pose (blue stick representation) of R05-2 in the active site of PA to validate the molecular docking study. The Mn ions are shown as purple sphere. Figure 3. The 2D presentation of ZINC codes of the best ten structures in molecular docking.

Figure 4. Schematic representations of the orientation of ZINC15340668 and ZINC17323185 in the active site of PA.

U

Figure 5. Temporal RMSD values of free PA (blue line), PA: ZINC15340668 (red line) and PA: R05-2 (green line).

A

N

Figure 6. Comparison of RMSF of Cα atoms per residue extracted from the R05-2-bound PA simulation and crystallographic data of 4AVG.

M

Figure 7. RMSF plots for the MD simulations. RMSF comparison between free PA (blue line), PA: ZINC15340668 (red line) and PA: R05-2 (green line).

TE

D

Figure 8. Schematic representations of the conformation of Lys 34 in PA in three states: free PA (A), PA: R05-2mcomplex (B) and PA: ZINC15340668 complex (C). data are extracted from the last frame of MD simulation.

A

CC

EP

Figure 9. DSSP analysis of free PA (A), PA: ZINC15340668 complex (B) and PA: R05-2 complex (C) during MD simulation.

17

18

D

TE

EP

CC

A

SC RI PT

U

N

A

M

19

D

TE

EP

CC

A

SC RI PT

U

N

A

M

20

D

TE

EP

CC

A

SC RI PT

U

N

A

M

21

D

TE

EP

CC

A

SC RI PT

U

N

A

M

22

D

TE

EP

CC

A

SC RI PT

U

N

A

M

23

D

TE

EP

CC

A

SC RI PT

U

N

A

M

24

D

TE

EP

CC

A

SC RI PT

U

N

A

M