Ensemble-based screening of natural products and FDA-approved drugs identified potent inhibitors of SARS-CoV-2 that work with two distinct mechanisms.

The recent outbreak of SARS-CoV-2 is responsible for high morbidity and mortality rate across the globe. This requires an urgent identification of drugs and other interventions to overcome this pandemic. Computational drug repurposing represents an alternative approach to provide a more effective approach in search for COVID-19 drugs. Selected natural product known to have antiviral activities were screened, and based on their hits; a similarity search with FDA approved drugs was performed using computational methods. Obtained drugs from similarity search were assessed for their stability and inhibition against SARS-CoV-2 targets. Diosmin (DB08995) was found to be a promising drug that works with two distinct mechanisms, preventing viral replication and viral fusion into the host cell. Isoquercetin (DB12665) and rutin (DB01698) work by inhibiting viral replication and preventing cell entry, respectively. Our analysis based on molecular dynamics simulation and MM-PBSA binding free energy calculation suggests that diosmin, isoquercetin, rutin and other similar flavone glycosides could serve as SARS-CoV-2 inhibitor, hence an alternative solution to treat COVID-19 upon further clinical validation.

Introduction

The outbreak of the novel Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) pandemic around the globe have called for an urgent research of drugs to treat the disease. Owing to the pressing health response, drug repurposing presents huge potential towards combating diseases. Drug repurposing is a strategy involved in identifying old approved drugs to treat new diseases [1,2]. This process is effective, efficient, and reliable in drug development compared to de novo approaches [[1], [2], [3]]. Such identified drugs have known pharmacokinetics profile and can be tested in clinical trials which fasten the processes of obtaining new drugs [4]. In an efforts to find new drugs to treat COVID-19 several approved antiviral drugs, [2,4,5] antimalarial drugs such as chloroquine and its derivative hydroxychloroquine, [3]; and [6] antiparasitic drugs [7] and antimicrobial agents such as azithromycin [8] just to mention a few, have been repurposed. Although current research on drug therapeutics and vaccines are yet to produce approved drugs or vaccine against SARS-CoV-2, some promising therapeutics such as remdesivir and galidesivir are in clinical trials [9].

Earlier works on drug repurposing against COVID-19 have targeted the main protease [1,10]. Polypharmacology, a process in drug design where on average a drug has 3 targets is highly recommended [1]. Identification of potential approved drugs with less side effects serve as a strong basis for drug repurposing [1]. Thus, analysis of SARS-CoV-2 molecular therapeutic targets forms an important process of identifying drugs to treat this disease [11]. In an attempt to identify new therapeutics against SARS-CoV-2, we resorted to use previously reported antiviral compounds from our research groups [12,13]. Such small molecules could be important hits towards finding new drugs. The selected hits molecules were used for similarity search for further molecular target evaluation. Molecular targets such as 3C-like protease (3CLpro), responsible for viral replication; RNA-dependent RNA polymerase (RdRp), responsible for viral transcription/replication; papain like protease (PLpro), responsible for cleavage of N-terminus of the replicase poly-protein; human angiotensin converting enzyme 2 (ACE2), responsible for regulating cardiovascular physiology in host, and spike receptor binding domain (RBD)-ACE2, responsible for viral attachment before entering host cells were considered important therapeutic targets in treating COVID-19. It is important to highlight that the SARS-CoV-2 main protease (Mpro) is a cysteine protease with an active site catalytic dyad HIS41 and CYS145 located between domain I and II (Fig. S1) [14,15]. The imidazole group of HIS polarizes and activates the SH group of CYS to form a highly nucleophilic CYSS−/HISH+ ion pair that reacts with the substrate [14].

Virtual screening via docking is an effective method and also a widely used approach in screening large libraries of drugs to determine potential drug candidates. To provide computation efficiency and accuracy, several free energy methods with different scoring functions have been developed and used in both hit identification and optimizations processes [1]. Free energy calculation methods such as thermodynamics integration (TI), [16] molecular mechanics Poisson Boltzmann/Generalized Born surface area (MM-PB(GB)SA), [17] free energy perturbation (FEP) [18] and Linear interaction energy (LIE) [19] are efficient methods in structure based drug design. In this study, we employed computational methods to identify new molecules against SARS-CoV-2. First, we screened antiviral natural products previously reported from Tanzanian medicinal plants against infectious bursal disease virus and newcastle disease virus [12,13]. The obtained hits prompted similarity search for FDA approved drugs that we further subjected to computational inhibitory assays against SARS-CoV-2 proteins. The potential molecules for drug repurposing as reported in this work provides an basis for scientists to develop agents against coronaviruses. Further in vitro and in vivo clinical validation of these molecules is necessary to provide useful information on COVID-19 treatment.

Materials and methods

Virtual screening and molecular docking

Ligand preparations

