key: cord-0898685-f7cqp80q
authors: Sargolzaei, Mohsen
title: Effect of nelfinavir stereoisomers on coronavirus main protease: molecular docking, molecular dynamics simulation and MM/GBSA study
date: 2020-12-04
journal: J Mol Graph Model
DOI: 10.1016/j.jmgm.2020.107803
sha: fb3907434488a97eb8734097145b1b0c2c246902
doc_id: 898685
cord_uid: f7cqp80q

In this study, the binding strength of 32 diastereomers of nelfinavir, a proposed drug for the treatment of COVID-19, was considered against main protease. Molecular docking was used to determine the most potent diastereomers. The top three diastereomers along with apo form of protein were then considered via molecular dynamics simulation and MM-GBSA method. During the simulation, the structural consideration of four proteins considered was carried out using RMSD, RMSF, Rg and hydrogen bond analysis tools. Our data demonstrated that the effect of nelfinavir RSRSR stereoisomer on protein stability and compactness is higher than the other. We also found from the hydrogen bond analysis that this important diastereomer form three hydrogen bonds with the residues of Glu166, Gly143 and Hie41. MM/GBSA analysis showed that the binding strength of RSRSR is more than other stereoisomers and that the main contributions to binding energy are vdW and electronic terms. The nelfinavir RSRSR stereoisomer introduced in this study may be effective in the treatment of COVID-19.

The COVID-19 has been a major cause of death in recent days. The worldwide number of deaths reached 431,192 by 15 June 2020, with 7,805,148 laboratory-confirmed cases [1] . Virus infection affects the function of the lungs, the digestive system, and the central nervous system [2] [3] [4] [5] . A novel coronavirus strain associated with fatal respiratory disease was reported at the end of 2019 [6] . This pathogen was temporarily named coronavirus by the World Health Organization (WHO) in 2019 [7] . Coronaviruses are single-stranded positive-sense RNA viruses with large viral RNA genomes [8] . Main protease is a major enzyme for the reproduction of the SARS-Cov-2.

Proteases play a key role in the replication of a number of viruses. These enzymes often serve as protein targets for antiviral therapy development [9] . The main protease of SARS-Cov-2 is similar to SARS-Covid with 96 per cent identity. Scheme 1. Nelfinavir and its five chiral centers. These five chiral centers lead to 32 possible stereoisomers for nelfinavir. Since each chiral center is labeled R or S label, each stereoisomer is named with five labels for chiral centers of 1 to 5, respectively.

Nelfinavir is an antiviral drug previously used to treat HIV infection [18] . Using homology modelling, molecular docking and binding free energy calculation, Xu et al., proposed nelfinavir as a SARS-CoV-2 main protease inhibitor [19] .

The experimental test of nelfinavir against SARS-CoV infected cells showed that the antiviral activity of this drug is high (EC50=0.048 µM) [20, 21] . Nelfinavir is a chiral molecule with molecular weight of 567.789 g/mol, 12 rotatable bonds and five chiral centers (see Scheme 1) . Since each chiral center may accept R or S stereoisomers, 32 configurations are produced for nelfinavir. Drug chirality is an important property because different stereoisomers have different pharmacodynamic, pharmacokinetic, and toxicological properties [22, 23] . In recent years, the pharmaceutical industry has shown a tendency to make new drugs in a single enantiomeric form. This approach, known as the chiral switch, has allowed the marketing of many drugs in a specific stereochemical configuration [24] [25] [26] [27] [28] [29] . Although nelfinavir, a 32-stromisomeric molecule, has been introduced as a candidate drug for treatment of Covid-19, no studies have been conducted to study the binding of nelfinavir stereoisomers to main protease. Structure of main protease of SARS-Cov-2 along with its binding cavity.

In this study, all possible nelfinavir stereoisomers are constructed and their binding strength to the protease enzyme (see Figure 1 and Figure 1S ) is studied via molecular docking. In addition, molecular dynamic simulation and MM/GBSA analysis are performed to consider dynamic behavior and binding analysis of the most important stereoisomer complexes, respectively.

