Molecular Docking and Dynamics Identifies Potential Repurposed Drug Candidates for COVID-19 Studies

1. Introduction

The novel coronavirus disease 19 (COVID-19) and the post COVID-19 condition called long COVID have affected lives and livelihoods across the world. As of 1 January 2023, we had over 664 million confirmed COVID-19 cases and 6.62 million deaths recorded globally [1] and an estimated 20 million excess deaths obtained by extrapolating the findings published recently in The Lancet [2]. According to the WHO, circa 144 million people suffered from long COVID by the end of 2021 when the recorded cases stood at around 278 million [3, 4], therefore we are likely to have over 330 million long COVID patients worldwide already. The recent resurgence of COVID-19 in China is predicted to result in up to 1 million deaths over the next few months [5].

COVID-19 is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which is part of the RNA virus family Coronaviridae and sub-family Coronavirinae. This sub-family includes four genera of viruses alpha coronaviruses, beta coronaviruses, gamma coronaviruses, and delta coronaviruses, which have been grouped primarily based on serology and phylogenetic clustering [6]. SARS-CoV-2 possesses four structural proteins (S, E, M, and N as defined below) and sixteen non-structural proteins (NSP1−16). The capsid formed from the nucleocapsid protein (N) is present outside the positive-sense, single-stranded RNA genome (+ssRNA) and the genome is further covered by an envelope that is mainly related to three structural proteins: membrane protein (M), spike protein (S), and envelope protein (E). Spike glycoproteins (S protein) form homotrimers that protrude from the host-derived viral envelope and furnish specificity for cellular entry receptors[7].

Vaccines have been useful to restore a level of normalcy but are unable to stop the transmission of the virus and its variants at population levels, therefore it is important to develop affordable therapies that are safe, effective and available, especially to low- and lower-middle income countries that have suffered from insufficient vaccination [8]. Various anti-viral drugs had been developed such as remdesivir, molnupiravir, favipiravir, nirmatrelvir, etc. but none of these meet all four requirements mentioned above. Remdesivir is a prodrug is available as a parenteral [9] while molnupiravir was not shown to have a good clinical outcome from the clinical trials data [10, 11]. Favipiravir was found to be not effective against disease progression [12] and nirmatrelvir, a medication developed by Pfizer is very potent in combination with ritonavir, but has a high price and low manufacturing capability to cater to global demand [13, 14].

Unfortunately, the mean and median capitalized research and development investment of getting a new drug to market were respectively US$1335.9 million and US$985.3 million for the period 2008-19 [15], therefore repurposing drugs approved for other indications to combat COVID-19 and long COVID is a worthwhile strategy [16]. This approach looks for new uses for approved or investigational drugs that are outside the scope of the original medical indication, and can dramatically shrink the costs and timeline involved. However, a major challenge is that there are 7,817 candidates in Compounds Australia open drug collection, so we recently developed a database to down-select the top candidates in a systematic manner [17, 18]. To validate our approach experimentally, 12 of the top 214 drugs, along with the current standards of care (remdesivir, molnupiravir and nirmatrelvir/ritonavir), were recently subjected to in vitro evaluation of anti-viral efficacy against SARS-CoV-2 Delta and Omicron variants of concern [19]. To decide which additional candidates from the list of 214 drugs should be evaluated next, this paper reports an in silico approach consisting of molecular docking studies followed by molecular dynamic studies.

2. Materials and Methodology