We started by building a small in-house library comprising of 19 natural products with antiviral properties previously obtained from Tanzanian medicinal plants [12,13]. The 2D and 3D structures of the selected natural products were built, then Gasteiger charges [20] were added and all structures energy minimized using universal force field (UFF) [21] employing steepest descents minimization in Avogadro [22]. The optimized geometric structures were converted to pdbqt file format using Open Babel [23] ready for screening by docking with AutoDock Vina [24].

Protein preparations and docking benchmarking

The study took advantage of the timely crystallized SARS-CoV-2 structural targets for drug design. The protein structures were obtained from the RCSB protein data bank (PDB) [25] with the following PDB IDs; 3C-like protease (3CLpro) (6LU7) [26], RdRp (6M71) [27], PLpro (6W9C) [28], nucleocapsid N-terminal RNA binding protein (6M3M) [29], Nsp9 RNA binding protein (6W4B) [30], native human ACE2 (1R42) [31] and spike RBD-ACE2 (6LZG) [1]. Missing hydrogen at pH 7.4 and Gasteiger charges [20] were added to all protein structures and were converted to pdbqt format as required by AutoDock Vina [24] for docking calculations.

Docking calculations were validated using a set of data with binding free energy as previously reported [32]. A set of experimental binding free energy reported in Refs. [33,34] were used to validate the docking scoring function implemented in our work. Docking calculation was performed to each complex and then the lowest binding free energy (more negative value) of each complex was compared with the experimental data and found to be in a reasonable agreement with a correlation value of r 2 = 0.83 (see Table S1). Besides using experimental data, a co-crystalized ligands of 3CLpro (PDB ID: 6LU7) and (PDB ID: 6W63) were extracted from the crystal binding pose and redocked again to the putative binding pocket to assess whether our method could reproduce experimental binding modes. Redocking gave acceptable RMSD values of ≤ 3 Å for all redocked ligands (Fig. S1). RMSD of ≤ 3 Å was found to be more reliable and can reproduce the experimental binding mode [35]. After validating the docking protocol, an in-house library of 19 natural products (see SI) were blindly screened against SARS-CoV-2 targets, the promising hit compounds with lower binding free energy of −7.5 kcal/mol were selected for further experiments and similarity search. Determination of the best binding energy of small molecules was based on the binding free energy of a co-crystallized ligand that was docked and found to be < −7.5 kcal/mol. In this work, we set the criterion that any ligands with binding free energy ≤ −7.5 kcal/mol similar to inhibition constant of Ki = 3.03 μM to be a potential candidate. Similarity search was performed to repurpose approved drugs for SARS-CoV-2 target proteins. All virtual screening docking calculations were done using AutoDock Vina tool [24] employing a python shell script customized for virtual screening.

Similarity search