The B3LYP method [30] and the 6-31G* basis set were used to optimize nelfinavir stereoisomers. Restrained Electrostatic Potential (RESP) fit method [31] was used to obtain atomic charges for each stereoisomer. The electrostatic potential of the optimized structures was calculated using the Hartree-Fock method and the 6-31G* base set using the Gaussian 03 package [32] . RESP charges were derived for optimized ligands using the AMBER program.

The 3D coordinate of Coronavirus main protease was derived from experimental pdb structure (PDB code: 5R81). We selected active site around the Z1367324110 ligand as ligand cavity for docking process. Molegro Virtual Docker (MVD) v6.0 [33] was used to perform docking. Piecewise Linear Potential (PLP) scoring functions were used to derive scoring function. For all compounds considered, the molecular docking was carried out with a grid resolution of 0.30 Å and a binding site radius of 17 Å. The crystallized ligand was used to define the search space as a reference ligand.

To search, 10 runs were applied using a maximum of 2000 iterations with a total population size of 50. Pose generation was done using threshold energy of 100. The simplex evaluation was completed with a maximum of 300 steps and the neighbor distance factor 1.

The top three stereoisomers were selected for the MD simulation based on the docking score calculated. Also, MD simulation was performed on the main protease apo form. MD simulation was performed using software package AMBER16 [34] . FF14SB force field [35] was used for protein. In addition, the GAFF force field [36] was applied for each stereoisomer.

Four complexes were dissolved into a rectangular box of TIP3P [37] water molecules and neutralized by the addition of an appropriate number of counter ions. Heating of all systems was done from 0 to 300K for 0.1ns after the minimization process. The density of the water molecules was relaxed for 0.1 ns at a constant pressure of 1atm and a temperature of 300K. Temperature control during the simulation was done using Langevin algorithm [38] with a collision frequency of 2ps -1 . All simulations were performed at a temperature of 300K for 100ns.

The analysis of the derived trajectories was carried out using the CPPTRAJ program [39] . Hydrogen bond analysis was also carried out during the simulation using the CPPTRAJ program. The MM-PBSA.py program [40] was used for the MM/GBSA [41, 42] and the pairwise free energy decomposition analysis. LIGPLOT was used to generate 2D ligand-protein interaction diagrams [43, 44] . J o u r n a l P r e -p r o o f binding site. The root-mean-square deviation (RMSD) values between the experimental structure of ligand and predicted ligand pose were then determined. The protocol used for our docking produced a ligand pose with RMSD of 1.15Å from experimental structure of ligand (see Figure 4 ). It has been shown that docking protocols that can return poses with RMSD value below 2Å are correct and reliable [45] . The comparison of our RMSD obtained with the threshold RMSD indicated above shows that we can continue research using our validated docking protocol. Using the above docking protocol, 32-stereoisomers shown in Figure 2 were docked against coronavirus main protease. The results of the moldock score obtained are shown in Table 1 . The results show that the moldock score varies from -178.6 for RSRSR to -147.2 for RSSRS stereoisomers. The docking results show that the residues of Glu166, Met165, Hie41, Pro168, Leu167, Hie164, Cys145, GLY143, Thr190 and Asp187 are important residues in protein binding site. Figure 2S in the supporting information shows the nelfinavir docked stereoisomers to the main protease binding site. It has been shown that an important part, called the anchor site (see Figure 5 ), exists at the binding site of the main protease. The binding stability of the ligand inside the main protease pocket can be significantly improved if part of the ligand is located at the anchor site [46] . The existence of the anchor site is also verified with consideration of the position of the crystallized ligand in the experimental pdb structure of the coronavirus main protease. As shown in Figure 4 , the position of ligand Z1367324110 is at the anchor site of the main protease.

