Computational Prediction of the Potential Target of SARS-CoV-2 Inhibitor Plitidepsin via Molecular Docking, Dynamic Simulations and MM-PBSA Calculations.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) replication depends on the interaction between the viral proteins and the human translation machinery. The cytotoxic peptide plitidepsin was found to inhibit CoV-2 up to 90 % at a concentration of 0.88 nM. In vitro studies suggest that this activity may be attributed to the inhibition of the eukaryotic translation elongation factor 1A (eEF1A). However, recent reports raised the potential for other cellular targets which plitidepsin may use to exert its potent antiviral activity. The lack of data about these potential targets represents a major limitation for its structural optimization. This work describes the use of a molecular modeling approach to rationalize the in vitro antiviral activity of plitidepsin and to identify potential cellular targets. The developed protocol involves an initial molecular docking step followed by molecular dynamics and binding free energy calculations. The results reveal the potential for plitidepsin to bind to the active site of the key enzyme SARS-CoV-2 RdRp. The results also highlight the importance of van der Waals interactions for proper binding with the enzyme. We believe that the results presented in this study could provide the grounds for the optimization of plitidepsin analogs as SARS-CoV-2 inhibitors.

Introduction

The world has been in a fierce fight against the SARS‐CoV‐2 (COVID‐19) infection, a pandemic that disrupted our normal life. By 4th November 2021, the virus has infected more than 248 million people and killed nearly 5 million as reported by the World Health Organization (WHO).[ , ] The positive‐sense virus belonging to the family of beta coronavirus is still spreading.[ , ]

Aiming to prioritize an effective therapy against the pandemic, the WHO has performed many clinical trials using existing drugs. Currently, the treatment protocol relies on oxygen therapy, the antiviral remdesivir, dexamethasone and broad spectrum antibiotics.[ , , ] Remdesivir demonstrated promising efficacy in many clinical trials,[ , ] however, it has no significant impact on the survival and infection duration. This created the urge to develop specific cure for sever and fatal COVID‐19 infection. The magnitude of infections promoted the repurposing of several approved drugs seeking with swift control of the morbidity, mortality, and spread of the virus. Repurposing clinically approved drugs has the advantage of eliminating the need for evaluating their safety and pharmacokinetic profiles.

Based on an in vitro study, SARS‐CoV‐2 the eukaryotic translation elongation factor 1A (eEF1A) was identified as a plausible target for potential SARS‐CoV‐2 drug discovery. It was previously identified as one of the host proteins implicated in SARS‐CoV‐2 replication cycle.[ , ] Interestingly, plitidepsin, a cytotoxic peptide was found to inhibit the eEF1A activity. It is hypothesized to have a similar binding mode to both the structurally related and structurally unrelated didemnin B and ternatin‐4, respectively.[ , ] In the work done by David E Gordon et al., 47 existing drug molecules were tested against SARS‐CoV‐2 with many showed substantial antiviral activity in cell cultures. Of these inhibitors, plitidipsin showed an activity of IC90 of 1.76 nM against SARS‐CoV‐2, which is 9 times as potent as ternatin‐4 in the same assay. In addition, in an assay that used the human cell lines (hACE2‐293T), plitidepsin was proven to be more potent, with an IC90 of 0.88 nM (exceeding remdesivir tested in the same cell lines by a factor of 27.5).[ , ]

