Effect of nelfinavir stereoisomers on coronavirus main protease: Molecular docking, molecular dynamics simulation and MM/GBSA study.

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.

Introduction

The COVID-19 has been a major cause of death in recent days. The worldwide number of deaths reached 431,192 by June 15, 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% identity.

Also, no mutations in these enzymes have been reported [10]. Several researchers have tried to discover new inhibitors for the main protease of the SARS-Cov-2 [[11], [12], [13], [15], [16], [17]].

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.

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–5, respectively.

In this study, all possible nelfinavir stereoisomers are constructed and their binding strength to the protease enzyme (see Fig. 1 and Fig. 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.

Fig. 1

Structure of main protease of SARS-Cov-2 along with its binding cavity.

Method

Optimization of stereoisomers and ligand atomic partial charge determination

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.

Docking

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.

Molecular dynamics simulation

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 300 K for 0.1ns after the minimization process. The density of the water molecules was relaxed for 0.1 ns at a constant pressure of 1 at m and a temperature of 300 K. 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 300 K for 100 ns

Molecular dynamics simulation analysis

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].

Results and discussion

The MVD score algorithm was first validated for the main protease experimental structure (PDB code: 5R81) to ensure that the ligands docked using the virtual docker of molegro represent a valid score and an accurate binding to the receptor. In order to validate the docking protocol, crystallized ligand was extracted from the protein structure and the ligand was docked into the 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 Fig. 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.

Fig. 4

The superposing of docked ligand on co-crystallized ligand of Z1367324110(PDB code: 5R81) for validation of docking protocol. The crystallized and docked ligands are shown in yellow and green colors, respectively. Cresset visualization software was used to plot surface. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Using the above docking protocol, 32-stereoisomers shown in Fig. 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. Fig. 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 Fig. 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 Fig. 4, the position of ligand Z1367324110 is at the anchor site of the main protease.

Fig. 2

Optimized structures of 32 stereoisomers of nelfinavir. Each stereoisomer was named based on label name of R or S on each five stereocenters as shown in scheme 1. The first position in nomenclature is the carbon 1 label, the second the carbon 2 label and so on.

Table 1

Results of molecular docking of 32 stereoisomers of nelfinavir against main protease.

Compound	Diastomeric Code	Moldock Score (kcal/mol)
1	RSRSR	−178.6
2	SRRSR	−174.1
3	SSRRS	−171.1
4	SSSRS	−169.5
5	RRSSR	−169.5
6	SRSSS	−168.5
7	RSRSS	−168.4
8	RSRRR	−167.4
9	RRSRS	−167.4
10	SRSRS	−165.0
11	RRRSR	−164.8
12	SSSSR	−164.6
13	SSRRR	−164.4
14	SSRRS	−163.7
15	SSSSS	−162.9
16	SRRSS	−162.6
17	RRSSS	−162.5
18	RRRRR	−161.9
19	SRSSR	−161.0
20	RSSSR	−160.4
21	RSSRR	−160.3
22	RRRRS	−159.6
23	SSRSR	−158.5
24	RRSRR	−158.0
25	RSSSS	−156.0
26	SRSRR	−155.6
27	RRRSS	−155.1
28	SRRRS	−155.0
29	SRRRR	−153.4
30	SSSRR	−153.3
31	RSRRS	−151.4
32	RSSRS	−147.2

Fig. 5

The position of the anchor site and electrostatic potential mapped to the molecular surface of coronavirus 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 Refs. [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 complex, SSRRS complex and apo protein are 2,367, 2,234, 2127 and 1,567, respectively. This result is supported by the findings of other researchers [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 (see Fig. 3a).

Fig. 3

(a) Root mean square deviation (RMSD), (b) root mean square fluctuation (RMSF), (c) gyration radius analysis for all three complexes along with apo-protein form of main protease. The related plots for RSRSR, SRRSR, SSRRS complexes and apo-protein are shown in brown, violet, green and orange colors, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Residual flexibility of a protein can be considered using RMSF data. Fig. 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. Fig. 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% of the simulation time. Also, Gly143 forms a hydrogen bond with RSRSR for 24% of simulation time. Hie41 is also a residue that is involved in ligand hydrogen bonding within 7% 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% of simulation time (see Fig. 6 ).

Table 2

Hydrogen bond analysis of stereoisomer of RSRSR, SRRSR and SSRRS during the molecular dynamics simulation.

	Acceptor	Donor H	Donor	Frac	AvgDist	AvgAng
RSRSR	Glu166@O	LIG305@H	LIG305@N	0.345	2.701	162.415
	Glu166@O	LIG305@H	LIG305@O	0.057	2.854	142.547
	LIG305@O	Glu166@H	Glu166@N	0.542	2.896	160.017
	LIG305@O	Gly143@H	Gly143@N	0.248	2.867	157.900
	LIG305@O	Hie41@H	Hie41@N	0.079	2.900	152.562
SRRSR	LIG305@O	Ser46@H	Ser46@O	0.572	2.693	163.440
	LIG305@O	Hie41@H	Hie41@N	0.018	2.863	153.069
SSRRS	Glu166@O	LIG305@H	LIG305@O	0.213	2.700	163.084

Fig. 6

Detailed view of the interactions between RSRSR, SRRSR and SSRRS with main protease (left panel). The carton representation of interaction between ligand and main protease (right panel).

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.

Table 3

Binding energy analysis in kcal/mol using MM/GBSA method.

	ΔEvdW	ΔEELE	ΔEGB	ΔESURF	ΔGGAS	ΔGSOL	ΔGBinding
RSRSR	−47.8 ± 0.20	−25.9 ± 0.38	39.7 ± 0.26	−5.8 ± 0.02	−73.7 ± 0.43	33.9 ± 0.25	−39.8 ± 0.29
SRRSR	−35.3 ± 0.24	−20.3 ± 0.38	30.9 ± 0.40	−4.7 ± 0.02	−55.6 ± 0.44	26.2 ± 0.38	−29.5 ± 0.22
SSRRS	−35.6 ± 0.21	−5.7 ± 0.26	17.5 ± 0.23	−4.1 ± 0.02	−41.3 ± 0.33	13.4 ± 0.22	−27.9 ± 0.24

ΔEvdW, van der Waals contribution; ΔEElE, electrostatic contribution; ΔGGAS represents gas-phase energy; ΔGGAS = ΔEvdW+ΔEElE;ΔESURF, nonpolar solvation free energy; ΔEGB, polar solvation free energy; ΔGSOL = ΔESURF + ΔEGB; ΔGBinding = ΔGGAS + ΔGSOL.

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 Fig. 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.

Fig. 7

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

Conclusions

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.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.