If we look at the orientation of 32 stereoisomers, we see that the RSRSR stereoisomer with the highest dock score has two tert-butyl and benzene groups within the main protease anchor site, while the other docking poses do not have those groups within the anchor site. On the other hand, only one portion of the ligand structure is inside the main protease anchor site for those stereoisomers with the lowest docked score.

Among the 32 stereoisomers considered in this study, according to our nomenclature, stereoisomer of SRRRR is the FDA-approved nelfinavir stereoisomer for the treatment of HIV disease. According to Table 1 , the SRRRR stereoisomer ranks twenty-nine in the docking score. Other researchers used only FDA-approved nelfinavir stereoisomer to study ligand interaction with main protease [47, 48] . Given that the structure of HIV protease differs from the main protease of coronavirus [49] , it is not expected that the FDA approved nelfinavir stereoisomer for HIV protease will be the only proper stereoisomer for the main protease of coronavirus. Anyway, our data shows that docked stereoisomer of SRRRR is bound inside the anchor site of main protease with its benzene group. This finding is consistent with the finding in ref [46] .

To further investigate the results of molecular docking, the top three nelfinavir stereoisomers, namely RSRSR, SRRSR and SSRRS, as well as the apo form of protein, were subjected to 100ns MD simulation. In order to investigate the stability of molecular dynamic simulations, RMSD data were plotted for three stereoisomer complexes as well as protein apo form of protein. In addition, mean RMSD values for RSRSR complex, SRRSR [50, 51] . They have shown that main protease is not unfolded during the simulation and its structure is stable. Our data shows that all of the considered complexes reach the equilibrium region after about 50 ns. According to this finding, we can use the equilibrium region of simulation for thermodynamics analysis. Another result that can be derived from the RMSD analysis is the effect of the stereoisomers studied on protein stability. This finding is derived by comparing RMSD of apo form of protein with RMSD of other three stereoisomers complexes. Our data show that the effect of RSRSR on protein stability is more than two others. On the other hand, the ligand of SSRRS has the lowest effect on protein stability.

Residual flexibility of a protein can be considered using RMSF data. Figure  3b shows the RMSF plot for apo form of protein and also complexes of RSRSR, SRRSR and SSRRS. The results of the RMSF analysis show that fluctuation of residues in three stereoisomer complexes increase with respect to apo protein. In some residue regions, such as 270-280, the value of RMSF increases slightly from 2Å, which is related to fluctuation in turn and coil structures. We can see that all the complexes remained stable throughout the simulations. Other researchers have also shown that main protease is not a high-flexible protein. [51] .

The compactness, stability and folding of protein are determined by gyration radius. Figure 3c shows Rg for apo protein and three other complexes versus time. Data shows that the protein remains compact for all structures during the simulation process. This finding is consistent with the results reported for Rg of coronavirus main protease [51] . Also, we can see that all ligands reduce the compactness of proteins slightly. The comparison of results shows that the effect of RSRSR on compactness of protein is more than other stereoisomers. Table 2 shows the hydrogen bond analysis of protein residues and ligands studied. In complex RSRSR, three hydrogen bonds are formed between Glu166 and ligand for 54, 34 and 5 percent of the simulation time. Also, Gly143 forms a hydrogen bond with RSRSR for 24 percent of simulation time. Hie41 is also a residue that is involved in ligand hydrogen bonding within 7 percent of the simulation. On the other hand, ligand SRRSR forms two hydrogen bonds with protein; one with Ser46 with a long lifetime of 57% during the simulation and the other with Hie41 with a very short life time. Ligand SSRRS forms a hydrogen bond with Glu166 for 21 percent of simulation time.