Compounds obtained from docking experiment (Fig. 1 a) were subjected to similarity search. A web tool, swisssimilarity (http://www.swisssimilarity.ch) which comprises of database with approved drugs (such as DrugBank [36], PubChem [37] etc) was used to perform similarity search [38]. An FP2 fingerprint with a Tanimoto score of 0–1 was used, with 0 depicting no similarity and 1 depicting complete identity. Compounds with higher Tanimoto score of ≥ 0.9 were selected for this study.

Fig. 1

(a) 2D structures of the hit natural product compounds obtained from virtual screening. (b) 2D structures of the FDA approved natural compounds obtained using similarity search, the drugs were similar to ligand 15.

Quantum mechanical calculations

Compounds obtained from similarity search were subjected to quantum mechanical calculations. The selected geometric structures were optimized using the density functional theory (DFT) method with the B3LYP (Becke three-parameter Lee Yang Parr) in conjunction with the 6-311+G(d,p) basis set [39,40]. The calculated vibrational frequencies ascertained that the structures were minima by the absence of imaginary frequencies. All the calculations were performed using Gaussian16 suite [41] and then the compounds were subjected to molecular docking against all targets to investigate their inhibition effects. Best hits from docking calculations were then selected for molecular dynamics (MD) simulation to evaluate their binding stability in the pockets. Molecules obtained from similarity search are shown in Fig. 1b.

Molecular dynamics (MD) simulation

All drug molecules with lower binding free energies obtained from docking screening were subjected to MD simulation to assess their time-lapse stability. MD simulation for all complexes was performed using GROMACS v2018 employing all-atom Optimized Potential for Liquid Simulation (OPLS-AA) force field [42,43]. Protein topologies were generated using GROMACS while ligand topologies were generated using LigParGen tool [44]. All systems were placed at the center of cubic box solvated with TIP4P water model [45]. The initial charges were; −5 for Mpro-DB01698 complex, −5 for Mpro-DB08995 complex, −5 for Mpro-DB12665 complex, −26 for RBD-ACE2-DB08995 complex and −5 for ACE2-DB08995 complex. Depending on the charge of each system, sodium ions were added to neutralize the system. All systems underwent two energy minimization steps to remove bad and unfavourable atom contacts from molecular docking experiments. Initially, the systems were energy minimized by using steepest descent algorithm followed by conjugate gradient algorithm and the systems converged at Fmax < 1000 kJ/mol in 50,000 energy minimization steps. Two equilibration phases with position restraints were carried out to relax the system. First, the systems were equilibrated in an NVT ensemble using the V-rescale thermostat [46] for 100 ps. The systems were further equilibrated in the second phase at an NTP ensemble using Parrinello-Rahman barostat at an isotropic pressure of 1 bar for 500 ps [47]. MD production was carried out without position restraints at an NPT ensemble for up to 100 ns depending on the system. During the MD production both pressure and temperature were maintained using Parrinello-Rahman barostat and V-rescale thermostat at 1 bar and 300 K, respectively. In all systems, periodic boundary conditions (PBC) were applied in all directions, Particle Mesh Ewald (PME) method was used to treat long-range electrostatic interactions with a cutoff distance on 11 Å for both electrostatic and van der Waals interactions [48]. Bond length were constrained using LINCS algorithm [49] and an integration time step of 2 fs was used for all system. An extended analysis was then performed using GROMACS built-in and in-house tools. The distances between protein and ligands was measured from specific groups which formed hydrogen bonds, and the free energy was computed as F = −k Tln(P) [50], where k is the Boltzmann constant, T is the temperature and P is the probability distribution of the distances. Visualization was done using PyMOL, [51] Chimera [52] and VMD [53].

Relaxed complex scheme (RCS)

To get more insight into the effects of protein conformation changes on ligand binding affinity, RCS or ensemble based docking was performed as previously described [4]. After performing an MD simulation of the apo and holo protein, a total of 50 snapshots were evenly selected from the trajectory. Docking calculations were performed on each structure and the binding energies were averaged and presented as mean ± SD. All docking calculations were performed using AutoDock Vina tool [24]. In each ensemble structure, the binding site was specified before docking.

MM-PBSA binding free energy calculations

Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) was calculated using g_mmpbsa tool [17]. The binding free energy (ΔGbind) for protein-ligand interaction is expressed in Equation (1) where, is the average molecular mechanics energy terms in vacuum which includes bonded and non-bonded interactions, is the free energy of solvation which includes Gpolar and Gnon-polar. γ is a coefficient related to surface tension, SASA is the solvent accessible surface area, b is the fitting parameter and TΔS is the entropic contributions to free energy. In g_mmpbsa an entropic contribution is not considered [17]. The binding free energy was calculated using a single trajectory, where a total of 200 snapshots were evenly extracted at a predetermined intervals, from the first 40 ns and last 40 ns of the trajectory. The solvent dielectric constant was 80, solute dielectric constant was 2, γ was 0.0226778 kJ/Mol/Å2 and b was 3.84928 kJ/Mol, the PB equation was solved by using the linear PBsolver.

Results and discussion

Before inferring to any docking calculation, it is very important to understand the limitations of the protocol and tools used. One of the limitations of many docking scoring functions is the ability to reproduce experimental results [54]. Thus, docking methods need to be validated before they can be used in any experiment. One of the most effective ways of validating docking scoring function is to assess its ability to reproduce the binding mode of the small molecule to its receptor determined by root mean square deviation (RMSD). An RMSD value of ≤ 3 Å is an acceptable value in assessing whether the binding mode is reproduced or not [35,55]. In addition to RMSD, comparison to experimental data is another approach used to validate docking methods. In this work, we have used our protocol [4] to validate the docking methods. Redocking of the co-crystalized ligand gave an RMSD value of 2.5 Å which is an acceptable value, this permitted us to continue with the experiments.

After validating the docking methods and identifying the binding sites of each target, ligand virtual screening was carried out against seven SARS-CoV-2 target proteins viz: 3CLpro, PLpro, Nsp9 RNA, Nucleocapside RNA, ACE2, RdRp and spike RBD-ACE2. The selected target proteins represent two broad categories of therapies, those acting to disrupt the virus represented by the first six targets and those blocking the virus from entering the human cell, represented by the last target. Ligand screening was carried out to identify a drug that can disrupt the virus and those that can prevent the virus from entering the host cell. Our in-house library of natural products were screened against each target protein and molecules with binding energies less than −7.5 kcal/mol were considered to have a higher binding affinity with the target protein and were selected for similarity search.

Targeting inhibition of 3CLpro