A study that used the time‐of‐addition assay, in which either remdesivir or plitidepsin was added to hACE2‐293T cells at different time intervals relative to infection, revealed that 20 nM of plitidepsin is able to strongly inhibit the expression of nucleocapsid protein even after 4 h from the infection. The finding highlights the possibility for a cytoplasmic replication‐stage inhibition that may be exerted by plitidepsin. Despite this study, confirmed the therapeutic potential of the host eEF1A, the same study also raised the potential for the existence of another target other than the host eEF1A that plitidepsin may be using to establish its effective anti‐corona viral activity. Interestingly, the overexpression of A399V mutant eEF1A in cells was shown to build resistance against didemnin B, and ternatin‐4, while plitidepsin retained its inhibition activity against cell proliferation. Unexpectedly, upon measuring the activity of plitidepsin against SARS‐CoV‐2 in transfected 293T cells with eEF1A‐A399V or eEF1A‐WT, it was found that the 293T‐A399V but not eEF1A‐WT cell lines is refractory to plitidepsin SARS‐CoV‐2 antiviral activity by a factor of >12 compared to the parental cell lines. This mutation did not have a similar impact on remdesivir inhibition. This observation further supports the argument that the antiviral mechanism of plitidepsin may not be only due to its impact on the eEF1A host protein, otherwise the plitidepsin would lose its activity completely and conferred a similar resistance gained upon using the Ala399→Val (A399V) mutant of eEF1A in cancer cells.

Pant et al, reported that peptide‐like molecules could be an alternative to small molecules as inhibitors to the COVID‐19 Mpro because of their flexibility. This was confirmed through the study performed on didemnin A, B, and C molecules to test their ability to inhibit the attractive anti‐viral target COVID‐19 viral protease (Mpro) which found that these molecules have a relatively high affinity to the active site of Mpro compared to nelfinavir, and diosmin. These results confirm that the anti‐corona viral activity of didemnin could not be attributed only to its effect on eEF1A confirming also our hypothesis.

Adding to the above discussion, the similar effect shown by plitidepsin to remidesivir in the reduction RNA genomic content at 8 and 12 h after infection and the additive effect shown by plitidepsin with remdesivir using the Synergy finder software, suggests the probability of the binding to the similar RNA‐dependent‐RNA polymerase site of remdesivir. Also, the much greater effect of plitidepsin on the accumulation of the N subgenomic RNA expression after 4 h of infection and maintained its significant effect throughout compared to remdesivir having no effect on N subgenomic RNA at 4 h, indicated that plitidepsin has one or more antiviral mechanism besides the eEF1A host protein inhibition.

Collectively, plitidepsin is suggested to be a promising candidate for COVID‐19. The magnitude of efficacy shown by plitidepsin and the observations stated here prompted our group to elucidate the probability of plitidepsin to act on another molecular target that affects the viral replication (in addition to the eEF1A host protein).

The SARS‐COV‐2 life cycle is well studied, and many enzymes were considered promising targets for the discovery of specific cure for the infection. Those targets are namely, RNA‐dependent‐RNA polymerase (RdRp), 3C‐like protease (3CLpro), papin like protease (PLPro), Helicase and 2′‐ O ‐methyltransferase. Each of the mentioned enzymes plays a crucial role in the virus life cycle. Determining if any of these enzymes can possibly be inhibited by plitidepsin could open a new era in the discovery of novel specific inhibitors with enhanced pharmacokinetic and pharmacodynamics over plitidepsin. Accordingly, in the current work we report the application of molecular modeling techniques into the identification of plitidepsin potential target in SARS‐COV‐2 inhibition.

Material and Methods

2.1. Molecular Docking

The x‐ray crystallographic structures of the SARS‐COV‐2 essential targets were retrieved from the PDB (protein data bank) using the following PDB entries 7BV2, 6Y2G, 6WKQ, 5RLW and 7JRN for SARS‐COV‐2 RNA‐dependent‐RNA polymerase (RdRp), main protease, methyl transferase, helicase and papin‐like protease, respectively. The binding sites dimensions were determined from the binding of each co‐crystalized ligand with its corresponding target. MOE 2019 was implemented in the entire docking studies and the results were visualized by Discovery Studio visualizer. All the potential targets were prepared according to the default settings, by removing crystallization water, adding hydrogens and minimization using AMBER:10 EHT force field. A pose retrieval step for each of the co‐crystalized ligand was conducted at the beginning. Plitidepsin was docked into the pre‐determined active site of the five selected targets. Finally, the most potential target for plitidepsin was selected based on docking score results and was carried on for molecular dynamics and MM‐PBSA calculations.