Based on the MD trajectories, the MM/GBSA method predicts binding free energy for ligand-receptor of RSRSR, SRRSR and SSRRS complexes (see Table 3 ). MM/GBSA analysis was done on the last 50 ns. According to the energy components of the binding free energies, the major favorable contributions to the stereoisomer binding are the van der Waals and electrostatic energies for all complexes. On the contrary, the term polar solvation free energy (∆EGB) is largely unfavorable for binding in five complexes. The terms van der Waals and electrostatic energies promote binding and can offset the negative effect of polar solvation free energy. The total calculated free binding energy (∆GBinding) against main protease for compound RSRSR is greater than other compounds. The order of binding free energy derived from the MM/GBSA analysis is consistent with the results of our docking study.

The vdW and electrostatic MM/GBSA values correlate with the number of hydrogen bonds. Thus, the more hydrogen bonding, the lower the electrostatic interaction energy. Additionally, because hydrogen bonding, in a molecular mechanics sense, is the interaction between positive and negative charges, the van der Waals interactions will also be more favorable resulting in lower vdW energies.

On the other hand, nelfinavir is very non-polar and is not very water soluble. Taking the hydrogen bonding analysis in Table 2 and the solvation energy values (∆GSOL and the components ∆EGB and ∆ESURF) in Table 3 into account, three nelfinavir stereoisomers have unfavorable solvation (∆Gsol) energies and the binding be driven by molecular interactions (vdW and electrostatic interactions).

MM/GBSA pairwise per-residue free energy decomposition calculations were used to find a detailed picture of the dominant residues involved in the binding process (see Figure 7) . Results show that a number of residues interact with ligand in RSRSR complex, and that their contribution of some of them to free binding energy is as follows: Glu166>Gly143>Hie41>Asp187>Met165>Hie164>Cys145>Thr190>Leu167> Pro39. Glu166 is located in the ligand cavity on extended secondary structure. The contribution of residues in the SRRSR complex to ligand interaction is as follows: Hie41> Ser46> Glu166> Met165> Met49> Thr25> Cys145> Leu27>Thr26>Gly143. The significant residue for interaction with SRRSR is Hie41, which is located on the secondary structure of the turn. The following residues interact with ligand in the SSRRS complex: Gln189 > Asp187 > Glu166> Hie 164>Val186>Leu167>Pro168> Pro168> Met165> Gln192. 

Gln189 is an important residue for interaction with ligand in the SSRRS complex.

Comparison of interaction energy based on the major nelfinavir stereoisomer-residue pair.

In this study, the effect of 32 nelfinavir stereoisomers was considered against coronavirus main protease. Our docking results showed that RSRSR, SRRSR and SSRRS stereoisomers are the top three stereoisomers for binding to protein. Molecular dynamics simulation was done for these three stereoisomer complexes and RMSD analysis demonstrated that main protease is a stable enzyme with small structural variations. Also, from RMSD analysis, we found that all complexes reach the equilibrium structure after approximately 50 ns and RMSD of RSRSR stereoisomer complex is more than two other complexes. RMSF analysis of the trajectory of all complexes showed that the main protease does not have a significant flexible region in its structure. In addition, RMSF data showed that the effect of RSRSR stereoisomer on flexibility of protein in some regions is more than other ones. Based on the Rg plot analysis, it was founded that main protease is a fold and compactness protein and that the effect of RSRSR ligand on compactness of protein is more than other ligands. Hydrogen bond analysis of the simulations demonstrated that RSRSR, SRRSR and SSRRS stereoisomers forms three, two and one hydrogen bond with main protease, respectively. MM/GBSA analysis demonstrated that van der Waals and electrostatic energies are the major favorable contributions to binding energy for three ligands and that binding free energy of RSRSR is more than two others. The calculated binding free energies from MM/GBSA analysis were well correlated with moldock score energies for three stereoisomers. MM/GBSA pairwise per-residue free energy decomposition analysis showed that Glu166, Gly143, Hie41 residues are important residues interacting with ligand RSRSR. Based on the findings of this study, further research is recommended for consideration of effect of nelfinavir RSRSR stereoisomer on COVID-19. 

World Health Organization. Coronavirus disease 2019 (COVID-19)