3CLpro represents an important therapeutic target in drug development to fight COVID-19. This protein is responsible for mediation of non-structural proteins (Nsps) required for the virus life cycle. The protein is cleaved from poly-protein to release Nsps4-Nsp16 [11]. Inhibition of this protein is essential in preventing viral replication and is a promising target for SARS-CoV-2. 3CLpro is made up of three domains: domain I is formed by amino acid residues 8–101, domain II formed by amino acid residues 102–184 and domain III formed by amino acid residues 201–303 and a loop which connects domain II and III formed by amino acid residues 185–200 [11,56]. The binding pocket of the 3CLpro is situated between domain I and II and is composed of the catalytic dyad CYS145 and HIST41 [11,56]. In order to gain more insight into protein plasticity and conformation changes at the active site, MD simulation of the free protein was carried out. RMSD and root mean square fluctuation (RMSF) analysis showed that, although the active site is composed of long loops it has few fluctuations (Fig. S2).

Screening of our in-house library of natural products against this target showed that four ligands: 1, 5, 15 and 17 possessed low binding free energy of ≤ −7.5 kcal/mol (Fig. 2 a). However, ligand 1 showed a higher affinity with the binding free energy of −8.2 kcal/mol corresponding to 0.93 μM. The four molecules were advanced to the next filter which employed the RCS. A RCS is an important filter as it represents protein flexibility which is a limitation of many docking algorithms, which does not accommodate protein flexibility. In this filter, ligands were re-ranked, some ligands which ranked high when docked to the crystal structure in the ensemble structure ranked least (Fig. 2a). Based on the ensemble holo structure, three ligands: 1, 5 and 15 were selected for similarity search to identify similar molecules. Using a Tanimoto index of 0.9, three glycosylated flavonoids (Fig. 1b) similar to ligand 15 were established and subjected to RCS docking to evaluate their binding affinity to the 3CLpro (Fig. 2a). All three drugs showed lower binding free energy compared to the first screened natural products. Binding mode analysis (Fig. S3) showed the three drugs bound to 3CLpro in the same fashion by forming at least two hydrogen bonds.

Fig. 2

(a) Binding free energy of the natural product compounds to 3CLpro (b) Probability distribution of RMSD (p(rmsd)) for both ligand and receptor during the simulation time of 100 ns. (c) RMSF differences of apo and bound 3CLpro. (d–e) Distances between specific atoms in drugs and protein amino acid residues that formed hydrogen bonds with DB12665, (f) Distances between specific atoms in drugs and protein amino acid residues that formed hydrogen bonds with DB01698.

Docking calculation of binding free energies is useful in predicting good and bad binders as well as screening large chemical libraries. However, they cannot provide details on the time-lapse stability and are mostly poor in estimating the binding affinity of the complex. To overcome these limitations, MD simulation was carried out to establish the time-lapse stability of the complex and binding affinity of the ligands which was evaluated using the MM-PBSA method. RMSD, RMSF, distances between drugs and protein, and the number of hydrogen bonds were used to analyse the stability of the complexes. RMSD analysis presented in Figure S4, clearly shows that, the system made by DB12665 exhibited some fluctuation during the first 40 ns and remained stable throughout the simulation time. The system made by DB01698 exhibited fluctuations throughout the simulation time, suggesting unstable interactions. However, the RMSD of DB08995 remained stable though it induced conformational fluctuation to the protein. The RMSD probability distribution for DB01265 exhibited single maximum for both ligand and protein (Fig. 2b) suggesting the presence of stable conformation. However, multiple maxima were observed for the protein and ligand in the system of DB01698 and DB08995, which suggests to have unstable interaction. On a closer look, the RMSD for DB08995 calculated after least square fit showed single peak. This analysis suggests that, although DB08995 induced conformation changes to the protein, it remained stable in its pocket. Atomic residues motion and fluctuations are highly related to protein functions. In order to provide information on how the binding of the drugs brought about dynamic changes and motion which influences their function, RMSF difference between apo and holo structures were analysed and presented in Fig. 2c. The binding of the three drugs clearly shows that they induced different motions to the protein. Interestingly, the binding of the DB08995 induced conformation changes at the regions of residues 100–250, whereas such induced fluctuations are not well pronounced in the complexes of DB01698 and DB12665 except at the region around 160 for DB01698 and 125 and 160 for DB12665. Another residue fluctuations for the later two complexes are observed at the region of 50. Such fluctuations were not observed in the former complex. This clearly shows that the binding of DB08995 induced protein motion and dynamics in a different manner when compared to other complexes which suggest that the drug retained its binding pose and changed protein functions. This analysis further suggests that DB12665 and DB08995 could therefore interfere with coronavirus viral RNA synthesis and replication.