Jain et al. [18] have developed the open access, user-friendly ‘CoviRx’ platform to display each of the 7,817 drugs along with physical and chemical properties, original indication, available data from multiple assays, COVID-19 clinical trials, any red flags such as pregnancy concerns, contraindications, etc. and drugs that are similar to the query on the basis of the Tanimoto coefficient [20]. This work made use of CoviRx.org along with the pharmacological down-selection methodology described by MacRaild et al. [18] involving a combination of 11 filters such as approval status, assay data against SARS-CoV-2, pharmacokinetic, pharmacodynamic and toxicity profiles. In addition, the pregnancy and pan-assay interference compounds (PAINS) categories were included to rule out drugs unsafe for pregnant women and those likely to produce false positive results. In this study, the 214 drugs which passed all the filters were subjected to virtual screening against various targets of SARS-CoV-2 followed by molecular dynamic studies of the best scoring ligands against each target (work flow depicted in Figure 1) [21]. The molecular targets include Spike receptor binding domain (RBD), Nucleocapsid protein RNA binding domain (NPRBD), and selected NSPs such as RNA-dependent RNA polymerase (RdRp), Papain-like protease (PL ^pro), 3CL protease/main protease (3CL ^pro/ M ^pro), Helicase, and Exoribonuclease domain (Table 1) [7].

Protein Preparation