The evidence of porcine hemagglutinating encephalomyelitis virus induced nonsuppurative encephalitis as the cause of death in piglets

Axonal Transport Enables Neuron-to-Neuron Propagation of Human Coronavirus OC43

Evolution of the novel coronavirus from the ongoing Wuhan outbreak and modeling of its spike protein for risk of human transmission

Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China, The Lancet

Cross-species transmission of the newly identified coronavirus 2019-nCoV

Emerging coronaviruses: Genome structure, replication, and pathogenesis

Antiviral Drug Discovery: Norovirus Proteases and Development of Inhibitors

Potential inhibitors for 2019-nCoV coronavirus M protease from clinically approved medicines, bioRxiv

Potential Inhibitor of COVID-19 Main Protease (Mpro) From Several Medicinal Plant Compounds by Molecular Docking Study

Statins and the COVID-19 main protease: in silico evidence on direct interaction

An investigation into the identification of potential inhibitors of SARS-CoV-2 main protease using molecular docking study

Prediction of the SARS-CoV-2 (2019-nCoV) 3C-like protease (3CL (pro)) structure: virtual screening reveals velpatasvir, ledipasvir, and other drug repurposing candidates

Structural basis of SARS-CoV-2 3CLpro and anti-COVID-19 drug discovery from medicinal plants

Computational Design of ACE2-Based Peptide Inhibitors of SARS-CoV-2

Binding site analysis of potential protease inhibitors of COVID-19 using AutoDock, VirusDisease

Conformations and Three-Dimensional Structures of Selected SARS-CoV-2 Drug Candidates

The epidemic of 2019-novel-coronavirus (2019-nCoV) pneumonia and insights for emerging infectious diseases in the future

Nelfinavir was predicted to be a potential inhibitor of 2019-nCov main protease by an integrative approach combining homology modelling, molecular docking and binding free energy calculation, bioRxiv

Synergistic antiviral effect of Galanthus nivalis agglutinin and nelfinavir against feline coronavirus

HIV protease inhibitor nelfinavir inhibits replication of SARS-associated coronavirus

The significance of chirality in drug design and development

Chiral Toxicology: It's the Same Thing…Only Different

Chirality as an Important Factor for the Development of New Antiepileptic Drugs

Intellectual property and chirality of drugs

Putting chirality to work: the strategy of chiral switches

The predicated demise of racemic new molecular entities is an exaggeration

The market of chiral drugs: Chiral switches versus de novo enantiomerically pure compounds

Chiral pharmaceuticals: A review on their environmental occurrence and fate processes

Density-functional thermochemistry. III. The role of exact exchange

A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model

Gaussian 03, Revision C.02

MolDock: A New Technique for High-Accuracy Molecular Docking

Simmerling, ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

Development and testing of a general amber force field

Comparison of simple potential functions for simulating liquid water

PTRAJ and CPPTRAJ: Software for Processing and Analysis of Molecular Dynamics Trajectory Data

Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models

The MM/PBSA and MM/GBSA methods to estimate ligandbinding affinities

MMPBSA.py: An Efficient Program for End-State Free Energy Calculations

LigPlot+: Multiple Ligand-Protein Interaction Diagrams for Drug Discovery

LIGPLOT: a program to generate schematic diagrams of protein-ligand interactions

Validation of molecular docking programs for virtual screening against dihydropteroate synthase

In Silico Exploration of the Molecular Mechanism of Clinically Oriented Drugs for Possibly Inhibiting SARS-CoV-2's Main Protease

Multidrug treatment with nelfinavir and cepharanthine against COVID-19

Unrevealing sequence and structural features of novel coronavirus using in silico approaches: The main protease as molecular target

Structural basis of SARS-CoV-2 3CLpro and anti-COVID-19 drug discovery from medicinal plants

Computational Design of ACE2-Based Peptide Inhibitors of SARS-CoV-2

Structural Analysis of COVID-19 Main Protease and its Interaction with the Inhibitor N3