In order to explain the dynamic nature of the ligands inside the active site, hydrogen bonds and their distances were measured over the simulation time. The number of hydrogen bonds (Fig. S6b-d) shows that, DB12665 and DB08995 formed a maximum of 8 and 7 hydrogen bond, throughout the simulation time, respectively. The number of hydrogen bonds for DB01698 decreased to 0 during the 40 ns, which suggests that the drug moved out of the pocket. From 50 ns a maximum of 3 hydrogen bonds were formed. Distances between groups of DB12665 which formed hydrogen bonds with residues of 3CLpro were analysed (Fig. 2d–e). DB12665 exhibited changes in the binding pose from the initial docking pose, and remained stable in the second pose. The system made by DB01698 showed higher fluctuations for both receptor and ligand RMSD suggesting to be unstable. Additional analysis based on the number of hydrogen bonds, distances between protein and ligand atoms and snapshots obtained from the trajectory (Fig. 2f, Fig. S7) indicated that DB01698 moved out the pocket after 40 ns. Further analysis on the distance and binding mode showed that DB08995 retained its binding pose (Fig. S6a).

In order to ascertain the observed stability, binding free energy based on MM-PBSA was calculated to evaluate the affinity of the drugs in their active site. MM-PBSA was calculated on a single trajectory using frames from the first 40 and last 40 ns. Based on the last frames, MM-PBSA analysis showed that complex of DB12665 and DB08995 had higher binding affinity than DB01698. Both van der Waals and electrostatic energy terms highly contributed to the interaction of the protein and drugs (Table 1, Table 2 ).

Table 1

Binding free energy (kJ/mol) for three complexes during the first 40 ns.

Energy	Complexes
	DB1265	DB01698	DB0895
ΔEvdW	−116.9 ± 24.4	−128.3 ± 31.8	−172.4 ± 17.6
ΔEele	−111.1 ± 56.2	−59.3 ± 28.3	−111.9 ± 31.4
ΔEpolar	161.6 ± 42.3	116.9 ± 32.9	157.4 ± 35.5
ΔEnon-polar	−16.6 ± 2.0	−16.5 ± 2.7	−19.3 ± 1.2
ΔEbinding	−82.9 ± 24.1	−87.2 ± 27.1	−146.3 ± 19.6

Table 2

Binding free energy (kJ/mol) for three complexes for the last 40 ns.

		Complexes
Energy	DB1265	DB01698	DB0895
ΔEvdW	−128.2 ± 17.4	−13.2 ± 23.9	−164.9 ± 13.2
ΔEele	−71.6 ± 16.7	−6.1 ± 14.2	−109.6 ± 31.2
ΔEpolar	184.5 ± 14.0	14.0 ± 28.9	165.5 ± 36.2
ΔEnon-polar	−18.6 ± 1.0	−1.8 ± 4.1	−18.9 ± 1.3
ΔEbinding	−33.9 ± 17.9	−7.1 ± 24.3	−127.9 ± 19.5

Per-residue energy decomposition analysis

To examine residue contributions to the binding energy of the three ligands, per-residue energy decomposition was performed. The analysis was done using the last 40 ns trajectories and results are presented in Fig. 3 . Interestingly, as observed from the MD simulation and MM-PBSA calculation, the per-residue analysis further showed that during the last 40 ns, residue TYR118 significantly contributed to the interaction with ligand DB01698 having the contribution energy of −3.7 kJ/mol, thus, indicating unstable interaction or weak inhibition of viral replication. On the other hand, a total of 8 residues: THR25 (−8.2 kJ/mol), THR26 (−4.1 kJ/mol), LEU27 (−3.1 kJ/mol), HIS41 (−8.4 kJ/mol), CYS44–4.1 kJ/mol, THR 45 (−3.1 kJ/mol), CYS145 (−4.5 kJ/mol) and MET165 (−6.7 kJ/mol) highly contributed to the interaction with DB08995. It is important to note that DB08995 interacted with the catalytic dyad residues (HIST41 and CYS145) which suggests its strong inhibition. When DB12665 interacted with 3CLpro, a total of 5 residue: SER46 (−4.2 kJ/mol), MET165 (−8.1 kJ/mol), PRO168 (−3.0 kJ/mol), ASP187 (−3.3 kJ/mol) and GLN189 (−5.9 kJ/mol), were found to have highly contributed to the interaction energy. However, interaction between DB08995 and MET165 contributed to lower energy than other residues. Further analysis of per-residue interactions implies that the two drugs could be a potent inhibitor of SARS-CoV-2 targeting viral replication by strongly binding to the catalytic dyad and residues required for substrate binding.

Fig. 3

Per-residue energy decomposition analysis of the three ligands bound to 3CLpro analysed during the last 40 ns of MD simulation.

Targeting inhibition of Nsp9 RNA and nucleocapsid RNA

This group of proteins enables the RNA molecules to adapt to their conformational function [57]. Virtual screening was done to identify molecules that can inhibit their functions. Screening inhibitors targeting Nsp9 RNA showed that ligand 7 with binding free energy of −8.3 kcal/mol could potentially inhibit its activity. However, after screening targeting nucleocapsid RNA, all the ligands did not show potential inhibition, except for ligand 7 which showed binding free energy of −7 kcal/mol. It is suggested that ligand 7 could bind selectively to the non-structural proteins RNA as it did not show favourable binding free energy to the main protease (Fig. S5).

