In Silico Design of New Inhibitors Against Hemagglutinin of Influenza

6 days ago - The RNA virus Influenza A is a serious public health problem with every year epidemics resulting in more than 250,000 deaths. A target zo...
0 downloads 0 Views 10MB Size
Subscriber access provided by University of Kansas Libraries

B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules

In Silico Design of New Inhibitors Against Hemagglutinin of Influenza Aurelie Perrier, Matthias Eluard, Michel Petitjean, and Anne Vanet J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b10767 • Publication Date (Web): 28 Dec 2018 Downloaded from http://pubs.acs.org on January 2, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

In Silico Design of New Inhibitors Against Hemagglutinin of Influenza Aurélie PERRIER,∗,†,‡,¶ Matthias ELUARD,†,§,¶ Michel PETITJEAN,k,¶ and Anne VANET∗,§,¶ †Chimie ParisTech, PSL Research University, CNRS, Institut de Recherche de Chimie Paris (IRCP), F-75005 Paris, France ‡Université Paris Diderot, Sorbonne Paris Cité, 5 rue Thomas Mann, F-75205 Paris Cedex 13, France ¶Epôle de Génoinformatique, Institut Jacques Monod, UMR7592, CNRS, F-75013 Paris, France §Pathologies de la réplication de l’ADN, Institut Jacques Monod, UMR7592, CNRS, F-75013 Paris, France kMTi, UMR-S 973, INSERM, University Denis Diderot, Paris 7, France E-mail: [email protected]; [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The RNA virus Influenza A is a serious public health problem with every year epidemics resulting in more than 250,000 deaths.

A protein cavity was identified

on the HA2 subunit of the hemagglutinin responsible for the entry of the virus into the host cell by endocytosis. The binding of a ligand in this zone rich in invariant residues and synthetic lethal couples, could prevent therapeutic escape and inhibit the conformational change at pH=5 initiating the membrane fusion in the endosome. Two pentapeptides, a linear peptide (EQRRS) and a cyclic peptide (DQRRD) have been proposed as potential ligands. Complex stability and the interactions between the ligand and the protein have been studied with the help of molecular dynamics and quantum chemistry methods. A high stability of the interactions has been obtained for these two ligands at both pH=7 and pH=5. Indeed, these two peptides present two cooperative modes of action that should prevent the conformational change at the origin of the spring-loaded mechanism at pH=5, (1) mechanical since they are docked on HA2, and (2) electronic since they modify the protonation states of key residues in the hairpin zone. This study thus paves the way towards the development of peptide ligands that can inhibit the membrane fusion process.

2

ACS Paragon Plus Environment

Page 2 of 34

Page 3 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Introduction Despite the efforts and success of antiviral research, the influenza virus continues to threaten public health because of its high rates of morbidity and mortality. Seasonal, it causes the death of 250,000 to 500,000 people. Four influenza pandemics were observed in the twentieth century: the Spanish flu (H1N1) of 1918 resulting in 20 to 40 million deaths, the Asian flu (H2N2) in 1957, Hong Kong influenza (H3N2) in 1968 and influenza A (H1N1) in 2009 (280,000 deaths). The search for new therapeutic molecules against the influenza virus is therefore a major issue and some molecules are currently under development or in clinical trials. 1 Potential targets for these molecules may be localised on host proteins but ligands that block their functions will lead to more rapid development of side effects. Thus, the most commonly sought targets are the functional proteins of the virus itself: entry inhibitors such as triterpenoids, neuraminidase inhibitors (NA) like the Zanamivir, 2 Oseltamivir, 3 Peramivir 4 and Laninamivir, 5 RNA-dependent RNA polymerase inhibitors, nucleocapsid protein inhibitors, and M2 ion channel inhibitors (Amantadine and Rimantadine 6 ). Despite the existence of these treatments, the influenza virus remains a major public health problem. Indeed, its high mutation rate allows it to quickly acquire resistance to these treatments. This is due to the RNA polymerase of influenza A virus, that have a high error rate (from 1 per 1000 to 1 per 100,000) and no error correction system. Thus, the treatments become ineffective in the long term (e.g. Amantadine and Rimantadine). In order to avoid a therapeutic escape, our team has developed a bioinformatic methodology 7–9 which consists in detecting invariant residues and synthetic lethals (SL) at the surface of viral proteins. 10 SLs represent pairs of residues for which either mutation is neutral but where the double mutation should make the virus unviable. Thus, these groups of residues constitute a potential therapeutic target, which should reduce therapeutic escape. This method allowed to define new targets on protease and reverse transcriptase of HIV 11,12 and on hemagglutinin of influenza. 13 The latter protein is a homotrimer (see Fig. 1) whose inactive monomer HA0 (566 amino acids, AA, with signal peptide) is cleaved in the extracellular medium 14,15 in 2 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

subunits HA1 (325 AA) and HA2 (222 AA). The HA1 subunit is structured as a globular head and has antigenic sites recognizing sialic acid on the surface of the cell to be infected. 16–18 The HA2 subunit is structured into an α helix stem and is strongly involved in the membrane fusion process thanks to the fusion peptide. 19–21

Figure 1: (A) 3D representation of the hemagglutinin H1 homotrimer (PDB code: 1RU7). Each monomer is composed of two polypeptides, HA1 (green) et HA2 (blue). (B) Structure of the HA1 polypeptide. The interaction zone with the sialic acid is represented in yellow. (C) Structure of the HA2 polypeptide. The fusion peptide which binds to the host cell membrane is part of the N-terminal end (in red) and the membrane anchor corresponds to the C-terminal end in blue When the virus enters into the endosome, the acidic pH (about 5) induces a conformational change: 22 the HA1 subunits deviate from each other and the HA2 subunit undergoes a conformational change called the spring loaded mechanism. This mechanism allows the hairpin of fusion peptide, to become linear to touch the endosome membrane, initiating the membrane fusion and the integration of the viral genome into the host cell. In a previous work, some of us found a group of SLs and invariant residues in this hairpin zone showing that this area would be a good potential target to avoid therapeutic escape. 13 This potential binding site, depicted on Fig. 2, is composed of nine residues. Two residues are SLs, the K395 and K398 lysines, while the residues E392, F393, N394, E397, E401, N402 and N405 are invariants. We have then checked if these identified functional residues correspond to a druggable pocket, that is to say that this target is able to accommodate small ligands. 23,24 4

ACS Paragon Plus Environment

Page 4 of 34

Page 5 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Among several softwares testing the druggability of potential cavities, 25,26 we have relied on the Fpocket code. 27 With a volume pocket of 221 Å3 , this protein cavity has been identified as a binding candidate pocket (four other binding sites were identified). In the same study, we proposed a tripeptide KFE consisting of two charged AAs (at a physiological pH) that can bind to the polar AAs of the target and a rather hydrophobic AA that could have affinity with the neighboring aromatic cycle of the protein residues.

Figure 2: Definition of the protein cavity (in green). The residues K395 and K398 are SLs (blue) and the residues E392, F393, N394, E397, E401, N402 et N405 are invariants (yellow). In the present study, we aim at proposing new ligands that could attach to this target and thus block the spring loaded mechanism essential to the pathogenicity of the virus. We decided to follow the strategy proposed in our previous study and focused our attention on the development of therapeutic peptides. Peptides have the advantage that their degradation products are non-toxic AA. 28 Compared with traditional drugs, peptides have the disadvantage of having a low capacity to cross membranes, as well as a short half-life (which can also be an advantage because the peptides do not accumulate in tissues). These drawbacks have stimulated the research in nanotechnologies with the final goal to produce lipidic and polymeric nanocarriers for peptide delivery. 29 Therapeutic peptides have already been 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

developed, for example against hepatitis C 30 or HIV. 31 Concerning the Influenza virus, peptides are also under development and we can cite for example, the entry block peptide, 32 FluPep, 33 or NDFRSKT peptide. 34 In this context, we have designed several polypeptides and tested their hairpin HA binding capacity with the help of docking simulations, molecular dynamics (MD) and calculations of the interaction energies within the MM/PBSA and Quantum Mechanics (QM) frameworks. Our study shows that two pentapeptides, one cyclic and one linear can attach to this target at pH=7 but also at pH=5 on both arms of the hairpin. Our simulations provide comprehensive insights on the interaction between the hemagglutinin and the polypeptide by assessing the contribution of the main amino acids to the binding energy. Additionally, a non-covalent interaction (NCI) analysis is presented to determine the nature of the interaction between the ligand and the binding pocket. Our theoretical study can thus provide useful insights to develop antiviral potential against the flu.

Methodology Hemagglutinin structures Over the 71 3D structures of hemagglutinin H1 available in the Influenza Research Database and in the Protein Data Bank (PDB), only 32 complete homotrimer structures were preserved (Table SI). The primary sequence alignments were performed with Clustal Omega 35 to observe the variations of the AA sequence in the area of interest. In order to choose the structures that will be used for the research of small therapeutic molecules, a clustering of these structures was carried out with the Modified Clustering Gain tool (MCG) 36 using the ascending hierarchical classification method. The classification criterion is the pairwise Root Mean Square Deviation (RMSD) computed with the MM-align algorithm 37 that computes the RMSD of two multimeric structures based on a sequence alignment. The representative structures (closest to the average structure) of each cluster are determined by the MCG tool. 6

ACS Paragon Plus Environment

Page 6 of 34

Page 7 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

As will be discussed later, the 1RU7 structure represents the largest cluster and will be used in the following computational steps.

Docking Simulations The structure of the 1RU7 protein used in the course of the docking simulation is determined from a clustering of 40 conformations extracted from the last 20 ns of the molecular dynamics of the isolated protein at pH = 7 (total simulation time: 50 ns, see next section). Clustering is performed with the MCG tool and is based on the pairwise RMSD calculated with ARMS 38 (adapted for calculating the RMSD between two conformations of the same protein). The conformation closest to the average structure of the largest cluster is chosen for docking. The interest of using a snapshot extracted from a MD simulation rather than a structure determined from X-ray diffraction experiments is to take into account the relaxation of the protein in the solvent, at T=310 K. For linear peptide ligands, minimized PDB structures are obtained with the PEP-FOLD3 online tool, 39 from the primary sequence of AAs. For the other cyclic peptides, since we could not use this tool, the structures were minimized at the Hartree Fock (HF)/6-31G level with Gaussian09 software. 40 We rely on this strategy because this QM calculation is able to provide both a starting structure for the docking simulation and the RESP charges for the force field setting (see next section). The files are prepared with Chimera’s DockPrep tool 41 in combination with Autodock Tools. 42 Flexible docking is carried out with AutoDock Vina 43 by defining the ligands and the side chains of the 9 residues of interest (E392, F393, N394, K395, E397, K398, E401, N402 and N405) as flexible (see Figure 2). We recall here that these nine residues correspond to the invariants and synthetic lethal couples identified in this protein cavity. 13 Taking into account the flexibility of these residues is a compromise between the rigid-body docking and the flexible model, that is to say between the exploration of the conformational space and the computational burden. 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

The grid box for ligand docking is defined around the area of interest with a dimension of 20 Å× 18 Å× 21 Å, i.e. 7560 Å3 with a spacing of 1 Å. For each ligand, twenty binding modes were generated with random starting positions for the ligand. For the linear and cyclic peptides under investigation, as explained in section "Ligand design strategy", the ligands were designed to obtain five specific (or target) interactions between a given amino acid (AAi ) of the ligand and a residue (Resi ) of the protein (with i = 1 to 5). To do so, we calculated the distances between all the atoms of Resi and AAi . If a distance between one atom of Resi and one of AAi is less than a given threshold (4 Å), then we consider that there is a specific interaction between the ligand AA and the protein residue. The number of "target" interactions (NTI) or contacts of interest for one specific ligand thus ranges between 0 and 5. The docking poses with a number of contacts of interest larger than a defined threshold (here, 3) were used as initial configurations for molecular dynamics simulations. We want to stress out that in this work, the NTI are computed by considering a unique criterion, namely the distance threshold. To refine our model, this distance criterion should be completed with the measure of the bond angle in the case of target H-bond interactions.

Molecular Dynamics Protocol Molecular Dynamics (MD) simulations were performed on the protein structure, in the absence and then in the presence of a ligand at pH=7 and 5. The simulations were performed using the Particle Mesh Ewald Molecular Dynamics program (PMEMD) 44,45 of the AMBER 14 code. 46 Force fields FF99SB-ILDN 47 and GAFF (General AMBER Force Field) 48 have been used for the protein and ligands, respectively. For the ligands, the RESP (Restrained Electrostatic Potential) atomic partial charges, calculated at the HF/6-31G level with Gaussian09, are used. The system was inserted into a periodic rectangular box of dimensions 98 Å× 102 Å× 153 Å, in the presence of an explicit solvent (water model TIP3P). 49 The systems were

8

ACS Paragon Plus Environment

Page 9 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

neutralized (addition of Na+ or Cl− counterions) in order to use the Ewald Particle Mesh method (PME) with a cut-off of 8 Å. MD were carried out after two minimization steps (the solvent alone and then the complete system). The minimized system was slowly heated to 310 K for 600 ps in the NVT ensemble. During the dynamics, the temperature was kept constant at 310 K with the help of the Berendsen thermostat. 50 NVT was followed by an equlibration step of 100 ps in the NPT ensemble. The production phase was also carried out in the NPT ensemble with the help of the Berendsen barostat (P =1 atm). The time constant is 2 fs (we used the SHAKE algorithm), the MD simulations were performed for 50 ns and coordinates were saved every 2 ps. The MD were carried out with the protonation states of the residues corresponding to pH = 7 and then to pH=5. Without ligand, we used the protonation state defined by previous studies. 51,52 At pH = 7, all the glutamic acids in the area of interest are negatively charged. The acidic pH of the endosome ( pH ' 5) changes the protonation state of the Glu, Asp and His residues. 22,51 These works have shown that the glutamic acids E392 and E397 (Fig. 2) are protonated at pH=5 while the residue E401 keeps its negative charge. In presence of a ligand, we have used the protonation states of the charged residues determined with the help of PROPKA 3.0. 53,54 More precisely, we have evaluated the effect of the ligands on the pKa values of the Glu residues inside the area of interest. We have extracted 40 conformations generated during the 20 last nanoseconds of the MD simulations at pH=7 for the protein 1RU7 with and without ligand. These 40 conformations were analyzed with PROPKA 3.0 through the on-line server PDB2PQR 55 in order to estimate the pKa of the residues at pH=5. A pKa value smaller than 5 indicates that the Glu residue is protonated. Hemagglutinin is an homotrimer, thus we analyzed 9 residues (3 per monomer): E392, E397. E401. E875. E880. E884, E1358, E1363 and E1367. Each residue is analyzed either in presence of the ligand or not, so that we analyzed a total of 18 populations of 40 pKa values. By using several statistical tests (mean pairwise comparison with Wilcoxon-

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Mann-Whitney test and one-sample Wilcoxon test) detailed in the Supporting Information section for the EQRRS ligand, 56,57 we have concluded that the pKa values of the Glu residues are lowered in presence of a ligand: the carboxy group of the Glu residues are charged. The results of molecular dynamics simulations are visualized with VMD (Visual Molecular Dynamics). 58 We were interested in the calculation of the Root Mean Square Deviation (RMSD) and the Root Mean Square Fluctuation (RMSF). The latter is the measure of the fluctuation of the atomic or residue position with respect to its average position, plotted as a function of the simulation time. The RMSF provides information on the flexible regions of a molecule.

Analysis of the interactions between the ligand and the protein To compute the protein-ligand binding free energy, we have first used the combined Molecular Mechanics/Generalized Born Surface Area method (MM / GBSA). 59 This approach, which offers an affordable methodology for predicting the binding energies in large biosystems, 60,61 is based on the difference between the free energies of the protein, the ligand and the complex in solution. At that stage, we have considered 400 conformations extracted from the last 20 ns of the simulation. Then, "ligand-protein" structures were extracted from the molecular dynamics trajectories (100 conformations for each trajectory, extracted every 200 ps during the last 20 ns of the dynamics). The pairwise interaction energies, ∆Epair , between each residue of interest (Res) of the protein and the ligand (Lig) were calculated using the following formula: ∆Epair = ERes−Lig −ERes −ELig . 62 The energies of the ERes−Lig Residue-Ligand complex, the residue of interest ERes and the ligand ELig are obtained with the help of Quantum Mechanics (QM) calculations relying on the Density Functional Theory (DFT). The functional M062X 63 which leads to a very good description of the non-covalent interactions was used in combination with the 6-311++G(d, p) basis set in presence of an implicit solvent described by the so-called Polarizable Continuum Model (PCM). 64 Finally, for some conformations, 10

ACS Paragon Plus Environment

Page 10 of 34

Page 11 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

an analysis of the non-covalent interactions was carried out with NCIplot code. 65,66 In the NCI model, non-covalent interactions are characterized by colored isosurfaces obtained by combining the calculated electron density and its reduced gradient.

Results and discussion Properties of the H1 hemagglutinin. The clustering of the 32 H1 hemagglutinin structures has yielded 4 clusters: 15 structures represented by 1RU7, 12 structures represented by 3UBN, 4 structures represented by 5UJZ and the 4M4Y structure alone. The comparison of the sequences is discussed in Supporting Information. The 1RU7 structure (of the H1N1 strain of 1934) represents the largest cluster and will be used in this study. The average RMSD calculated for the 32 structures of H1 is 1.34 Å, which shows the homology of the 3D structures of H1 hemagglutins. We have thus carried out MD simulations on the 1RU7 structure at pH=7 and pH=5 without ligands and we have focused our attention on the structural modifications of the protein pocket represented on Fig. 2. The calculation of the RMSD (Fig. S2) shows that the 1RU7 protein reaches an equilibrium state after 30 ns. The last 20 ns of the dynamics can therefore provide representative conformations of 1RU7 that will be used for ligand docking. We obtained the same conclusions for the MD simulation at pH=5. In the course of the first 20 ns, the protein structure is particularly destabilized by the modification of the protonation states, in particular those of the glutamic and aspartic acids E397 and E392 in the protein cavity, and by the appearance of new repulsive interactions with the protonation of some histidines. The simulation time is too short to visualize the conformational change described in the introduction. Fig. 3(A) shows the RMSF obtained at pH=7 and pH=5 for one monomer and a zoom of the protein cavity (defined in Fig. 2) is provided on Fig. 3(B). Theprotein cavity is not very flexible in the course of the MD simulation, the RMSF value being systematically 11

ACS Paragon Plus Environment

The Journal of Physical Chemistry

HA1

HA2

(B)

RMSF (Å)

(A)

RMSF (Å)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

Residue

Residue

Figure 3: (A) RMSF calculated for the Cα of the AAs of the protein 1RU7 at pH=7 (blue) and pH=5 (red). Only one monomer of Fig. 1 is represented. The HA1 and HA2 subsunits are respectively identified with orange and turquoise lines. (B) Zoom: RMSF calculated for the residues of interest in the protein pocket. smaller than 2.0 Å. As shown by Figure 3(A), this region corresponds to one of the most stable regions of the monomer. A slightly higher flexibility (RMSF =1.7 Å) is observed for the Lys395 whose position corresponds to the bend of the loop. The modification of the pH does not lead to significant variations in the flexibility. The analysis of the RMSF for the complete homotrimer (see Fig. S3) comforts these conclusions.

Ligand design strategy and docking results As depicted on Fig. 4, we have selected five residues in the therapeutic protein pocket previously described (Fig. 2) that are distributed on both sides of the peptide fusion loop. The design strategy to select these residues relies of the formation of electrostatic interactions between the ligand and the K395, E392, E401 and K398 residues. Besides, the ligand should interact with the N394 residue through the formation of a hydrogen bond. The aim is to (i) preserve a sufficient number of contacts along the MD simulation in order to maintain the stability of the complex, (ii) target a couple of synthetic lethals (K395 and K398) and three

12

ACS Paragon Plus Environment

Page 13 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

invariants (N394, E392 and E401), and (iii) prevent the change of conformation at acidic pH by binding with residues localized on both sides of the peptide fusion loop.

Figure 4: Representation of the five targeted residues (K395, N394, E392, E401 et K398) in the zone of interest. The residues K395 and K398 are SLs (blue) while N394, E392 and E401 are invariants (yellow). In this framework, 16 linear pentapeptides have been proposed. In order to establish the desired contacts with the five residues of interest, the peptides possess the following sequence of five AAs: (1) Asp or Glu to create an electrostatic interaction with K395; (2) Asn or Gln to form a hydrogen bond with N394; (3) and (4) His or Arg to provide interactions with E392 and E401 (that are deprotonated at pH=7), either with a hydrogen bond with the short side chain of His, or through electrostatic interaction with the longer side chain of Arg; (5) Ser (short chain) to create an interaction between K398 and the positive charge in C-ter. Compared to our previous study where a tripeptide was assessed, 13 we directly study pentapeptides. The advantage of using larger peptides is to target five different key residues which constitute as many strong anchors, the inconvenient is that the control of the larger number of degrees of freedom compared to middle size peptides (like tetrapeptides) can be more difficult (larger entropy contributions).

13

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The structures of the 16 possible peptide combinations were docked on 1RU7 and the results are presented in Table 1. To analyze the docking results, we have focused our attention on two criteria: (1) the ligand should present a large binding affinity with the HA (which is measured by the docking score) and (2) it should target the key residues we are interested in. This second criterion can be evaluated thanks to the calculation of the number of target interactions (NTI) defined in the methodology section. Concerning the first criterion, Table 1 shows that the differences between the computed average docking scores are smaller than the standard error of the Vina docking scores (2.8 kcal.mol−1 ). Indeed, the average docking score ranges between -7.03 and -5.09 kcal.mol−1 for the 16 peptides. It is thus impossible to discriminate the ligands on the basis of this first criterion. One can also note that for each ligand, the standard deviation of the docking scores is small (0.20 kcal.mol−1 on average) which prevents the discrimination of one docking pose among the others. For each peptide, we thus focused our attention on the second criterion, namely the NTI for the "best" configuration (i.e. the docking pose with the largest NTI). We also consider for each ligand the average number of NTI calculated from the 20 poses. As shown by Table 1, the linear peptides with the AA pair RR or HR in positions 3 and 4 provide a higher NTI. For these peptides, some docking poses present more than 3 target contacts over the 5 desired interactions. In particular, the docking of the peptide EQRRS revealed a conformation with 4 target contacts. The more flexible side chain of the arginine allows this AA to adapt to the movement of the residues E392 and E401. Second, the peptides that can interact more favorably with N394 present a glutamine at the 2-position, the contact being facilitated by the longer side chain of the glutamine compared to the asparagine one. Finally, over the 16 peptides, we observe no docking poses with the 5 desired contacts. In order to reduce the flexibility of the peptides and strengthen the interactions, the linear peptides are cyclized through the formation of a covalent bond between N-ter and Cter. Relying on the identification of the AA combinations providing stronger interactions for

14

ACS Paragon Plus Environment

Page 14 of 34

Page 15 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Table 1: Results of fhe flexible molecular docking for the 16 linear peptides and 8 cyclic peptides. The sequence of the peptide, the type of peptide, the average docking score on the 20 poses (Av. Score, in kcal.mol−1 ), the standard deviation of the score (Score Sd.), the average number of target interactions (Av. NTI) over the 20 poses (for each pose, the number of interaction is between 0 and 5) and the number of target interactions for the "best" docking pose (NTI) are given. Sequence Type Av. Score Score sd. Av. NTI NTI DNRRS linear -5.59 0.21 1.10 3 DNHRS linear -5.56 0.23 1.05 3 EQRRS linear -5.49 0.18 1.05 4 ENRRS linear -6.40 0.19 1.00 3 DQRRS linear -5.34 0.18 0.90 3 DNRHS linear -5.16 0.24 0.80 1 EQRHS linear -5.98 0.13 0.75 2 ENHRS linear -5.86 0.20 0.70 3 ENHHS linear -5.67 0.13 0.70 3 EQHRS linear -5.41 0.20 0.70 3 ENRHS linear -6.07 0.13 0.65 1 DQHRS linear -5.58 0.20 0.65 3 DNHHS linear -5.48 0.17 0.65 2 DQHHS linear -5.46 0.24 0.60 2 EQHHS linear -5.82 0.21 0.45 1 DQRHS linear -5.89 0.19 0.40 2 DQRRE cyclic -6.43 0.24 1.20 3 EQRRE cyclic -6.28 0.18 1.20 4 DQHRE cyclic -5.32 0.21 0.90 3 EQHRD cyclic -6.94 0.23 0.75 3 DQRRD cyclic -7.03 0.26 0.65 3 EQRRD cyclic -5.09 0.20 0.65 2 EQHRE cyclic -5.53 0.13 0.60 3 DQHRD cyclic -5.87 0.32 0.45 2 the linear peptides, we have proposed 8 cyclic pentapeptides listed in Table 1 and the average docking score for these ligands is more favorable than the linear peptides ones. However, the significance of this increase could not be demonstrated due to the standard error of the Vina docking scores. For six cyclic peptides, we obtain docking poses with at least 3 out of 5 contacts (4 contacts for EQRRE).

15

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

According to these criteria, 6 peptides were retained: the linear EQRRS represented on Fig. 5(A) and the cyclic EQRRE (with 4 contacts out of 5), and 4 cyclic peptides DQRRE, DQRRD depicted on Fig. 5(B), EQHRD and DQHRE (with 3 contacts out of 5). MD simulations were performed at pH=7 for these 6 peptides, the starting configuration corresponding to the selected pose. Only the linear EQRRS and cyclic DQRRD ligands were able to form a stable complex with the 1RU7 protein in the course of the MD simulation and hereafter, we will focus our discussion on these two pentapeptides.

Figure 5: (A)docking pose selected for the linear peptide EQRRS in interaction with 1RU7. Four out of the five target contacts are obtained (with K395, E392, N394 and E401). (B) docking pose selected for the cyclic peptide DQRRD with 1RU7. Three out of the five target contacts (with E392, N394 and E401) are obtained.

Interaction of the EQRRS linear peptide with 1RU7 From the docking pose provided in Fig. 5(A), a MD simulation was carried out at pH=7. To follow the evolution of these interactions, the distances characterizing the target contacts were measured over the 50 ns (see Fig. 6). These results show that the interactions between the ligand and the residues N394 (Interaction 2), E392 (Interaction 3) and E401 (Interaction 4) are stable with contact distances smaller than 5 Å. The interaction between K398 and the C-ter of the peptide (Interaction 5) is also stable with an average distance of 3 Å, the 6 Å steps corresponding to a rotation of the C-N-Cα -Cb bond of the serine. A gradual stall of the contact 1 (between the Glu AA, EQRRS, and K395) is observed with a distance ranging 16

ACS Paragon Plus Environment

Page 16 of 34

Page 17 of 34

from 7 to 11 Å. For the interaction 2, there is a significant change in the contact distance at around 48-49 ns. Before 48 ns and after 49 ns, the closest contact corresponds to the interaction between the oxygen atom of the Q lateral chain and one of the H atoms of the N394 NH2 terminal group. Around 48-49 ns, there is a conformational change and the closest contact corresponds to a H-bond interaction between the amino groups of Q and Asn. Due to steric constraints, this conformation is not stable.

Distance

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 0

Interaction 1 Interaction 2 Interaction 3 Interaction 4 Interaction 5

5

10

15

20

25

30

35

40

45

50

Time (ns) Figure 6: Distances (in Å) characterizing the 5 target conctacts between one AA of the EQRRS ligand and one of the 5 residues of 1RU7 (pH=7). Interaction 1: E-K395, Interaction 2: Q-N394, Interaction 3: R-E392, Interaction 4: R-E401, Interaction 5: S-K398. According to the pairwise energies ∆Epair depicted in Fig. 7, the EQRRS ligand is mainly coordinated to the E392 residue with an average interaction energy of ∼-20.0 kcal.mol−1 over the last 20 ns of the MD simulation. The pairwise energy of the ligand with E401 (respectively N394) is constant with a value of ∼-9.0 kcal.mol−1 (resp. ∼-15.0 kcal.mol−1 ) while a small 17

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

decrease of the ∆Epair absolute value is observed for the K395 residue. This destabilization can be directly linked to the increase of the distance between the ligand and this residue. Finally, the rotation of the C-N-Cα -Cb bond of the ligand serine observed in Fig. 6 leads to a decrease of the pairwise energy of the EQRRS peptide with the K398 residue.

NCI

Figure 7: Evolution of the QM Pairwise energies ∆Epair in kcal.mol−1 of the EQRRS peptide with the five residues of interest (K395: black, N394: red, E392: green, E401: blue, K398: orange) along the MD trajectory. The structure extracted at t=45.5 ns has been analyzed with NCIplot (see Fig. 8). The NCI analysis of the conformation extracted at 45.5 ns depicted in Fig. 7 provides a complementary qualitative picture of the ligand binding mode with the five targeted residues. The H-bonds are displayed as blue while weak attractive interactions are represented as green regions. Indeed, three out of the five target contacts correspond to the expected interactions, namely S with K398 (interaction 5), R (EQRRS) with E401 (interaction 4) and R (EQRRS) with E392 (interaction 3). More precisely, hydrogen bonds are present 18

ACS Paragon Plus Environment

Page 19 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

between (i) the COO− terminal group of the serine and the residue K398, (ii) an amine of the second arginine EQRRS of the peptide and E401 as well as (iii) two H-bonds between this same arginine and E392. An attractive electrostatic interaction is also present between the arginine EQRRS and the residue E392, which corresponds to a target interaction. The latter observations can thus rationalize the most favorable calculated pairwise energy ∆Epair of the ligand with the residue E392. In addition, a H-bond is present between the NH of the backbone of the N394 residue and a CO group of the peptide. The desired H-bond between the two side chains of the Q AA and the N394 residue is not present because the ends are too far apart (4.8 Å). Finally, there is no strong interaction between the glutamic acid (E) charges of the peptide and K395, only a weak Van der Waals interaction between the aliphatic chains is present. Therefore, despite the flexibility of the linear peptides, this MD simulation shows a high stability of the complex at pH=7 due to the formation of 6 H-bonds. This is confirmed by the calculation of the interaction energy with the MM / GBSA method. The calculated binding energy, -27.9 kcal.mol−1 (see Table SII), proves that the EQRSS-1RU7 complex is stable. For information, this assembly is much more stable than the KFE-1RU7 complex proposed in a previous study. 13 Using the same computational MD protocol for the latter system, we found a MM /GBSA binding free energy of -3.6 kcal.mol−1 . Since hemagglutinin is a homotrimer, we considered the possibility to bind three different EQRSS ligands with the 1RU7 protein. Namely, each ligand was positioned in the same region of interest depicted in Fig. 2 and a MD simulation of 50 ns was conducted. The variations of the RMSD as well as the distances characterizing the target contacts for each couple of ligand / monomer are respectively given in Figures S6 and S7. These results show that the three ligands stay in the defined binding region during the MD simulation. For all the ligand-monomer complexes, the shortest distance between EQRSS and the protein corresponds to the interaction 3 (R-E392) and the largest one to the interaction 1 (E-K395).

19

ACS Paragon Plus Environment

K-395

E

Q

!H-bond "

S

N-394

R

H-bond

K-398 H-bond

R

H-bonds E-392

E-401

Figure 8: NCI analysis of the non-covalent interactions between the ligand EQRRS and the five residues of interests (K395: black box, N394: red, E392: green, E401: blue, K398: orange). Isosurfaces correspond to s =0.5 a.u. With this more realistic model, this simulation thus confirms the ability of the EQRRS ligand to bind to the hemagglutinin protein. Finally, we have investigated the stability of the EQRRS-1RU7 complex at pH=5. To this purpose, we have first evaluated the effect of the EQRRS ligand on the pKa values of the Glu residues E392, E397 and E401 inside the area of interest, following the procedure detailed in the Methodology section. We have shown that there is an effect of the peptide binding on the pKa of the Glu residues interacting with this ligand (see Supporting Information). More precisely, the pKa values of E392 and E491 decrease in presence of the ligand and these two 20

Page 21 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

residues are thus charged at pH=5. We remind here that the Glu E392 was protonated for 1RU7 (without ligand) at acidic pH. In this framework, at pH=5, the calculation of the RMSD shows that the complex reaches an equilibrium state after 30 ns (see Fig S4). Besides, the analysis of the distances characterizing the target contacts (Fig. S8) reveals that the EQRSS-1RU7 complex is stable at pH=5. Over the MD run, the closest distances correspond to the interactions 5 (E-K398),4 (R-E401) and 3 (R-E392) as illustrated on a representative snapshot on Fig 9. Therefore, there is no decomplexation of the EQRSS peptide at both pH=7 and pH=5. Overall, three target contacts (with K398, E401 and E392) are maintained in the course of the MD simulations. Interestingly, these three contacts are on both sides of the loop which may inhibit the conformational change at the origin of the membrane fusion process. Moreover, the modification of the pKa values and charges of the Glu residues due to the presence of the ligand should also impact the conformational change and the overall spring-loaded mechanism. The EQRSS ligand thus presents two cooperative modes of action, through steric and electronic effects.

Interaction of the DQRRD cyclic peptide with 1RU7 We then considered the interaction of the DQRRD cyclic peptide with 1RU7 and we started the MD simulation with the docking pose, represented in Fig. 5, that shows the presence of 3 out of the five 5 target contacts (with E392, N394 and E401). At pH=7, the distances characterizing the target contacts (Fig S.9) shows that the ligand is in contact with the 3 residues: E392 (stable distance of ∼ 2Å with some rare fluctuations at 7Å) K395 and K398 (very stable interaction at ∼ 3Å). Nevertheless, there is no contact with the residues N394 (at ∼ 12Å) and E401 (at ∼ 9Å). The visualization of the trajectory reveals that the side chains of the residues Q and R of the ligand are folded over the cycle, not allowing a contact with N394 and E401.

21

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Representative snapshot of the MD simulation of the EQRRS interacting with 1RU7 at pH=5. The calculation of the QM pairwise energies ∆Epair (see Fig. 10) shows an equivalent contribution for the 3 contacts with E392, K395 and K398, with a greater energy stability for K398 (at ∼-12 kcal.mol−1 ) and some energy fluctuations for E392 (between -20 and - 7 kcal.mol−1 ), which is correlated with the distance variations previously discussed. Moreover, the ∆Epair energies calculated for the conformations extracted between 34.5 ns and 36.5 ns show the existence of some interactions with the residues N394 (-5 kcal.mol−1 ) and E401 (-16 kcal.mol−1 ). The NCI analysis of the conformation extracted at t=41.5 ns shows that the stability of the 3 interactions with E392, K398 and K395 arises from the formation of 3 H-bonds with these residues. In addition, 3 additional interactions with residues that are outside the area of interest contribute to the stability of the complex: (1) an H bond between the Q587 residue of the HA1 subunit and the COO− group of a D amino acid (DQRRD) close to the contact 5 with K398, (2) a H-bond between the K717 residue of the HA1 subunit and the

22

ACS Paragon Plus Environment

Page 22 of 34

Page 23 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

NCI

Figure 10: Evolution of the QM Pairwise energies ∆Epair in kcal.mol−1 of the DQRRD peptide. See Figure 7 for more details. The structure extracted at t=41.5 ns has been analyzed with NCI (see Fig. 11). COO− group of the same D amino acid (close to K395) and (3) a H-bond between the E654 residue of the HA1 subunit and the NH2 group of the second arginine (DQRRD), close to the contact 4. Finally, the stability of the DQRRD-1RU7 complex was investigated at pH=5 with the help of MD simulations. At that stage, due to the presence of the ligand, the two E392 and E491 residues are deprotonated. The analysis of distances characterizing the target contacts (Fig S10) reveals that during the first 38 ns of the simulation, the ligand mainly interacts with E392 (stable distance of ∼ 2Å with some rare fluctuations) and with the residues K395 and K398 with an average distance of ∼ 4Å and larger fluctuations. Contrary to the simulation at pH=7, the side chains of the DQRRD residues are no longer folded over the

23

ACS Paragon Plus Environment

Q587 K395 H-bond H-bond

C

K398

K717

E392 H-bond

H-bond

E654 Figure 11: NCI analysis of the non-covalent interactions between the ligand DQRRD and the five residues of interests. See Fig. 8 for more details. cycle, allowing an interaction with E392. After 38 ns, there is a conformational change: the ligand remains in the binding pocket and presents interactions with the residues K395 and E392 and, to a lesser extent, with N394 and K398. Therefore, the DQRRD-1RU7 complex is stable over the MD simulations at pH=7 and pH=5 thanks to the existence of stable contacts with the target residues. Interestingly, these different contacts are systematically on both sides of the loop and always involve a couple of SLs (K398 and K395).

24

Page 25 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Conclusion The aim of this study was to propose new ligands potentially active against the Influenza A virus. To this purpose, we have relied on molecular docking, molecular dynamics simulation at pH=7 and 5, and the analysis of the interaction energies with MM / GBSA and QM approaches. We have defined five target residues in the HA2 subunit of the hemagglutinin, with 3 invariants and a couple of synthetic lethals, and proposed linear and cyclic pentapeptides ligands with possible affinity towards the area of interest. Molecular docking has enabled us to select two ligands, the linear peptide EQRRS and the cyclic peptide DQRRD. The MD simulations have then revealed the high stability of the ligand-protein complex for these two peptides. At pH=7 and 5, for the EQRRS peptide, three stable and strong interactions contacts are maintained in the course of the MD runs (with E392, E401 and K398). For the DQRRD peptide, a minimum of three out of the five target contacts are conserved along the simulations. These 2 complexes are stabilized by the formation of H-bonds that strenghten the electrostatic interactions, as shown by the NCI analysis. Importantly, for EQRRS, the 3 residues that are interacting with the ligand are localized on both sides of the peptide fusion loop, which may inhibit the conformational change at the origin of the membrane fusion process at pH=5. In the same vein, the DQRRD peptide is systematically bound to a couple of synthetic lethals which is of high interest to make the virus unviable. It is also important to note that these two ligands have an influence on the pKa of the protein residues, particularly concerning the glutamic acid E392 one. The binding of the ligand complicates the protonation of the E392 residue and this phenomenon should maintain the stability of the complex at pH=5. These two peptides thus present two cooperative modes of action that should prevent the spring-loaded mechanism at pH=5, (1) mechanical since they are docked on HA2, and (2) electronic since they modify the protonation states of key residues in the hairpin zone. Hence, they represent valuable candidates to constitute new therapeutic molecules against the influenza virus.

25

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

As a perspective, the study of these 2 peptides will also be carried out with the 3UBN PDB structure (the representative structure of the second most important cluster for the HA proteins). In the target region, compared to 1RU7, 3UBN presents a a K395H mutation and the impact of this mutation on the stability of the complex at pH=7 and pH=5 will be investigated. Besides, simulations on larger time scale in presence of the endosome membrane should be considered to confirm that the spring-loaded mechanism is efficiently blocked by the ligand. Finally, in vitro tests are planned to experimentally validate the efficacy of these peptides in blocking the entry of the virus into the cell.

Supporting Information Available H1 structures. Analysis of the MD simulation of 1RU7 without ligand. Complementary analysis of the interaction of the EQRRS linear peptide with 1RU7. Effect of the ligand on the pKa values: results of the statistical tests. Complementary analysis of the interaction of the DQRRD cyclic peptide with 1RU7.

Acknowledgements This work was granted access to the HPC resources of CINES and IDRIS under the allocation 2017-A0010810135 made by GENCI (Grand Equipement National de Calcul Intensif). We thank the experts of the genoinformatics epole as well as G. Baldacci who allowed its creation. F. Fauchereau and L. Le Bras are acknowledged for helpful scientific discussions. We thank P. Le Chien for his constant encouragements since 1993. This work was done in memory of Lou, who died of AIDS at the age of 5.

26

ACS Paragon Plus Environment

Page 26 of 34

Page 27 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

References (1) Wu, X.; Wu, X.; Sun, Q.; Zhang, C.; Yang, S.; Li, L.; Jia, Z. Progress of small molecular inhibitors in the development of anti-influenza virus agents. Theranostics 2017, 7, 826–845. (2) Holzer, C. T.; Von Itzstein, M.; Jin, B.; Pegg, M. S.; Stewart, W. P.; Wu, W.-Y. Inhibition of sialidases from viral, bacterial and mammalian sources by analogues of 2-deoxy-2,3-didehydroN-acetylneuraminic acid modified at the C-4 position. Glycoconjugate Journal 1993, 10, 40–44. (3) Kim, C. U.; Lew, W.; Williams, M. A.; Liu, H.; Zhang, L.; Swaminathan, S.; Bischofberger, N.; Chen, M. S.; Mendel, D. B.; Tai, C. Y. et al. Influenza neuraminidase inhibitors possessing a novel hydrophobic interaction in the enzyme active site:? design, synthesis, and structural analysis of carbocyclic sialic acid analogues with potent anti-influenza activity. J. Am. Chem. Soc. 1997, 119, 681–690. (4) Kohno, S.; Yen, M.-Y.; Cheong, H.-J.; Hirotsu, N.; Ishida, T.; Kadota, J.; Shimada, J. Phase III randomized, double-blind study comparing single-dose intravenous peramivir with oral oseltamivir in patients with seasonal influenza virus infection. Antimicrobial Agents and Chemotherapy 2011, 55, 5267?–5276. (5) Yamashita, M.; Tomozawa, T.; Kakuta, M.; Tokumitsu, A.; Nasu, H.; Kubo, S. CS-8958, a prodrug of the new neuraminidase inhibitor R-125489, shows long-acting anti-Influenza virus activity. Antimicrobial Agents and Chemotherapy 2009, 53, 186–192. (6) Pielak, R. M.; Chou, J. J. Influenza M2 proton channels. Biochim. Biophys. Acta 2011, 1808, 522 – 529. (7) Vanet, A.; Muller-Trutwin, M.; Valere, T. Method for identifying motifs and/or combinations of motifs having a boolean state of predetermined mutation in a set of sequences and its application. 2010; US7734421B2. (8) Vanet, A.; Muller-Trutwin, M.; Valere, T. Method for identifying motifs and/or combinations of motifs having a boolean state of predetermined mutation in a set of sequences and its application. 2011; US8032309B2.

27

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(9) Vanet, A.; Muller-Trutwin, M.; Valere, T.; Brouillet, S.; Ollivier, E.; Marsan, L. Method for identifying combinations of motifs that do not mutate simultaneously in a set of viral polypeptide sequences comprising a putative drug binding site. 2011; US7917303B2. (10) Bazin, C.; Coupage, R.; Middendorp, S.; Vanet, A. Between compensatory mutations and synthetic lethals: genetic mutations, a new challenge for tomorrow’s medicine. Science postprint 2014, 1, e00035. (11) Brouillet, S.; Valere, T.; Ollivier, E.; Marsan, L.; Vanet, A. Co-lethality studied as an asset against viral drug escape: the HIV protease case. Biology Direct 2010, 5, 40. (12) Petitjean, M.; Badel, A.; Veitia, R. A.; Vanet, A. Synthetic lethals in HIV: ways to avoid drug resistance. Biology Direct 2015, 10, 17. (13) Lao, J.; Vanet, A. A new strategy to reduce influenza escape: detecting therapeutic targets constituted of invariance groups. Viruses 2017, 9 . (14) Clancy, S. Genetics of the influenza virus. Nature Education 2008, 1, 83. (15) Steinhauer, D. A. Role of hemagglutinin cleavage for the pathogenicity of influenza virus. Virology 1999, 258, 1 – 20. (16) Schmidt, A. G.; Therkelsen, M. D.; Stewart, S.; Kepler, T. B.; Liao, H.-X.; Moody, M. A.; Haynes, B. F.; Harrison, S. C. Viral receptor-binding site antibodies with diverse germline origins. Cell 2015, 161, 1026 – 1034. (17) Gamblin, S. J.; Skehel, J. J. Influenza hemagglutinin and neuraminidase membrane glycoproteins. J. Biol. Chem. 2010, 285, 28403–28409. (18) Gamblin, S. J.; Haire, L. F.; Russell, R. J.; Stevens, D. J.; Xiao, B.; Ha, Y.; Vasisht, N.; Steinhauer, D. A.; Daniels, R. S.; Elliot, A. et al. The structure and receptor binding properties of the 1918 influenza hemagglutinin. Science 2004, 303, 1838–1842. (19) Carr, C. M.; Chaudhry, C.; Kim, P. S. Influenza hemagglutinin is spring-loaded by a metastable native conformation. Proc. Natl. Acad. Sci. USA 1997, 94, 14306–14313.

28

ACS Paragon Plus Environment

Page 28 of 34

Page 29 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(20) Skehel, J. J.; Wiley, D. C. Receptor binding and membrane fusion in virus entry: the influenza hemagglutinin. Annual Review of Biochemistry 2000, 69, 531–569. (21) Garcia, N. K.; Lee, K. K. Dynamic viral glycoprotein machines: approaches for probing transient states that drive membrane fusion. Viruses 2016, 8 . (22) Xu, R.; Wilson, I. A. Structural characterization of an early fusion intermediate of influenza virus hemagglutinin. J. Virol. 2011, 85, 5172–5182. (23) Trosset, J.-Y.; Vodovar, N. In Target identification and validation in drug discovery: methods and protocols; Moll, J., Colombo, R., Eds.; Humana Press: Totowa, NJ, 2013; pp 141–164. (24) Mangiatordi, G. F.; Alberga, D.; Siragusa, L.; Goracci, L.; Lattanzi, G.; Nicolotti, O. Challenging AQP4 druggability for NMO-IgG antibody binding using molecular dynamics and molecular interaction fields. Biochim Biophys Acta Biomembr. 2015, 1848, 1462 – 1471. (25) Halgren, T. A. Identifying and characterizing binding sites and assessing druggability. J. Chem. Inf. Model. 2009, 49, 377–389. (26) Cruciani, G.; Pastor, M.; Guba, W. VolSurf: a new tool for the pharmacokinetic optimization of lead compounds. Eur. J. Pharm. Sci. 2000, 11, S29 – S39. (27) Le Guilloux, V.; Schmidtke, P.; Tuffery, P. Fpocket: An open source platform for ligand pocket detection. BMC Bioinformatics 2009, 10, 168. (28) Albert, L. Peptides as drugs: is there a market? J. Pept. Sci. 2002, 8, 1–7. (29) Santalices, I.; Gonella, A.; Torres, D.; Alonso, M. J. Advances on the formulation of proteins using nanotechnologies. Journal of Drug Delivery Science and Technology 2017, 42, 155 – 180. (30) Cheng, G.; Montero, A.; Gastaminza, P.; Whitten-Bauer, C.; Wieland, S. F.; Isogawa, M.; Fredericksen, B.; Selvarajah, S.; Gallay, P. A.; Ghadiri, M. R. et al. A virocidal amphipathic helical peptide that inhibits hepatitis C virus infection in vitro. Proc. Natl. Acad. Sci. USA 2008, 105, 3088–3093.

29

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(31) Cooper, D. A.; Lange, J. M. Peptide inhibitors of virus?cell fusion: enfuvirtide as a case study in clinical discovery and development. Lancet Infect. Dis. 2004, 4, 426–436. (32) Jones, J. C.; Turpin, E. A.; Bultmann, H.; Brandt, C. R.; Schultz-Cherry, S. Inhibition of influenza virus infection by a novel antiviral peptide that targets viral attachment to cells. J. Virol. 2012, 80, 11960–11967. (33) Nicol, M. Q.; Ligertwood, Y.; Bacon, M. N.; Dutia, B. M.; Nash, A. A. A novel family of peptides with potent activity against influenza A viruses. J. Gen. Virol. 2012, 93, 980–986. (34) Rajik, M.; Jahanshir, F.; Omar, A.; Ideris, A.; Hassan, S.; Yusoff, K. Identification and characterisation of a novel anti-viral peptide against avian influenza virus H9N2. Virol. J. 2009, 6, 74. (35) Sievers, F.; Wilm, A.; Dineen, D.; Gibson, T. J.; Karplus, K.; Li, W.; Lopez, R.; McWilliam, H.; Remmert, M.; Söding, J. et al. Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Molecular Systems Biology 2011, 7 . (36) Meslamani, J. E.; André, F.; Petitjean, M. Assessing the geometric diversity of cytochrome p450 ligand conformers by hierarchical clustering with a stop criterion. J. Chem. Inf. Model. 2009, 49, 330–337. (37) Mukherjee, S.; Zhang, Y. MM-align: a quick algorithm for aligning multiple-chain protein complex structures using iterative dynamic programming. Nucleic Acids Res. 2009, 37, e83. (38) Petitjean, M. On the root mean square quantitative chirality and quantitative symmetry measures. ?J. Math. Phys 1999, 40, 4587–4595. (39) Lamiable, A.; Thévenet, P.; Rey, J.; Vavrusa, M.; Derreumaux, P.; Tufféry, P. PEP-FOLD3: faster de novo structure prediction for linear peptides in solution and in complex. Nucleic Acids Res. 2016, 44, W449–W454. (40) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09 Revision D.01. Gaussian Inc. Wallingford CT 2009.

30

ACS Paragon Plus Environment

Page 30 of 34

Page 31 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(41) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera?A visualization system for exploratory research and analysis. J. Comput. Chem. 2004, 25, 1605–1612. (42) Morris, G. M.; Huey, R.; Lindstrom, W.; Sanner, M. F.; Belew, R. K.; Goodsell, D. S.; Olson, A. J. AutoDock4 and AutoDockTools4: automated docking with selective receptor flexibility. J. Comput. Chem. 2009, 30, 2785–2791. (43) Trott, O.; Olson, A. J. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J. Comput. Chem. 2010, 31, 455–461. (44) Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh ewald. J. Chem. Theory Comput. 2013, 9, 3878–3888. (45) Le Grand, S.; Gotz, A. W.; Walker, R. C. SPFP: Speed without compromise?A mixed precision model for GPU accelerated molecular dynamics simulations. Comput. Phys. Commun. 2013, 184, 374 – 380. (46) Case, D.; Cerutti, D.; Cheatham, T. I.; Darden, T.; Duke, R.; Giese, T.; Gohlke, H.; Goetz, A.; Greene, D.; Homeyer, N. et al. AMBER 2014. 2014; University of California, San Francisco. (47) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Struct., Funct., Bioinf. 78, 1950–1958. (48) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general AMBER force field. J. Comput. Chem. 2004, 25, 1157–1174. (49) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Comput. Chem. 1983, 79, 926–935. (50) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684–3690.

31

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(51) Yu, Z.; Chao, W.; Lifeng, Z.; Niu, H. Exploring the early stages of the pH-induced conformational change of influenza hemagglutinin. Proteins: Struct., Funct., Bioinf. 82, 2412– 2428. (52) Guan, S.; Wang, T.; Kuai, Z.; Qian, M.; Tian, X.; Zhang, X.; Yu, Y.; Wang, S.; Zhang, H.; Li, H. et al. Exploration of binding and inhibition mechanism of a small molecule inhibitor of influenza virus H1N1 hemagglutinin by molecular dynamics simulation. Sci. Rep. 2017, 7, 3786. (53) Sondergaard, C. R.; Olsson, M. H. M.; Rostkowski, M.; Jensen, J. H. Improved treatment of ligands and coupling effects in empirical calculation and rationalization of pKa values. J. Chem. Theory Comput. 2011, 7, 2284–2295. (54) Olsson, M. H. M.; Søndergaard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions. J. Chem. Theory Comput. 2011, 7, 525–537. (55) Dolinsky, T.; Nielsen, J.; McCammon, J.; Baker, N. PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. Nucleic Acids Res. 2004, 32, W665– W667. (56) Mann, H. B.; Whitney, D. R. On a test of whether one of two random variables is stochastically larger than the other. The Annals of Mathematical Statistics 1947, 18, 50–60. (57) Shapiro, S. S.; Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 1965, 52, 591–611. (58) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J Mol Graph. 1996, 14, 33–38. (59) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W. et al. Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models. Acc. Chem. Res. 2000, 33, 889–897.

32

ACS Paragon Plus Environment

Page 32 of 34

Page 33 of 34 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(60) Greenidge, P. A.; Kramer, C.; Mozziconacci, J.-C.; Wolf, R. M. MM/GBSA Bbinding energy prediction on the PDBbind data set: successes, failures, and directions for further improvement. J. Chem. Inf. Model. 2013, 53, 201–209. (61) Cerón-Carrasco, J. P.; Perez-Sanchez, H.; Zuñiga, J.; Requena, A. Antibodies as carrier molecules: encapsulating anti-inflammatory drugs inside herceptine. J. Phys. Chem. B 2018, 122, 2064–2072. (62) Cerón-Carrasco, J. P.; Cerezo, J.; Requena, A.; Zuñiga, J.; Contreras-García, J.; Chavan, S.; Manrubia-Cobo, M.; Pérez-Sánchez, H. Labelling herceptin with a novel oxaliplatin derivative: a computational approach towards the selective drug delivery. J. Mol. Model. 2014, 20, 2401. (63) Zhao, Y.; Truhlar, D. G. Screened-exchange density functionals with broad accuracy for chemistry and solid-state physics. Theo. Chem. Acc. 2008, 120, 215–241. (64) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999–3094. (65) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; Contreras-García, J.; Cohen, A. J.; Yang, W. Revealing noncovalent interactions. J. Am. Chem. Soc. 2010, 132, 6498–6506. (66) Contreras-García, J.;

Johnson, E. R.;

Keinan, S.;

Chaudret, R.;

Piquemal, J.-P.;

Beratan, D. N.; Yang, W. NCIPLOT: a program for plotting noncovalent interaction regions. J. Chem. Theory Comput. 2011, 7, 625–632.

33

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Graphical TOC Entry Peptide inhibitors

34

ACS Paragon Plus Environment

Page 34 of 34