key: cord-0336068-609tdqpx authors: Buthelezi, Ndumiso M; Amoako, Daniel G.; Somboro, Anou M.; Khan, Rene B.; Kumalo, Hezekiel M. title: A molecular dynamic investigation of Human Rhinovirus 3C Protease Drug Target: Insights towards the design of potential Inhibitors date: 2022-03-07 journal: bioRxiv DOI: 10.1101/2022.03.06.483152 sha: 742789ec5ffedb708a25f1053a41ca40642a8a4b doc_id: 336068 cord_uid: 609tdqpx The 3C protease is distinguished from most proteases due to the presence of cysteine nucleophile that plays an essential role in viral replication. This peculiar structure encompassed with its role in viral replication has promoted 3C protease as an interesting target for therapeutic agents in the treatment of diseases caused by human rhinovirus (HRV). Herein we present a comprehensive molecular dynamics study of the comparison of two potent inhibitors, sg85 and rupintrivir complexed with HRV-3C protease. The binding free energy studies revealed a higher binding affinity for sg85 −58.853 kcal/mol than for rupintrivir −54.0873 kcal/mol and this was found to be in correlation with the experimental data. The energy decomposition analysis showed that, residues Leu 127, Thr 142, Ser 144, Gly 145, Tyr 146, Cys 147, His 161, Val 162, Gly 163, Gly 164, Asn 165, Phe 170 largely contributed to the binding of sg85, whereas His 40, Leu 127 and Gly 163 impacted the binding of rupintrivir. It further showed that His 40, Glu 71, Leu 127, Cys 147 Gly 163 and Gyl 164 are crucial residues that play a key role in ligand-enzyme binding; amongst these residues are residues of the conserved active site (His 40, Glu 71 and Cys 147). These findings provide a comprehensive understanding of the dynamics and structural features and will serve as guidance in the design and development of potent novel inhibitors of HRV. Graphical Abstract Human rhinovirus (HRV), also known as the causative agent of common cold [1, 2] , are nonenveloped viruses that contain a single-strand ribonucleic acid (ssRNA) genome enclosed in an icosahedral (20-sided) capsid [3] . It is estimated that 18.8 billion and 150 million people were infected with upper respiratory and lower respiratory infections, respectively [4] . HRV is also reported by the World Health Organization (WHO) to be responsible for approximately 2 million children deaths necessitating the search and improvement of new potent drugs. HRV is categorized into 3 species, namely, RV-A, RV-B, and RV-C [5] . RV-C being recently identified and studies show that HRV-3C-is most prevalent HRV identified in hospitalized children [6] . HRV is a positive-sense ssRNA virus of approximately 7200 bp [7] . Upon infection, the positive-strand RNA genome of HRV-3C is translated into a large poly-protein that is required for the production of new infectious virion and is dependent on two virally encoded proteases 2A and 3C. [8, 9] . The 3C protease plays an indispensable role in viral replication in proteolytic cleavage of large polyproteins to functional proteins, and enzymes required for viral RNA replication both structurally and enzymatically [10] . The enzyme possess a strictly conserved catalytic active site as shown in Fig.1 in blue, thus, making it an attractive target for therapeutic agents in the treatment of diseases caused by human rhinovirus [11, 12] . Page | 5 Several small molecule HRV inhibitors have been identified, with 3C protease (3CP) emerging as a promising therapeutic target [12] . It has recently been established that sg85 and rupintrivir inhibits HRV-3C protease. Rupintrivir (formerly AG7088) is an irreversible inhibitor of the 3C protease of the human rhinovirus (HRV) that was discovered using structure-based drug design techniques [13] . Clinical trials reports showed that therapeutic application of rupintrivir was able to reduce viral load and moderate the severity of illness in a human experimental HRV challenge trial, hence providing a proof of concept for the mechanism of 3C protease inhibition. It was later discovered that rupintrivir had no significant effect on virus replication or disease severity, as such rupintrivir's clinical development was halted. [14] . Sg85 is a peptidomimetics with Michael acceptor warheads permanently disable the protease by covalent binding to its catalytic site and is the most promising candidate [15] . In cell-based tests, Sg85 suppresses the replication of enteroviruses by inhibiting the EV68 3C-protease. [16] . It is believed that inhibition of picornavirus replication reduces the severity and shortens the duration of cold symptoms while rupintrivir is an irreversible inhibitor of HRV-3C protease and is regarded as the most potent novel peptidomimetic with inhibition activity against HRVs [17] . However, in terms of binding efficiency at the molecular level, there is no proof or comparison between these two inhibitors. Hence, the molecular understanding of HRV-3C protease structure and ligand-binding mechanisms remains unclear. Molecular dynamics simulations are one of the useful tools that have the capacity to provide a concise understanding at atomic level of complex systems [18] . To study the binding of the two inhibitors with HRV-3C protease, we employed molecular dynamics simulations (MD). In recent years, different post-dynamics analyses such as root mean square fluctuation, principal component analysis, Mechanics/Generalized Born Surface Area and dynamic cross-correlation analyses have been applied extensively in understanding the dynamics of biological complex systems at the atomic level. Therefore, this work provides understanding of molecular dynamics of HRV-3C protease complexed with sg85 and rupintrivir from a computational point of view. Furthermore, we aim at providing insight into the binding-landscape of the two inhibitors in HRV-3C protease by conducting a broad comparative MD analysis of the HRV-3C-Free, HRV-3C-Sg85 and HRV-3C-Rupintrivir systems. We hope that this work could serve as a cornerstone in drug design and developments of more potential drugs against HRV protease, which could lead the fight against the virus. Page | 7 The X-ray crystal structure of human rhinovirus-3C protease, PDB code 5FX6, bound to rupintrivir was extracted from Protein Bank Database (www.rcsb.org) [19] .The Michael acceptor, Sg85, was obtained from the PDB 2YNB where it was in complex with corona virus HKU4 [20] . Prior to molecular docking all non-standard residues including ions, Na+, cl-, H2O were removed from the structure of human rhinovirus-3C protease. Subsequently, Sg85 was removed from the 2YNB complex and optimized using Avogadro software before docking into the binding site of human rhinovirus 3C protease using Autodock vina integrated in UCSF Chimera. All missing residues were remodeled using the modeler plugin incorporated in UCSF Chimera graphical user interface [21] . In preparation for MD simulation, hydrogen atoms and amber charges were then added to the ligand and saved separately in mol2 format whiles the protein was saved in pdb format. Molecular visualization of the Ligand and receptor were conducted in Chimera and Avogadro softwares [22] . The simulation setup was done using the Graphical Processing Unit (GPU) version of Particle Mesh Ewald Molecular Dynamics (PMEMD) integrated in the Amber 14 suite [23, 24] . Atomic partial charges were generated for the ligands using the General Amber Force Field GAFF force field during the ANTECHAMBER stage on the AMBER 14 [25] . Likewise, the FF14SB force field was used to parametize the receptor [25] [26] . This was then followed by the LEAP phase, where additional hydrogen atoms were added to the system. Also, neutralization of the system was carried out through the addition of counter ions such as NA + and Cl -. The system was then subjected to solvation by immersing in an orthorhombic box with TIP3P [27] water molecules at a distance of 10 Å. In our previous work, system setup, minimizations, heating and equilibration steps were thoroughly explained [28, 29] . Subsequently the system was subjected to a 200ns MD simulation time, and all trajectories generated during the simulation run were then saved for every 1ps and analyzed. The CPPTRAJ and PTRAJ modules [30] incorporated in the Amber 14 suite were used for analyzing the trajectories: Parameters such as RMSD, RMSF, RoG, DCC, PCA were determined. The Origin data analysis tool (www.originlab.com) was utilized for all graphical plots while Chimera [21] was used for visualizations. Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) was used to calculate the binding free-energy profiles of rupintrivir and sg85 bound to HRV-3C protease. [31, 32] . The binding free energy calculation gives insights into the association of the protein-ligand in a complex and also gives an account on the end point energy calculation. Diffirential binding free energy was calculated taking into account 1000 snapshots from 200 ns trajectories. The calculation of the binding free energy is described by the following set of equations: Gsol = GGB/PB + GSA (4) Where Egas represents gas-phase energy, internal energy is represented by Eint, whiles Eele and Evdw represent the electrostatic and Van der Waals contributions, respectively. Egas is the gas phase, which is directly elevated from the FF14SB force terms. Gsol, which stands for solvation energy, can be divided into polar and nonpolar contribution states. The polar solvation contribution, GGB, is calculated by solving the GB equation, whereas the nonpolar solvation contribution, a water probe radius of 1.4 Å is used to determine GSA, is estimated using a. T and S represent temperature and total solute entropy, respectively. A per-residue decomposition analysis of the interaction energy for each residue was performed using the MM/GBSA method to determine the contribution of each amino acid to the total binding free energy profile between rupintrivir and sg85. The dynamic cross correlation (DCC) method has been widely used to calculate the correlation coefficients of motions between atoms in a protein. [33] . The CPPTRAJ module was used to calculate the residue-based fluctuations during the simulation [30] . DCC is represented by the following formula given below: With Cij being the cross-correlation coefficient, which varies with a range of -1 to +1. A fully correlated and anti-correlated motion represents the upper and lower limits during the simulation, respectively. Where i and j represent the i th and j th residues and ri and rj symbolize displacement vectors correspond to ith and jth residues, respectively. In this study, DCC calculations were carried out considering the backbone C- atomic fluctuations. Principal component analysis (PCA) also referred to as essential dynamics (ED) and is one of the utmost unconventional methods for trajectory analysis. PCA has proven to be powerful and robust, opening up new avenues for visualizing and exploring the dynamics of protein cavities. [34, 35] . PCA describes the eigenvectors and eigenvalues, which represent the direction of motions and the amplitudes in those directions of the protein, respectively [36] . The root-mean-square deviation (RMSD) was plotted to identify the stability of studied enzyme- Furthermore, Fig.S2 gives information about the β-hairpin in the active site residues in the regions 120-138 of all systems. It has been proposed that the β-hairpin acts as launch sites in initial events of protein folding [38] . Herein, the evolutions of the β-hairpin during a MD simulation for all conformations were examined and snapshots along the trajectories of MD simulations are shown in Fig.2 . Page | 12 HRV-hairpin. Experimental results have shown Sg85 to be active against a wide range of HRV isolates suggesting that it has a higher binding affinity for the binding site compared to Rupintrivir [42] . We can conclude from the results that location of Sg85 in the active site leads to a conformational rigidity and β-hairpin stability [15] . Radius of gyration (Rg) was computed to measure the compactness of the protein structure and additionally to give details into stability of the complexed systems [43, 44] . An Rg plot in Fig.S3 was used to analyze the overall protein dimensions of the HRV systems. From Fig.S3 the similarities between the three systems are obvious, as evidenced by the similarity in residual alignment within the secondary and tertiary structures. HRV-3C-Free and HRV-3C-Rupintrivir displayed a slighter higher Rg compared to HRV-3C-Sg85. These results are in accordance with RMSF, which justified a substantial increase in biomolecular flexibility of the HRV-Free and HRV-3C-Rupintrivir when compared to the HRV-3C-Sg85. This indicates that HRV possesses a rigid structural stability when it is bound to Sg85 as compared to. A theoretical explanation for this observation is that binding of Sg85 to the active site of HRV causes conformational rigidity, which inhibits HRV's conformational flexibility. The dynamic cross correlation was applied to quantify the correlation coefficients of motions of Anti-correlated residual motion in Fig.3 , for and sg85 occur between residues 50 -70 relative to 110, while a strong correlated residual motion in case of rupintrivir occur between 20 -80 relative to each other. The correlation map displays a strong correlation of residual motion in the HRV-3C-Rupintrivir conformation, while in the HRV-3C-Sg85 conformation there was an observed greater anti-correlated motions (negative correlations) during simulation time. This is supported by the RMSF results, which also shows a strong interaction with residues compared to HRV-3C-Rupintrivir. conformations of HRV-3C. Fig.4 illustrates that the Sg85 (in black) showed to have the largest surface area followed by Rupintrivir (in red) when compared to Free (in green). This could imply that the Sg85 is less flexible allowing a more easy binding toward the HRV enzyme compared to Rupintrivir. As evident in Fig.S2 and Fig.2 , the HRV-3C-Free showed to be more flexible in contrast to Rupintrivir and Page | 17 Sg85. The evidence in Fig.4 unveils basic insights of the dynamic behavior of biological systems at molecular level. The overall binding free energy for both Sg85 and Rupintrivir was calculated using the MM/PBSA method. Per-residue energy decomposition was computed to gain insight to the ligand-residue interaction and this was accomplished by 1000 snapshots of the 200 ns simulation trajectories. Rupintrivir was estimated to be -58.853 kcal/mol and -54.087 kcal/mol. The binding free energy of Sg85 is higher by (~4 kcal/mol) than that of Rupintrivir, which further suggests that Sg85 has a more intimate interaction with the active site residues compared to sg85. This is in great accordance with recent publications that have rated Sg85 as the most potent inhibitor of HRV enzyme [15] . As revealed from the per-residue energy contribution plots shown in Fig.5 The current study highlights the most important amino acids involved in inhibitor binding. The generated pharmacophore library contains only highly contributing amino acid residues, as determined by free binding energy contributions obtained from molecular dynamic (MD) simulations. [45, 46] . We trust that the findings from this study serve as a foundation towards the design of new inhibitors that will interact with catalytic triad residues to discover novel compounds to inhibit HRV using a pharmacophore model. Ligand-protein interaction and binding at molecular level is an imperative biological process. In this report, we have provided insights into the binding mechanism of the rupintrivir and sg85 on 3C protease with the application of advance molecular dynamics simulation tools. It has been experimentally reported that sg85 is a potent inhibitor when compared to rupintrivir but to this end, hence this comparative study was conducted using various advanced MD simulation and post-analysis tools to gain understanding into the binding landscape and dynamic structural features associated with the binding of rupintrivir and sg85 on 3C protease. From the results obtained, sg85 induced a more stable protein structure thus suggesting a more intimate binding and interaction of sg85 with the 3C protease when compared to rupintrivir as seen from the RMSD, RMSF, RG, DCC and PCA plots. The binding free energy calculation revealed a higher binding affinity for sg85 -58.853 kcal/mol than for rupintrivir -54.087 kcal/mol and this is in correlation with the experimental data. This thus suggests that sg85 can be used as scaffold to design new inhibitors. The energy decomposition analysis showed that, residues Leu 127, Thr Leu 127, Cys 147 Gly 163 and Gyl 164 are crucial residues that play a key role in ligand-enzyme binding; amongst these residues are residues of the conserved active site (His 40, Glu 71 and Cys 147). Therefore, from these results we believe that the residue interactions and the binding free energy analysis must be considered for a pharmacophore model and structure-based design of novel and potent inhibitors of HRV. Authors declare no conflict of interest. How does rhinovirus cause the common cold cough? Antigenic Diversity, and Advancements in the Design of a Human Rhinovirus Vaccine Classification and Evolution of Human Rhinoviruses Recurrent rhinovirus infections in a child with inherited MDA5 deficiency Molecular Epidemiology of Human Rhinoviruses and Enteroviruses Highlights Their Diversity in Sub-Saharan Africa. Viruses Viral proteases as targets for drug design Picornaviruses and nuclear functions: targeting a cellular compartment distinct from the replication site of a positive-strand RNA virus Roles of the Picornaviral 3C Proteinase in the Viral Life Cycle and Host Cells Proteases of Human Rhinovirus: Role in Infection BT -Rhinoviruses: Methods and Protocols Activity of the Human Rhinovirus 3C Protease Studied in Various Buffers, Additives and Detergents Solutions for Recombinant Protein Production Protease with Potent Antiviral Activity against Multiple Rhinovirus Serotypes Source : Proceedings of the National Academy of Sciences of the United States of Amer Vitro Resistance Study of Rupintrivir, a Novel Inhibitor of Human Rhinovirus 3C Protease The Enterovirus 3C Protease Inhibitor SG85 Efficiently Blocks Rhinovirus Replication and Is Not Cross-Resistant with Rupintrivir The enterovirus 3C protease inhibitor SG85 efficiently blocks rhinovirus replication and is not cross-resistant with rupintrivir The Enterovirus Protease Inhibitor Rupintrivir Exerts Cross-Genotypic Anti-Norovirus Activity and Clears Cells from the Norovirus Replicon Towards exact molecular dynamics simulations with machine-learned force fields Design and structure-activity relationships of novel inhibitors of human rhinovirus 3C protease Structure and inhibition of EV-D68, a virus that causes respiratory illness in children UCSF Chimera--a visualization system for exploratory research and analysis Avogadro: an advanced semantic chemical editor, visualization, and analysis platform Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald The General AMBER Force Field (GAFF) Can Accurately Predict Thermodynamic and Transport Properties of Many Ionic Liquids Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB Comparison of simple potential functions for simulating liquid water Possible allosteric binding site on Gyrase B, a key target for novel anti-TB drugs: homology modelling and binding site identification using molecular dynamics simulation and binding free energy calculations A comparative molecular dynamics study on BACE1 and BACE2 flap flexibility Software for Processing and Analysis of Molecular Dynamics Trajectory Data Calculating structures and free energies of complex molecules: combining molecular mechanics and continuum models Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding mDCC_tools: characterizing multimodal atomic motions in molecular dynamics trajectories Principal Component Analysis reveals correlation of cavities evolution and functional motions in proteins Automatic Identification of Mobile and Rigid Substructures in Molecular Dynamics Simulations and Fractional Structural Fluctuation Analysis From Principal Component to Direct Coupling Analysis of Coevolution in Proteins: Low-Eigenvalue Modes are Needed for Structure Prediction Measuring and comparing structural fluctuation patterns in large protein datasets Formation and stability of beta-hairpin structures in polypeptides Understanding β-Hairpin Formation by Molecular Dynamics Simulations of Unfolding Long range Trp-Trp interaction initiates the folding pathway of a pro-angiogenic β-hairpin peptide Correct folding of an α-helix and a βhairpin using a polarized 2D torsional potential In vitro activity of pleconaril and AG7088 against selected serotypes and clinical isolates of human rhinoviruses Effect of molecular weight and temperature on physical aging of thin glassy poly(2,6-dimethyl-1,4-phenylene oxide) films Molecular dynamics study of Zn(abeta) and Zn(abeta)2 Per-Residue Energy Footprints-Based Pharmacophore Modeling as an Enhanced In Silico Approach in Drug Discovery: A Case Study on the Identification of Novel β-Secretase1 (BACE1) Inhibitors as Anti-Alzheimer Agents Per-residue energy decomposition pharmacophore model to enhance virtual screening in drug discovery: a study for NMB acknowledges School of Health Sciences, UKZN, for financial assistance and the Center of High Performance Computing (CHPC), for computational facility (www.chpc.ac.za).