Targeting inhibition of PLpro and RdRp

All the three proteins; 3CLpro, PLpro and RdRp are involved in viral RNA synthesis and replication. PLpro is concerned with cleavages of the N-terminus of the poly-protein replicase to release Nsp1-3 which is required for correcting viral replication [58]. This protein is also responsible for antagonizing hosts innate immune system, [58] making it an important therapeutic target for identification of coronavirus inhibitors. Screening of our in-house natural products and drug obtained from a similarity search showed lower binding energy of −7 kcal/mol (Fig. S8a). The binding free energy suggested that these compounds would not be good inhibitors of PLpro. This motivated us to further investigate other targets responsible for viral RNA synthesis such as RdRp.

RNA-dependent RNA polymerase (RdRp) is an important enzyme of coronavirus responsible for replication and transcription of the virus. The research shows that RdRp is an important target as it is involved in the synthesis of Nsp12-RdRp RNA primer which activates the binding of Nsp12. Thus, targeting RdRp is crucial in the treatment of coronavirus. Virtual screening of natural products from our library and DrugBank showed that drug DB08995 could have potential inhibition with binding free energy −9.1 kcal/mol slightly corresponding to 0.2 μM (Fig. S8a). Natural product (ligand 1) from our library and drug CID442428 also demonstrated inhibition with the binding free energy of −8 and −8.5 kcal/mol, respectively. Binding mode analysis of the crystal structure revealed that DB08995 exhibited multiple binding modes. The binding mode with the lowest binding free energy as shown in Figure S8b did not form any hydrogen bonds with the target protein.

Targeting inhibition of host specific receptor ACE2

Previous studies have shown that enzymatic inhibition of ACE2 has a potential role in treatment of cardiovascular diseases such as hypertension, heart attack, and diabetes, as well as preventing the fusion of corona virus [59,60]. Therefore, we performed screening to target inhibition of ACE2 activity.

Virtual screening targeting inhibition of ACE2 revealed that both natural products from in-house library; ligand 1 and ligand 7, and those from FDA approved drugs; CID442428, DB01698 and DB8995 resulted to protein inhibition with the binding free energy of ≤ −9 kcal/mol (Fig. S8a). When all five compounds were subjected to RCS to evaluate the effects of the protein conformation on the binding energy, it turned out that only DB01698 and DB08995 showed strong inhibition with the binding free energy of −9.3 and −9.4 kcal/mol, respectively (Fig. 4 a). In order to establish the binding stability of DB08995 complex, MD simulation was carried out. RMSD measurement suggested that DB08995 was stable and did not induce conformation to the receptor, although it showed some changes inside the pocked as shown by ligand RMSD (Fig. 4b). Changes in ligand RMSD could be attributed to the large pocket size of the ACE2 receptor. The distance between ligand atoms and receptor residues that formed hydrogen bonds was measured to establish the origin of these changes. As shown in Fig. 4c, the sugar moiety 2 (red colour) exhibited more fluctuation within the pocket and it interacted with residues ASP350, GLU401 and GLU402 between the two helices as shown in Figure S9b. However, the middle ring remained stable with fewer fluctuations within the hydrogen bond distance cut-off. This was further confirmed by the number of hydrogen bonds formed during the simulation time. The number of hydrogen bonds formed is consistent with changes in distances (Fig. S9a).

Fig. 4

(a) Binding free energies of the six molecules against ACE2. (b) RMSD of ligand and receptor for ACE2 complex. (c) Hydrogen bond distances between specific groups that formed hydrogen bonds in the complex. (d) RCS binding free energy of DB01698 and DB08995 bound to spike RBD-ACE2. (e) RMSD for native and bound spike RBD-ACE2 (e) Hydrogen bonds formed for DB08995-spike RBD-ACE2 complex.

To investigate the stability, binding affinity based on MM-PBSA was calculated and results are presented in Table 3 . Binding free energy suggested that DB08995 strongly inhibited ACE2 activity. Both van der Waals and electrostatics energy terms favourably contributed to the interaction while polar energy terms unfavourably contributed to the interaction. These results suggest that, BD08995 could be potent inhibitor of ACE2 and thus it can prevent coronavirus from entering the host cells.

Table 3

Binding free energy (kJ/mol).

Energy	ACE2-DB08995
ΔEvdW	−185.2 ± 24.5
ΔEele	−92.4 ± 42.9
ΔEpolar	134.8 ± 29.5
ΔEnon-polar	−23.0 ± 1.9
ΔEbinding	−165.9 ± 31.9

Targeting inhibition of spike RBD-ACE2

BD08995 preferentially binds to ACE2 residues at the interface