2.2. Molecular Dynamics and MM‐PBSA Calculations

GROMACS 5.1.1 software,[ , ] was employed to perform all the molecular dynamics (MD) simulations in the current work. To support the outcomes from SARS‐CoV‐2 RdRp docking study, three MD simulation experiments were performed. Two experiments were performed on the enzyme bound to the RNA in complex with either plitidepsin or the crystal reference remdesivir triphosphate, while the third one was conducted on the enzyme bound only to RNA. The same protocol applied from our previous work was typically used in the current study.[ , , , ] The enzyme and ligands topologies were first generated and joined. Then, solvation of the three systems was achieved using single point charge (SPC) water model. Neutralization of the three systems was achieved by adding suitable sodium and chloride counter‐ions. To get energy minimized structures, all the three systems were subjected to not more than 50,000 steps of steepest descent minimization algorithm and a cutoff<10.0 kJ/mole using GROMOS96 43a1 force field. Two consecutive equilibration process was conducted starting with NVT ensemble for 2 ns followed by NPT ensemble for 8 ns. NVT ensemble maintains particles numbers, volume and temperature (310 K) at constant values. Similarly, NPT ensemble maintains particles numbers, pressure and temperature at constant level. The long range electrostatics were described using the standard Particle Mesh Ewald (PME) method. All the equilibrated structures were involved in a production stage for 150 ns using a time step of 2 fs, and saving the average structural coordinates for each 10 ps. Finally, the root means square deviation (RMSD) for the C alpha of the entire system was calculated from the saved trajectories retrieved from the production stage. Also, the distances of the formed hydrogen bonds between the SARS‐CoV‐2 RdRp and plitidepsin were monitored throughout the entire simulation using in‐built GROMACS commands.

2.3. MM‐PBSA Calculation of Binding Free Energy

The calculations of the binding free energy were performed using following equation:

G(Complex) refers to the free energy of the corresponding entity (total, complex, receptor or ligand). The total free energy for each mentioned entity (complex, receptor and ligand) was calculated from all the generated trajectories using the g_mmpbsa. Every type of the calculated energies was provided along with the values of standard deviations before summing them together to yield the average total free energy of each entity. Finally, the binding free energy was calculated by applying the mentioned equation. The calculations were applied for both the two complexes of SARS‐CoV‐2 polymerase‐Plitidepsin and SARS‐CoV‐2 polymerase‐Remdesivir.

The g_mmpbsa package computes four types of energies, including Vander Waal's, electrostatic interaction energies polar energy of solvation and non‐polar solvation energy (SASA model).

The binding free energy for every component (ligand, receptor or complex) was computed by summing all the mentioned four types of energies. At last, the binding free energy was calculated by subtracting the total free energy of both the receptor and ligand from the total free energy of the complex. generally, the negative values of the binding free energy indicate a stable binding and the lesser the free energy values, the more favorable the binding between the ligand and the target.

Results and Discussion

3.1. Docking Studies

In order to predict the potential target for plitidepsin as SARS‐CoV‐2 inhibitor, we performed a docking scan across four key enzymes within the SARS‐CoV‐2 life cycle: RNA dependent RNA polymerase (RdRp), 3C‐like protease (3CLpro), papin like protease (PLPro), Helicase and 2′‐O‐methyltransferase. An initial docking method validation step was performed by re‐docking each native ligand into its corresponding enzyme which in addition provides basis for comparing plitidepsin score with the docking score for each co‐crystalized reference.

Generally, the methodology used produced docked structures with acceptable RMSD values (0.65–1.2) between the co‐crystalized ligand and the docked poses (Supporting Information, Table S1). Using the validated docking method, plitidepsin has shown good scores against all the tested targets and achieved docking scores of −9.7, −7.2, −7.4, −6.5 and −6.0 Kcal/mole for RdRp, SARS‐CoV‐2 main protease, 2′‐ O ‐methyltransferase, SARS‐CoV‐2 helicase and SARS‐CoV‐2 papain‐like protease, respectively. Plitidepsin was able to achieve comparable docking scores with that of the co‐crystalized ligands in all the targets (Supporting Information, Table S1). The results predict high affinity of plitidepsin towards SARS‐CoV‐2 RdRp (−9.7Kcal/mol), even higher than the calculated affinity for remdesivir (−8.2 kcal/mol).