The virus-based targets with good resolution were downloaded from the protein data bank (PDB) (https://www.rcsb.org/). These proteins were prepared by adding hydrogens, refining the protein by optimization using PROPKA at pH 7.0, removing water molecules beyond 3.0Å, followed by structural minimization. PROPKA predicts the pKa values of ionizable groups in proteins and protein-ligand complexes based on the 3D structure. This process was performed using protein preparation wizard of Maestro 22 (Schrodinger, LLC, New York, NY, 2022) [22]. The details of virus-based protein targets and their respective PDB IDs are shown in Table 1.

Table 1. Targets used for virtual screening.

S. No	Target protein	PDB ID	Resolution
1.	Spike Receptor Binding Domain	6M0J	2.45 Å
2.	Nucleocapsid protein RNA Binding domain	6VYO	2.2 Å
3.	Papain Like protease (NSP3)	6W9C	2.7 Å
4.	Main protease (3CL protease) (NSP5)	6W63	2.1 Å
5.	RNA dependent RNA polymerase (NSP12)	6M71	2.9 Å
6.	Helicase (NSP13)	7NNG	2.38 Å
7.	Exoribonuclease domain (NSP14)	7R2V	2.53 Å

Table 1. Targets used for virtual screening.

S. No	Target protein	PDB ID	Resolution
1.	Spike Receptor Binding Domain	6M0J	2.45 Å
2.	Nucleocapsid protein RNA Binding domain	6VYO	2.2 Å
3.	Papain Like protease (NSP3)	6W9C	2.7 Å
4.	Main protease (3CL protease) (NSP5)	6W63	2.1 Å
5.	RNA dependent RNA polymerase (NSP12)	6M71	2.9 Å
6.	Helicase (NSP13)	7NNG	2.38 Å
7.	Exoribonuclease domain (NSP14)	7R2V	2.53 Å

Ligand Preparation

The short-listed drugs (214) were downloaded from Pubchem database in “sdf” format (https://pubchem.ncbi.nlm.nih.gov/). These drugs were imported in Maestro, and LigPrep module was used to prepare the ligands by applying optimize potentials for liquid simulation (OPLS4) forcefield. The maximum ligand size was set to 500 atoms to remove larger molecules as they do not obey “Lipinski rule of five”. In addition, ionization was set to neutralize using the Epik module of Schrodinger suite.

Grid Generation

Receptor grid generation wizard was used to select the native position of the co-crystal ligands to validate the docking protocol. If a protein which do not possess co-crystal ligand, then an extensive literature search was conducted to look for the co-ordinates to be given as input for x, y and z co-ordinates to the grid wizard. On the other hand, for the spike protein (6M0J), as it does not possess co-crystal ligand, the key interacting residues with Angiotensin converting enzyme 2 (ACE2) were considered as the center of the grid.

Table 2. Grid co-ordinates of the target proteins with no co-crystal ligand.

S. No	PDB ID	X-axis	Y-axis	Z-axis
1.	6M71 [23]	120.535	116.572	140.185
2.	6W9C [24]	-46.0103	14.2997	29.9464
3.	6VYO [25]	-19.843	10.273	-6.002

Table 2. Grid co-ordinates of the target proteins with no co-crystal ligand.

S. No	PDB ID	X-axis	Y-axis	Z-axis
1.	6M71 [23]	120.535	116.572	140.185
2.	6W9C [24]	-46.0103	14.2997	29.9464
3.	6VYO [25]	-19.843	10.273	-6.002

Molecular Docking Studies

The ligands prepared using ligprep module, the minimized protein and the generated grid files were used for performing flexible ligand docking studies in extra precision (XP) mode by adding Epik state penalties. The number of poses to be included and to write were set to ten and one respectively. All other default parameters were retained.

Molecular Dynamics (MD)

The workstation used to perform the dynamic studies was of Tyrone, equipped with 120GB HDD and 12GB of NVIDIA graphic card. The solvent model used was of transferable intermolecular potential 3P (TIP3P) and the temperature was set to 310.15K. Firstly, the protein was preprocessed to identify the problems if any, and rectified by adding hydrogens and filling the missing residues using prime module of the software. Further, water molecules beyond 5 Å were removed and at pH of 7.0 ± 2.0, heteroatoms state was generated using the software’s Epik module [26]. Secondly, system builder wizard was run by choosing the TIP3P solvent model using OPLS4 force filed and generating an orthorhombic box with dimensions of 10x10x10 Å [27, 28]. The number of ions required to neutralize the system was calculated and added. In addition, salt with concentration of 0.15M was checked in. Thirdly, the output of system builder was used to carry out minimization step for the entire system of complex and solvent for 100ps. Lastly, MD simulation was run using isothermal-isobaric ensemble with constant number of particles (NPT) mode for 100ns.

3. CL Protease/Main Protease

3CL protease, a C30 endopeptidase is found in almost all coronavirus variants and is involved in cleaving polyprotein at different conserved sites [35]. It possesses three major domains (D1, D2 and D3) and the active site remains between D1 and D2 [36]. From our docking analysis, a total of 14 drugs were found to have docking score close to the co-crystal ligand (-6.89 kcal/mol). Among these, two drugs were found to have docking score > -9.0 kcal/mol. These two drugs are treprostinil, a prostacyclin vasodilator, and ETC-1002 (bempedoic acid), an anti-hyperlipidemic agent, with docking scores of -10.76 kcal/mol and -9.65 kcal/mol, respectively. Since, treprostinil had the better score, we analyzed its binding pattern further. Treprostinil formed a total of four H-bonds each with Cys44 and Glu166, and two interactions each with Thr190 and Gln192. The hydroxyl group at 3^rd position on the octyl side chain interacted with the carbonyl group of Thr190 and amino group of Gly192. In addition, the hydroxyl group at 2^nd position of the cyclopentanaphthalene interacted with Glu166 and the carboxylic terminal of treprostinil formed H-bond with Cys44 (Figure 6a).

To understand the behavior of treprostinil, the docked pose of co-crystal ligand as well as treprostinil along with the protein was taken forward for molecular dynamics. During the simulation process, the fluctuations were low as the RMSD was below 3Å. The RMSD of the protein corresponding to treprostinil is relatively less compared to the co-crystal, indicating greater stability of the complex. Further, the ligand RMSD plot shows that co-crystal ligand had less RMSD than treprostinil until 50ns, however, huge deviations was observed until the end of the simulation. Treprostinil had RMSD ranging between 2Å-4Å until 35ns, beyond which the RMSD increased to 4Å from 37ns and attained equilibrium until the end of the study (Figure 6b). The residues that predominantly interacted with the ligand during the simulation process include Asp187, Gln189, Thr190 and Gln192. Additionally, these active site residues have fewer fluctuations during the simulation time as their random mean square fluctuation (RMSF) is below 2Å (Figure 6c). Among these residues, Asp187, Thr190 contributed for >40% of interaction while Glu166, Arg188, Gln189 and Gln192 contributed for interactions ranging between 10% and 20% (See supplementary file 3.4). The confirmations between 34ns-37ns were also captured as the ligand had slight increased RMSD during this interval, and it was observed that the drug retained within the active site. (Figure 6d). As treprostinil has far better docking score than the co-crystal ligand, and the complex is more stable, it can be considered a potential candidate to target main protease experimentally.

RNA Dependent RNA Polymerase (RdRp)

RNA dependent RNA polymerases are highly conserved proteins and catalyze the replication process of an RNA template. Therefore, inhibiting this target could help in preventing the replication of the viral genetic material, thereby reducing the viral load [37]. When we performed the docking studies with the shortlisted ligands, alendronate, a bisphosphonate drug, had the least binding score of -7.86 kcal/mol, better than the control drug remdesivir has lower score of -3.27 kcal/mol. Alendronate was able to form a total of eight interactions each with Asp452, Arg553, Arg555, two with Asp623, and three with Thr556. All the interactions are due to the functional groups present on the bisphosphonate molecule. On the other hand, remdesivir had only one H-bond interaction with each of Asp618, Pro620 and Cys622. Tyr455 also formed pi-pi stacking interaction with the drug due to aromaticity (Figure 7a).

MD studies of the alendronate complex with RdRp reveals a stable RMSD plot. Although it varied from its docking pose at 1ns, the interacting residue was still in contact with alendronate resulting in escalated RMSD to 6Å and greater stability with low 4Å - 6Å fluctuations until 80ns (Figure 7b). Asp623 and Asp760 had constant interactions with alendronate while Arg553, Arg555, Asp618, Lys621 and Cys622 had partial interactions. Beyond 80ns, the ligand further changed its conformation resulting in higher RMSD. The ligand for the last 20ns showed partial interactions with Arg553, Lys621, and Cyc622 while multiple strong interactions were seen with Asp618 and Asp623. The interactions with Asp760 slightly decreased during the end of the study. RMSF plot shows minor fluctuation in the active site residues (<2.5Å), indicating protein stability (Figure 7c). The LP contact plot shows Asp623 and Asp760 contributed 87% and 77% respectively, to the total interactions. Contributions of other residues such as Arg553 (36%), Arg555 (22%), Asp618 (21%), Lys621 (34%) and Cys622 (24%) were moderate (c.f. supplementary file 3.5, alendronate). The time frame analysis also shows the interaction of Asp623 and Asp760 during the fluctuating intervals (Figure 7d).

To compare our results with the control drug, MD analysis was carried out. In our analysis, remdesivir was stable for the initial 10ns with RMSD 4Å - 4.5Å, followed by fluctuation, especially at 20ns - 45ns up to 9Å, then a drop in RMSD of 6Å-7Å during 45ns-55ns, and equilibrium for the rest of the simulation (Figure 7b). The protein was also stable as the interacting residues were within the RMSF of 2.0Å (Figure 7c). The LP contact plot shows minor contributions of the residues during the simulation study, unlike alendronate complex (c.f. supplementary file 3.5, remdesivir). While the key interacting residues of RdRp were present in the MD study of remdesivir, the interactions with alendronate were much stronger. The residue that participated in the direct interactions with remdesivir include Ala550 (20%), Arg555 (11%), Asp623 (13%) and Asp760 (20%) with moderate to low contributions. Alendronate had multiple stronger interactions with the active site and the drug remained within the active site during the entire simulation.

Figure 7. (a) Docking image of alendronate and remdesivir in the active site of RdRp (PDB_id: 6M71) (b) RMSD plot showing alendronate and remdesivir (c) RMSF plot showing fluctuations in protein residues (d) time frame analysis of alendronate showing its poses a at different intervals (e) time frame analysis of remdesivir showing its poses at different time intervals.

Figure 7. (a) Docking image of alendronate and remdesivir in the active site of RdRp (PDB_id: 6M71) (b) RMSD plot showing alendronate and remdesivir (c) RMSF plot showing fluctuations in protein residues (d) time frame analysis of alendronate showing its poses a at different intervals (e) time frame analysis of remdesivir showing its poses at different time intervals.

Helicase

Helicases are the enzymes responsible for unwinding of helical structure of genetic material producing ssDNA/ssRNA that is essential for 3R’s (replication, recombination and repair) [38, 39]. The unwinding process occurs by breaking the H-bonds sequentially between two strands of the double helix and this whole process of bond breaking is catalyzed by ATPase that consumes a molecule of ATP and produces sufficient energy to break the bond [40]. Hence, preventing the bond breaking could halt the unwinding process and thus the 3R’s. As these enzymes are highly conserved across all coronaviruses, it is an interesting target to block the viral replication [41]. Upon docking of the shortlisted 214 FDA/TGA approved drugs, the molecule with best binding score was found to be cromolyn (-5.87 kcal/mol) which is slight better than the co-crystal ligand (-4.02 kcal/mol). The residue that participated in interactions with cromolyn were Lys320, Arg443 and Gly532, while only Lys320 was involved in interacting with the co-crystal ligand (Figure 8a).

RMSD plot from our MD studies reveal initial fluctuations followed by stability. Although Glu261, Ser289, Thr440, Arg442, and Glu540 were absent in the docked complex, they were found to interact with the ligand during the first 5ns of MD studies leading to RMSD of 6Å. Between 10ns -35ns the ligand fluctuated to 10Å after which it attained equilibrium (Figure 8b). During this, the protein active site residues had minor fluctuations (<2Å) as shown in the RMSF plot (Figure 8c). Beyond 35ns, Glu261 contributed significantly (64%) and other residues moderately (15%-20% for Asp260, Asn268, Phe291, Phe437, Thr440, Lys462 and His464) to ligand bonding. Two pi-cation interactions were also seen with His290 (10%) and His464 (11%) (see supplementary file 3.6, cromolyn). Despite fluctuations, the ligand never left the active site during the entire simulation (c.f. time frame analysis, Figure 8d). When we performed the MD studies to the co-crystal ligand, the RMSD trend was similar to that of cromolyn as RMSD escalated between 7Å-8Å for the first half of the simulation time and then declined to 6Å (Figure 8c). All the interacting residues with the co-crystal ligand were within the active site region and the major contributing residues include His290 (33%), Lys320 (76%), Tyr324 (30%), Arg443 (94%) and Lys569 (15%). Since cromolyn and the co-crystal ligand have similar behaviour during the simulation, both molecules could have similar potency.

Exoribonuclease Domain

NSP14 has two terminals with different activities; the N-terminal has 3’-5’ exoribonuclease activity, while the C-terminal has N7-methyl transferase (MTase) activity [42]. While the C-terminal preserves the stability of the viral RNA, the N-terminal role is crucial for maintaining the RNA genome [43]. To exhibit MTase activity, NSP14 uses S-Adenosyl methionine (SAM) as a methyl donor that causes methylation on the Nitrogen at the 7th position of guanine (5’ end), which is necessary for mRNA stability as well as the translation process in human cells [44]. Therefore, by introducing a molecule with more affinity to the active site residues than the SAM, it is possible to block the activity of MTase. When we performed docking studies with this protein (PDB: 7R2V), three drugs had docking score closer to the co-crystal (-11.75 kcal/mol). These drugs were cromolyn (-11.95 kcal/mol), allopurinol (-11.60 kcal/mol) and inosine (-11.59 kcal/mol). Since cromolyn had a better score, we further analyzed it in terms of interactions. SAH has interactions with Arg310, Asp352, Ala353, Tyr368 and Trp385. Unlike the SAH, cromolyn had different interacting residues such Tyr368, Asn386 and Asn388. However, the binding site residues were the same in all three drugs (Figure 9a). As cromolyn had the best score, MD studies were conducted to analyze its behavior and compare its binding properties with that of the co-crystal ligand.

MD analysis of cromolyn shows RMSD ranging between 1.5Å - 6.4Å, and interactions with Asn306, Leu366, Tyr368, Asn388 and Phe426 until 12ns. After 16ns, water-mediated interactions were mainly observed and the RMSD was between 3Å-6Å. (Figure 9b). RMSF plot shows that active site residues in cromolyn complex had similar pattern to that of the co-crystal (Figure 9c). During the simulation, key contributing residues were Asn306 (62%), Trp385 (34%), Tyr420 (58%) and Phe426 (24%) (c.f. supplementary file 3.7, cromolyn). Time frame conformations of the ligand were captured to see the binding pose at different fluctuating intervals. From Figure 9d, it is observed that the change in the conformation was due to cromolyn’s rotatable bonds. When we compared the results with co-crystal ligand, the RMSD reveals fluctuations for the initial 5ns (RMSD 2.4Å - 5.6Å) which started to decline and attain equilibrium (fluctuations between 2.5Å-3Å) towards the end of the simulation (Figure 9b). The residue that contributed in the study include Gln313, Asp331, Ile332, Asp352, Tyr368, Trp385, Asn386 and Asn388. Among these, Asp352 (92%), Tyr368 (98%), Trp385 (74%) and Asn388 (41%) are especially important (c.f. supplementary file 3.7, co-crystal).

4. Conclusions and Future Studies

From the molecular docking and molecular dynamic studies, it is clear that the drugs, treprostinil, natamycin, cromolyn and alendronate have shown promising results against various targets of SARS-CoV-2. RMSD plots and time frame analysis have revealed the behaviour of these drugs during the simulation. In case of main protease, treprostinil had a better RMSD plot when compared to the co-crystal ligand, with limited fluctuations. On the other hand, the fluctuations of natamycin and the residues in papain like protease were similar, with interactions maintained until the end of the simulation. Cromolyn also shows interesting results with three different targets (N protein RNA binding domain, NSP13 and NSP14). It maintained constant interactions during the entire molecular dynamic simulation. Alendronate was able to form strong interactions with two targets (spike receptor binding domain and RNA dependent RNA polymerase). Its interaction with RdRp was far better than remdesivir, which is a current standard of care. From this study, we propose that these four drugs (treprostinil, alendronate, natamycin and cromolyn) could have anti-viral effects and could be a good candidate to treat SARS-CoV-2 infections. As treprostinil promotes vasodilation of pulmonary and systemic arterial vasculature, it might have beneficial effects for COVID-19 patients through improving oxygenation and nitric oxide levels in the blood vessels, inhibiting platelet aggregation, and its anti-inflammatory property [45-47]. Our study further shows that this antihypertensive drug has better binding pattern than the co-crystal ligand. Alendronate is an amino-bisphosphonate indicated for osteoporosis, and patients administered with this drug have reportedly low incidences of COVID-19 testing and hospitalisation [48, 49], which could be due to their immunomodulatory effects and active site interactions discerned in our study. Natamycin is an anti-fungal drug used in eye infections and we could find only one publication pertaining to COVID-19; an in silico study has suggested its use in combination with Baloxavir marboxil and RU85053 as a cocktail as they inhibit 3CL^pro, PL ^pro and RdRp [50]. Our simulation study shows significant binding interactions with PL ^pro and better dynamic behaviour, thus strengthening its evidence base for further investigation. Our study shows that cromolyn, a mast cell stabilizer used as an asthma drug, is promising due to better binding score with multiple interactions and similar dynamic behaviour to that of the co-crystal; unsurprisingly, it is currently in phase III COVID-19 clinical trial [51, 52]. Further in vitro and in vivo studies must be carried out to verify our hypothesis and quantify the anti-viral activities of these four drugs. We are currently implementing the approach presented in this paper to our broader shortlist of 1029 repurposed drug candidates, and improving our filtering methodology with better alternatives to the Tanimoto coefficient [20, 53, 54]. We are also extending our approach to the key variants of SARS-CoV-2 [8].