Many studies have showed that inhibition of the ezymatic activity of ACE2 effectively blocks the interaction of spike SARS-CoV-2 and thus prevents the coronavirus fusion in the human cells [11,[61], [62], [63]]. It has been discovered that the virus spike RBD binds to ACE2 before entering the host [11]. Research also shows that SARS-CoV-2 binds strongly to ACE2 than SARS-CoV [61,64].

Screening results against the spike RBD-ACE2 showed that two drugs DB01698 and DB08995 could be potent inhibitors. RCS binding free energy showed that the two drugs DB01698 and DB08995 had a mean binding free energy of −9.33 and −9.59 kcal/mol, respectively (Fig. 4d). Interestingly, DB08995 exhibited two binding poses as compared to DB01698 (Fig. 5 a–b) and interacted better with ACE2 residues than DB01698.

Fig. 5

(a–b) Interaction and binding modes of DB08995 at the spike RBD-ACE2 interface. (c–f) Measured distances between DB08995 with spike RBD-ACE2 with their respective free energies to the right.

Binding of spike RBD to ACE2 does not abrogate the activities of DB08995

MD simulation was carried out to assess the binding stability of DB08995 and its ability to disrupt the fusion of the coronavirus into the host cells. DB08995 in its binding site B was selected as it showed lower binding free energy. The stability of the complex was assessed by measuring the RMSD, distances between the spike RBD and ACE2, drug and ACE2 hydrogen bonds as well as the binding energy at every 2 ns for the complex.

The measured RMSD for native spike RBD-ACE2 equilibrated to stability within 5 ns with an average RMSD of 0.24 nm. The RMSD further increased up to 0.29 from 25 ns and remained stable until the end of the simulation time. Although the RMSD of the bound spike RBD-ACE2-DB08995 equilibrated quickly, it showed persistent fluctuations until the end of the simulation (Fig. 4e). Comparing RMSD of native and bound spike RDB-ACE2, one can easily observe that DB08995 induced constant fluctuation to the native RBD-ACE2 throughout the simulation time, which suggests that binding of DB08995 could destabilize the spike RBD-ACE2 interaction, and therefore, prevent coronavirus membrane from stable binding and fusion into host cells. Hydrogen bonds formed between RBD-ACE2 and DB08995 are presented in Fig. 4f, which shows that DB08995 exhibited different binding poses. During the first 10 ns the maximum number of hydrogen bonds was 4 and increased up to 9 during 20 ns. The formation of hydrogen bond is important in stabilization of the interaction between the drug and its receptor. The number of hydrogen bonds later decreased to 5 during the end of the simulation time. In order to explain these changes, the distance between the drug and some residue was measured to monitor the stability and changes in binding pose (Fig. 5c–f).

To provide quantitative information, the free energy of the distances was calculated using F = −k Tln(P(dist)), where F is the free energy, k is Boltzmann constant, T is temperature and P(dist) is the probability. The free energy showed that DB08995 was bound to spike RBD forming hydrogen bond with residue GLU484, and it exhibited three minima at 0.28, 0.49 and 0.69 nm. Two deep minima at 0.28 and 0.49 nm suggests that, DB08995 bound to GLU484 disrupts its interaction with ACE2. Such inhibition is vital in preventing/altering molecular recognition between spike RBD-ACE. On the other hand, the measured distance between ACE2-DB08995 shows a deep minimum at 0.4 nm which suggests that the drug strongly bound and interacted with ACE2, hence preventing its recognition with the spike RBD.

DB08995 disrupt molecular recognition at the RBD-ACE2 interface

To obtain crucial insight into the disruptive nature of the DB08995 in relation to the spike RBD-ACE2 interface, the distance between RBD and ACE2 was measured for both native and bound protein (Fig. 6 a–b). Binding of DB08995 induced changes in distance between the two proteins. An increase in distance between the measured residues suggests weak interaction induced by the DB08995 at the interface which implies that spike RBD could weakly interact with ACE2 in the presence of DB08995. Time-lapse binding affinity of DB08995 into its cavity and spike RBD-ACE2 dissociation constant were also calculated so as to gain further understanding about the nature of the drug (Fig. 6c–d). The binding affinity of DB08995 in its cavity remained stable. However, the dissociation constant between spike RBD-ACE2 increased upon binding of DB08995 which suggests unstable interaction between spike RBD and ACE2.

Fig. 6

Interaction between DB08995 with spike RBD-ACE2. (a–b) Distance between RBD and ACE2 for native shown in (a) and bound shown in (b). (c) The binding affinity of DB08995 at the interface at every 3 ns (d) the dissociation constant between RBD and ACE2 for native and bound. (e) RMSF difference between ACE2 and RDB for bound and unbound structures.