As depicted from Figure  1, plitidepsin was able to engage in a variety of hydrogen bonding as well as hydrophobic interactions at the binding site. It formed many 10 hydrogen bond interactions with the key residues Arg836, Gln815, Lys593, Arg555, Ser759 and Asn691). In addition, plitidepsin engaged in many hydrophobic interactions with LYS551, ARG555, CYS622, LYS593, TRP598, ALA688, LEU758, PRO832 and CYS813. Comparing the binding mode of plitidepsin with that of remdesivir, reveals the superiority of plitidepsin over remdesivir in having multiple interactions with the active site residues Figure  1. Unlike SARS‐COV‐2 RdRp, plitidepsin had relatively little interactions and binding affinities with the other targets. The docking results are consistent with the expected fitting of the large molecule plitidepsin into the wider pocket of SARS‐COV‐2 RdRp as opposed to the relatively narrow binding sites of the other tested targets.

Figure 1

(A) Binding mode of plitidepsin with SARS‐COV‐2 RNA‐dependent‐RNA polymerase (B) Binding mode of remdesivir with SARS‐COV‐2 RNA‐dependent‐RNA polymerase.

3.2. Molecular Dynamics Calculations

Molecular dynamics approach was used to further validate the potential target identified from the docking calculations and to test the effect of the bound ligand on the flexibility of SARS‐COV‐2 RdRp enzyme. MD has proven its worth in many drug discovery studies where running such calculations provides higher level of accuracy especially when studying the flexible nature of macromolecules. Moreover, MD has the power to estimate the binding stability between of ligand‐target complexes. The MD experiment was performed using the free SARS‐COV‐2 RdRp and for the same enzyme in complex with plitidepsin and remdesivir separately.

3.3. RMSD Analysis and Hydrogen Bond Monitoring

The SARS‐COV‐2 RdRp enzyme relies on the flexibility of its wide active site to conduct its function though accommodating both the original template and the generated replicate during viral genome replication. Thus, this degree of active site dynamicity is crucial for the enzyme to exert its activity. Part of conducting our MD simulations was to examine the potential effect of the bound ligand on the overall flexibility of the enzyme. As expected, the RMSD for the C‐alpha for all the residues in free enzyme reached 4 Å in the free enzyme highlighting the dynamicity of the enzyme, Figure  2. Interestingly, the RMSD for the C‐alpha for all the residues and ligands reached 1.1 and 1.9 Å for the enzyme in complex with plitidepsin and remdesivir, respectively Figure  2. We propose that the imposed enzyme rigidly after binding to either ligand may contribute to the activity of these ligands. Here again, plitidepsin has shown superior result over remdesivir as an inhibitor for the SARS‐COV‐2 RdRp.

Figure 2

RMSD for all the trajectories of the MD simulations.

We also took the advantage of the GROMACS ability to measure the stability of the formed hydrogen bond interactions to validate the predicted binding mode of plitidepsin with its potential target. After running specific built‐in commands in GROMACS, the average distance, as well as the standard error for each hydrogen bond formed between plitidepsin with its potential target was calculated. A hydrogen bond is considered valid when the distance between the hydrogen bond donor and acceptor is maintained at less than 3.5 Å. As detailed in Table  1, all the formed hydrogen bonds kept acceptable lengths throughout the entire MD simulations indicating their stability and reliability. Also as depicted from Figures  3 and 4, plitidepsin was able to form various stable hydrogen bonds with its target through the entire MD simulations.

Table 1

The average distances of all the hydrogen bonds formed between plitidepsin and SARS‐COV‐2 RNA‐dependent‐RNA polymerase through the entire 150 ns MD simulation.