The stability of the proteins was further investigated by analysing the atomic residue fluctuations. In addition, the dynamics of the residues were analysed so as to understand how the binding of DB08995 altered the protein fluctuations. In Fig. 6e we present the RMSF of the difference between the native and bound spike RBD-ACE2, which shows that the binding of DB08995 at the interface inhibited residue fluctuations at the region of 404–505 of RBD and induced fluctuations in the region of 370–390. The presence of the spike RBD to the ACE2 further altered atomic fluctuations to the protein. Previously, this region was reported to be involved in spike glycoprotein oligomerization for other coronavirus like murine hepatitis which forms dimeric domain required for receptor binding and entering the cells [65]. These results suggest that DB08995 may change the spike RBD dynamics which influence RBD oligomerization, receptor binding as well as recognition and hence treating the associated disease.

Water is known to play an important role in protein-ligand interaction, protein function and molecular recognition [66]. The role of water in mediating the interaction of drug and protein was also investigated. Figure S9 (left panel) shows that binding of DB08895 into ACE2 and spike RBD-ACE2 (Figure S9 left panel) was further mediated by water. Although water is observed at the interface between spike RBD-ACE2, it is not present in crystal. Such water formed a network of hydrogen bond between DB08995 and the spike RBD. The formed networks were important for DB08995 activity that may facilitate the loss of interaction and molecular recognition between spike RDB and ACE2.

In an effort to identify new drug to fight SARS-CoV-2 which causes COVID-19, a combination of computational methods were applied. Three FDA approved glycosylated flavonoids showed potent inhibition against SARS-CoV-2 targeting different proteins. DB12665 is an investigational drug for treatment of kidney cancer, renal cell carcinoma, advanced renal cell carcinoma, and thromboembolism of vein in pancreatic cancer [36]. This drug has also shown promising antiviral activities against influenza A & B viruses [67] and, therefore, it could be a potential candidate to be repurposed for SARS-CoV-2 targeting viral replication. DB01698 is an approved drug used to decrease capillary fragility [68]. This drug is also reported to have different pharmacological profiles including antiretroviral and antiviral activities [69]. It could also be repurposed against SARS-CoV-2 targeting inhibition of viral fusion into the host cells. DB08995 is an approved drug for venous disease and it also treats various disorders of blood vessels such as hemorrhoids, varicose veins, poor circulation in the legs, and bleeding in the eye or gums. In addition it serves as a food supplement [70]. The present work reveals that DB08995 target SARS-CoV-2 in two different ways; inhibiting viral replication and inhibiting viral fusion into host cells. This drug could further be repurposed for SARS-CoV-2. Recently, DB08995 was reported by other groups to have antiviral properties targeting inhibition on SARS-CoV-2 viral replication in silico [[71], [72], [73]].

Natural sources of the identified FDA approved drug potential for repurposing

The identified FDA approved drug molecules potential for repurposing are natural products belonging to the group of glycosylated flavonoids. This section provides the sources of these natural products for easy access and remedy against COVID-19. DB08995 is an aglycones flavone derivative [74] mostly found in green tea (Camellia sinensis; Theaceae) and herbal of Rosemary (Rosmarinus officinalis; Lamiaceae) with a concentration < 1–5 mg/serving. DB08995 is also found in juices and orange wines (Citrus sinensis; Rutaceae) in a concentration of 250 mL serving size. It is found in 0.2–0.5 mg/100 g of fresh weight and 1 mg/serving [74]. DB01698 is a citrus glycoside flavonoid found in many plants including green tea, passion flowers, apple and buckwheat [69]. It is an important component of food stuff [69]. DB12665 is highly found in mango (Mangifera indica L.; Anacardiaceae) [75], tea (Camellia sinensis; Theaceae) [76] and in the leaves of custard apple (Annona squamosa L.; Annonaceae) [77]. These computationally identified SARS-CoV-2 inhibitors may serve as cheap alternative remedy to combat COVID-19 once their efficacy is proved since they are plenty and easily available from natural sources.

Conclusion

Due to the urgent need to identify drug necessary in fighting COVID-19, we have carried out computational work to quickly identify FDA approved small molecules that could be considered for immediate use. To balance the computational accuracy, a combination of docking screening, molecular dynamics simulation and free energy calculations based on MM-PBSA methods were applied which identified DB12666, DB01698 and DB08995 as potential inhibitors of SARS-CoV-2. The later drug was discovered to work with two distinct mechanisms; targeting the disruption of viral replication/synthesis and inhibition of viral fusion into the host cells. The identified drugs are suggested for quick clinical evaluation to provide an immediate solution against this deadly disease.

Funding

This work was jointly supported by Data-Driven Innovation (DDI) project through CREATES scheme at NM-AIST, Arusha; the University of Dar es Salaam through its Directorate of Research Grant, and the Kenya Education Network (KENET) Mini-Grant 2018/2019.

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.