Hydrogen bond name	Average distance (Å) +/− SD
Hydrogen bond with Arg836	1.95 +/− 0.07
Hydrogen bond with Gln815	2.66 +/− 0.13
Hydrogen bond with Lys593	2.71 +/− 0.36
Hydrogen bond with Lys593	2.51 +/− 0.17
Hydrogen bond with Arg555	2.52 +/− 0.08
Hydrogen bond with Arg555	2.25 +/− 0.05
Hydrogen bond with Arg555	2.45 +/− 0.05
Hydrogen bond with Arg555	2.95 +/− 0.40
Hydrogen bond with Ser759	2.33 +/− 0.09
Hydrogen bond with Asn691	2.65 +/− 0.09

Figure 3

The percentage of existence for each formed Hydrogen bond.

Figure 4

Number of hydrogen bonds between plitidepsin and SARS‐COV‐2 RNA‐dependent‐RNA polymerase.

3.4. Binding Free Energy Calculations Using MM‐PBSA Approach

The binding free energy for plitidepsin with its potential target SARS‐COV‐2 RdRp was calculated using the MM‐PBSA approach which enables the calculation of binding free energy for all the conformations in the saved trajectories. The g_mmpbsa package generated by Kumari et al., was used to calculate all the MM‐PBSA binding free energy forms (van der Waal energy, Electrostatic energy, Polar solvation energy and SASA energy) for the two complexes of SARS‐COV‐2 RdRp bound plitidepsin and remdesivir. The binding free energy of interaction of plitidepsin and remdesivir with the SARS‐COV‐2 RdRp was decomposed into van der Waals, electrostatic and polar components (Table  2). The calculated binding energy show a significantly higher affinity for plitidepsin compared to remdesivir which is consistent with the results from the docking and MD calculations. Interestingly, the van der Waals interaction energy was calculated to be 2–3x that of the electrostatic interaction energy (Table  2) which highlights the importance of the different alkyl and benzyl moieties of plitidepsin in binding to SARS‐COV‐2 RdRp binding site.

Table 2

Summary of the interaction energies and the binding free energy for both the complexes.

Complex	ΔEbinding (kj/mol)	ΔEElectrostatic (kj/mol)	ΔE Van der Waal (kj/mol)	ΔEpolar solvation (kj/mol)	SASA (kJ/mol)
Plitidepsin	−332.2±18.9	−119.2±15.7	−295.1±22.1	111.6±17.7	−29.5±1.5
Remdesivir	−317.2±18.7	−109.8±14.9	−283.5±23.2	102.9±16.9	−26.8±1.5

ΔE

Conclusion

In the current study, we have used a combination of various molecular modeling techniques in order to rationalize the reported activity of plitidepsin by identifying potential targets among SARS‐COV‐2 key enzymes. Applying a protocol of molecular docking, MD and MM‐PBSA free energy calculations we proposed that plitidepsin can possibly bind to the four tested targets with COV‐2 RNA‐dependent RNA polymerase (RdRp) being the most potential. Analysis of the docking calculations identified potential binding site interactions and results from MD simulation revealed the importance of restraining the enzyme flexibility for its inhibition. In addition, analyzing the data from MM‐PBSA free energy calculations reveals the significance of hydrophobic contacts to the enzyme active site compared to the electrostatic interactions. This emphasizes the importance of the alkyl and benzyl moieties for proper plitidepsin binding. These modeling outcomes described in this article may provide the grounds for a lead optimization program for the introduction of plitidepsin‐like compounds as COV‐2 RdRp inhibitors.

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

Author Contribution Statement

Mahmoud A. El Hassab, Loah R. Hemeda, Zainab M. Elsayed and Mohammed K. Abdel‐Hamid conducted the computational work of the manuscript. Wagdy M. Eldehna, Sara T. Al‐Rashood and Hatem A. Abdel‐Aziz wrote the manuscript and approved the final version.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.