Deep learning and docking

Interestingly, from S1 Table, we notice that DFCNN, DeepBindBC and Autodock Vina predict different drug compounds as top hits. This indicates that these three methods may be potentially complementary to each other through judging protein-ligand interactions from different perspectives. For instance, Amenamevir and Azithromycin have better Autodock Vina scores whereas Odanacatib and Nitisinone are found to have high DFCNN scores.

We first rejected those drugs that had poor prediction by any of the methods. The 22 drugs with DFCNN score above 0.9 and docking score bellow -7 kcal/mol were firstly selected and presented in S1 Table. Among the 22 drugs, we further excluded the drugs with a DeepBindBC score below 0.7, resulting in 14 drugs shown in bold in S1 Table. Overall there are 5 nucleoside analogues, 2 antibiotic drugs, 3 antivirus drugs, 2 anticancer drugs and 3 other drugs (S1 Fig) selected by the molecular vector-based and structure-based screening process. Sofosbuvir is both a nucleoside analogue and an antivirus drug. These 14 drugs are subject to force field-based screening in the next stage.

Among 2 antibiotic drugs, Azithromycin, a drug used to treat a variety of bacterial infections, showed top Autodock Vina score of -8.6 kcal/mol, good DFCNN score and DeepBindBC score (0.9093 and 0.8589), respectively. Gautret, P. et al claimed that combined with Hydroxychloroquine, Azithromycin can have good efficiency in treating COVID-19 with significant viral load reduction [17]. However it should be noticed that currently there is no evidence of the effectiveness of Azithromycin in the treatment of COVID-19 and have many debates about effective of Azithromycin on COVID-19 [18].

The top two predicted molecules by DeepBindBC are nucleotide analogues. Sofosbuvir is a nucleotide analogue inhibitor of hepatitis C virus (HCV) NS5B polymerase to treat infectious liver disease, [19] whereas Clofarabine is a purine nucleoside antimetabolite used for treating refractory acute lymphoblastic leukaemia [20]. More nucleotides analogues in the candidate list were selected by our method, such as Adenosine, Vidarabine, and Gemcitabine, indicating some RdRp-nucleotides interaction patterns have been implicitly recognized by the proposed hybrid drug screening method.

Force field-based simulation

To further screen the 14 selected drugs (2D structures shown in S1 Fig) and understand their interactions and stability, we have performed MD simulations on RdRp-drug complexes. The structure stability is estimated by Root Mean Square Deviation (RMSD) over 100ns time scales, shown in S2 Fig. The drugs such as Azithromycin, Pralatrexate, Romidepsin, Teriflunomide and Vidarabine are found stable indicated by the minimum RMSD fluctuations. On the other hand, the drugs such as Adenosine, Amenamevir, Fipronil, Gemcitabine and Sofosbuvir have high RMSD fluctuations. The number of hydrogen bonds formed between RdRp and the drug in the 100 ns MD simulation was also computed and shown in S3 Fig. Among the 14 RdRp-drug complexes, Pralatrexate clearly showed more hydrogen bonds with the RdRp than others. We found that there is no direct correlation between RMSD fluctuation and size (or number of initial contacts). For instance, the Romidepsin has very large size (molecular weight: 544.7 g/mol), but shows less binding stability compared to Pralatrexate (molecular weight: 479.5 g/mol). Also, Vidarabine have molecular weight of 271.3 g/mol, shows more stable binding compared to Sofosbuvir (molecular weight: 527 g/mol) and Amenamevir (molecular weight: 484.6 g/mol). In MD simulation, small size ligands can also have very tight binding depending on the ligand and binding site (pocket) properties [21]. Whether the binding is stable is more related to the complementary between pocket and ligand in both geometric and physio-chemical features. So the RMSD calculations can help to provide useful information about the binding stability.

The binding free energy vs coordination number (CV: collective variable) from metadynamics simulations is shown in S4A Fig. The lowest energy conformations of protein-drug complexes for Amoxicillin, Azithromycin, Pralatrexate and Sofosbuvir showed more contacts in the interface region, as indicated by the high coordination numbers (S4A (iii) Fig), while most other compounds favor smaller coordination number (close to zero) indicating no or weak interactions (S4A (i) Fig).

Detailed interaction patterns between RdRp and the four most optimal compounds (Azithromycin, Pralatrexate, Amoxicillin and Sofosbuvir) are shown in Figs 1B (i)(ii) and S4B (i)(ii), whose structures are taken from the last frame of the 100 ns MD simulation. Azithromycin and Pralatrexate interact with 16 amino acid residues of RdRp to form a stable complex. The RdRp-Azithromycin complex is mainly dominated by van der Waals interactions, whereas Pralatrexate involves more polar and charge interactions. According to the calculated free energy difference (ΔG) values from the metadynamics simulations between the unbound state and the binding state for Amoxicillin, Azithromycin, Pralatrexate and Sofosbuvir (S2 Table), Azithromycin and Pralatrexate (-305.8 kJ/mol, -128.6 kJ/mol) show more favorable binding energy than Amoxicillin and Sofosbuvir (-67.3 kJ/mol and -89.9 kJ/mol). Remdesivir has covalent bond with the RNA primer, our method can only estimate binding free energy between non-colvalent-binding protein and ligand. Furthermore, our system does not contain RNA premier, so the interaction mechanism and pattern of the Pralatrexate and Azithromycin would be different from Remdesivir.

It is noted that all the nucleoside analogues highly recommended by deep learning-based screening methods were excluded from the force field based screening process. The possible explanation is that our protein-drug systems do not contain the RNA primers during the MD simulation, and without RNA the nucleoside analogues would not form base pair like interaction, hence showing no binding in MD simulation.

Interaction patterns (taken from the last frame of 100 ns MD simulation) between RdRp and Azithromycin, Pralatrexate are shown in Fig 1B (i)(ii), respectively. The interaction between Remdesivir in its monophosphate form and RdRp (PDB ID: 7BV2) is also given for comparison in Fig 1B (iii). Azithromycin forms 2 hydrogen bonds with GLN573 and ILE494 through keto and hydroxyl groups respectively, and many hydrophobic related interactions with the RdRp binding site (e.g. LYS577, LEU576, ALA685) through alkyl groups, whereas Pralatrexate shows enhanced and more stable interactions with RdRp binding site, including 6 hydrogen bonds with GLN573, ARG569, ASN496, ASN497, LYS500 and GLY590. Pralatrexate also forms Alkyl or Pi-Alkyl interaction with LYS577 and LEU576, and salt bridges with ARG569 and LYS500. Azithromycin and Pralatrexate share 10 common neighbor residues (62.5%) of RdRp, as shown in Fig 1B(i)(ii) and S3 Table, indicating a similar binding cavity.

In a recent study, the authors have shown that the monophosphate active form of Remdesivir interacts with the RdRp and covalently incorporated into partial double stranded RNA template of RdRp at the +1 position [6]. Shown in Fig 1B (iii) and S3 Table. The key interacting contacts between Remdesivir and RdRp include ARG555. Remdesivir only shares 2 common close contact residues with Pralatrexate, and no common close contact residue with Azithromycin, suggesting a different interaction pattern.

This screening procedure relies on non-covalent bond interactions and assumes RdRp is in apo form (without RNA primer), Fig 1B(i)(ii) shows Azithromycin and Pralatrexate interacting with the region of RdRp which consists of residues such as GLN573, ARG569, ALA685. Comparing our modeled structure with Remdesivir in its monophosphate-RdRp complex, it is possible Azithromycin and Pralatrexate occupy part of the RdRp cavity with the non-covalent binding which may interfere the entry of the RNA primer strand to the cavity.

Pralatrexate and Azithromycin inhibit the replication of SARS-CoV-2 in vitro

To further confirm the efficiency of the hits from the virtual screening, we tested the antiviral activity of Azithromycin, Pralatrexate, Amoxicillin and Sofosbuvir in vitro. Experiments were performed in a biosafety level 3 laboratory where regulation requires. Vero cells were infected with SARS-CoV-2 (BetaCoV/Shenzhen/SZTH-003/2020, GISAID No. EPI_ISL_406594) at a MOI(multiplicity of infection, which represents the ratio of the numbers of virus particles to the numbers of the host cells in a given infection medium.) of 0.02 (the cytopathic effect was mild at 48 hours post-infection with this MOI) in the presence of varying concentrations of the tested drugs, and the inhibition rates were evaluated by quantification of viral copy numbers in the cell supernatant via quantitative reverse transcription polymerase chain reaction reverse transcription polymerase chain reaction (qRT-PCR) and confirmed with immunofluorescence assay (Fig 2). The results showed that Pralatrexate and Azithromycin could efficiently inhibit the replication of SARS-CoV-2, with half-maximal effective concentration (EC₅₀) values of 0.008 and 9.453 μM (Fig 2A), whereas Remdesivir achieved an inhibitory activity with EC₅₀ value of 8.777 μM within the same experimental system (S5 Fig). Indirect immunofluorescence assay (IFA) showed similar results with qRT-PCR assay (Fig 2B). CCK-8 assay of the two drugs showed that the half-cytotoxic concentration (CC₅₀) values of Pralatrexate and Azithromycin on Vero cells were 0.167 μM and > 100 μM, respectively, and the calculated the selectivity indexes (SI) of Pralatrexate and Azithromycin were 20.878 and >10.579, respectively. Whether the two drugs worked at the stage of viral entry or post entry was analyzed using time-of-addition assay as previously reported [5]. The results showed that Pralatrexate functioned at a stage post virus entry, while Azithromycin functioned at both entry and post-entry stages of the SARS-COV-2 infection in Vero cells (Fig 2C). Furthermore, surface plasmon resonance (SPR) experiments were performed to test the in vitro binding of Pralatrexate and Azithromycin with immobilized RdRp protein of SARS-CoV-2. Both drugs showed expected binding response in S6 Fig.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The antiviral activities of the test drugs against SARS-COV-2 in vitro.

A, Vero cells were infected with SARS-COV-2 at an MOI of 0.02 in the presence of the indicated concentrations of the tested drugs for 48 hours. The viral yield in the cell supernatant was then quantified by qRT-PCR. Meanwhile, cytotoxicity of these drugs to Vero cells was measured by CCK-8 assay. B, Immunofluorescence microscopy of virus infection upon treatment of Pralatrexate and Azithromycin at the indicated concentrations. IFA was performed at 48 hours post-infection. Scale bar, 100 μm. Cells were immunostained for the Viral protein (green) and DNA (blue). C, Time-of-addition experiment of Pralatrexate and Azithromycin. Cells were infected with SARS-COV-2 at an MOI (multiplicity of infection, represents the ratio of the numbers of virus particles to the numbers of the host cells in a given infection medium.) of 0.02 with different treatment, and virus yield in the infected cell supernatants was quantified by qRT-PCR. For the group of “Entry”, the drugs were added to the cells for 1 hour before viral attachment, and at 2 hours post-infection, the virus-drug mixture was replaced with fresh culture medium. For the group of “Post-entry”, drugs were added at 2 hours post viral infection, and maintained until the end of the experiment. For the group of “Full-time”, Vero cells were pre-treated with the drugs for 1 hour, and the virus was then added to allow attachment for 2 hour. Afterwards, the virus-drug mixture was removed, and the cells were cultured with drug-containing medium until the end of the experiment.

https://doi.org/10.1371/journal.pcbi.1008489.g002

Structural modeling of RdRp and drug compound dataset

The RdRp sequence and its modelled structure were obtained from https://zhanglab.ccmb.med.umich.edu/C-I-TASSER/2019-nCov/. The RdRp-ligand model was constructed by I-TASSER [26]. The ligand was taken from the template protein (PDB ID: 3BR9) [27] by COFACTOR algorithm [28] within I-TASSER using structure comparison and protein-protein networks. We extract the amino acids within 1 nm of the ligand as the binding pocket. RMSD between the modeled structure and the recent experimental RdRp structure (PDB ID 6M71) is calculated (~0.516 Å), shown in S11A Fig [10]. RNA polymerase superfamily region is also very similar between these two structures (RMSD = 0.456 Å, S11A Fig).

TargetMol-Approved_Drug_Library, which contains 1906 compounds, was used as virtual screening library. These 1906 compounds collected by TargetMol are drugs approved by Food and Drug Administration (FDA), the European Medicine Agency (EMA), or China Food and Drug Administration (CFDA), or included in the US Pharmacopeia (USP) Dictionary, the British Pharmacopoeia (BP), the European Pharmacopoeia (EP), the Japanese Pharmacopoeia (JP), or Chinese Pharmacopoeia (CP) Dictionary.

Machine-learning based drug screening

Pocket molecular dynamics and metadynamics

Further drug screening was carried out by force field based molecular dynamic (MD) simulations. The initial protein-drug complexes was from top score conformation Autodock Vina docking, the ligand was edited by pymol software [32] to make it in correct protonation state at pH 7. In this study, we selected 14 drug binding complexes for MD simulation, including Adenosine, Amenamevir, Amoxicillin, Azithromycin, Clofarabine, Fipronil, Gemcitabine, Nitisinone, Pralatrexate, Raltegravir, Romidepsin, Sofosbuvir, Teriflunomide and Vidarabine, respectively.

We also refined a pocket molecular dynamics simulation (pocket MD, S11B Fig) to facilitate the simulation process by only keeping the binding pocket region for simulation. Binding free energy calculation can be estimated by metadynamics simulations to explore whether protein-ligand will bind in solution. Metadynamics relies on addition of a bias potential to sample the free energy landscape along a specific collective variable of interest [33,34]. Note that the binding free energy calculations from Metadynamics may only be suitable for detect the general trend of binding in virtual screening.

The pocket MD is same as the classical MD simulation, except that we only using the pocket region to reduce system size for simulation (S11B Fig), which is inspired by a previous dynamic undocking (DUck) method [35]. An in-house script was used to extract the pocket region of the protein (1nm within the binding ligand), the N terminal and C terminal ends were capped with the ACE and NHE terminals, respectively. We applied a position restrain to the ACE and NHE terminals to maintain the relative conformation of the pocket. MD simulation was carried out by Gromacs with AMBER-99SB force field [36,37]. The topology of ligand and the partial charges of ligand was generated by ACPYPE [38], which relies on Antechamber [39]. Firstly, we created a dodecahedron box and put the target-ligand complex at the center. A minimum distance from the protein to box edge was set to 1 nm. We filled the dodecahedron box with TIP3P water molecules [40], the counter ions were added to neutralize the total charge using the Gromacs program tool [41]. The long-range electrostatic interactions under the periodic boundary conditions was calculated with Particle Mesh Ewald approach [42]. A cutoff of 14 Å was used for van der Waals non-bonded interactions. Covalent bonds involving hydrogen atoms were constrained by applying the LINCS algorithm [43].

We performed the energy minimization steps with a step-size of 0.001ns, 100 ps simulation with isothermal-isovolumetric ensemble (NVT), and 10ns simulation with isothermal-isobaric ensemble (NPT) for water equilibrium. After that, a 100ns NPT production run (step size 2 fs) was carried out. The Parrinello-Rahman barostat and the modified Berendsen thermostat were used for simulation with a fixed temperature of 308 K and a pressure of 1 atm. RMSD and hydrogen bond number of the trajectory were calculated using Gromacs tools.

The simulation was continued using the metadynamics approach for exploring the free energy landscape. The interface coordination number of atoms of protein ligand complex was used as collective variable (CV). The protein-ligand interface coordination numbers correlate with the numbers of atom contact, and larger coordination number usually indicates that protein-ligand is in binding state.

The coordination number C is defined as follows by Plumed: (1) and (2)

In the simulation, n was 6, m was 12, d₀ was 0 nm and r₀ was 0.5 nm. d₀ is a parameter of the switching function. r_ij is the distance between atom i and atom j. The degrees of contacts between two groups of atoms can be estimated by above function(1) [44]. Metadynamics simulation for each protein-ligand system was performed for 100 ns (except protein-Azithromycin, which was extended to 300ns in order to reach the 0 Coordination Number and achieve convergences). During the metadynamics simulation, Gaussian values were deposited every 1 ps with a height of 0.3 kJ/mol. The widths of the Gaussians were 5 for the coordination number. The free energy landscapes of the metadynamics simulations along the CV were generated by the Plumed program and plotted using Gnuplot [45].

Tools used in analysis

The USCF Chimera, VMD, ICM-browser-Pro (http://www.molsoft.com/icm_browser_pro.html) and Discovery Studio Visualizer 2019 were used to generate the structure and to visualize the 2D protein-ligand interactions [46–48]. Clusfps (https://github.com/kaiwang0112006/clusfps) which depends on RDKit [49] was used to cluster the drugs in the dataset. The drug fingerprint was used as inputs with algorithm of Murtagh [50] being used for clustering 1906 drugs into 20 groups.

Cell culture

Vero cell (ATCC, CCL-81) was cultured at 37°C in Dulbecco’s modified Eagle’s medium (DMEM, Gibco) supplemented with 10% fetal bovine serum (FBS, Gibco) in the atmosphere with 5% CO₂. Cells were seeded in 96-well plates and cultured overnight with a density of 5 × 10⁴ cells/well prior infection or drug feeding. Remdesivir, Azithromycin, Pralatrexat, Sofosbuvir and Amoxicillin were obtained from Selleck Chemicals. All drugs were dissolved in DMSO to prepare 50 mM stock solutions, and stored at -20°C. DMSO was used in the controls.

Viral stock titration by 50% tissue culture infective dose (TCID₅₀)

TCID₅₀ was measured as previously reported [51]. In brief, Vero cells in 96-well plates were grown to 80% confluence and infected with 10-fold serial dilutions of the stock SARS-CoV-2 (BetaCoV/Shenzhen/SZTH-003/2020, GISAID No. EPI_ISL_406594) for 1 h at 37°C. The inoculum was removed, and cells were overlaid with fresh DMEM plus 2% FBS. At 5 days post infection (d.p.i), plates were assessed for the lowest dilution in which 50% of the wells exhibited cytopathic effects. The values of TCID₅₀ were calculated according to the Reed-Muench method [52].

Evaluation of antiviral activities of the drugs in Vero cells

Firstly, the cytotoxicity of the five drugs on Vero Cells were determined by CCK8 assays (Sangon). Then the antiviral activities of the drugs were evaluated as previously reported with some modification [5]. Vero cells seeded in 96-well plates were pre-treated with the different doses of the indicated drugs for 1 h, and then virus was subsequently added at multiplicity of infection (MOI) of 0.02 to allow infection for 2 h. Then, the virus-drug mixture was removed and cells were further cultured with fresh DMEM with 2% FBS and the indicated concentrations of drugs. At 48 hours post infection (h.p.i), the cell supernatant was collected and viral RNAs were extracted using the QIAamp RNA Viral Kit (Qiagen, Heiden, Germany) for further quantification analysis. The cells were collected for indirect immunofluorescence assay (IFA). All the experiments involving infectious SARS-CoV-2 were handled in BSL-3 facilities at the Shenzhen Third People's Hospital.

Quantitative reverse transcription polymerase chain reaction

This assay was carried out as described previously [53]. Viral RNAs were extracted from the samples using the QIAamp RNA Viral Kit (Qiagen, Heiden, Germany), and quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using a commercial kit (Genrui-bio) targeting the S and N genes. The specimens were considered positive if the Ct value was ≤ 38.0, and negative if the results were undetermined. Specimens with a Ct higher than 38 were repeated. The specimen was considered positive if the repeat results were the same as the initial result and between 38 and 40. If the repeat Ct was undetectable, the specimen was considered negative.

Indirect immunofluorescence assay (IFA)

IFA was carried out as previously reported [54,55]. Vero cells were fixed in 4% formaldehyde at 48 hours post infection. Then cells were permeabilized in 0.5% Triton X-100, blocked in 5% BSA in PBS, and then probed with the plasma of this patient or healthy control at a dilution of 1:500 for 1 h at room temperature. The cells were washed three times with PBS and then incubated with either goat anti-human IgG conjugated with Alexa fluor 488 at a dilution of 1:500 for 1 h (Invitrogen). The cells were then washed and stained with hoechest-33342 (Invitrogen) to detect nuclei. Fluorescence images were obtained and analyzed using EVOS FL Auto Imaging System (Invitrogen).

Protein expression and purification

The genes for nsp12 of SARS-CoV-2 isolate BetaCov/Wuhan/WH01/2019 (EPI_ISL_406798) was chemically synthesized with codon optimization for insect cells (Spodoptera frugiperda) by Synbio Technologies. The sequence was fused with a C-terminal thrombin cleavage site, a 6×His-tag and a 2×Strep-tag, and incorporated into pFastbac-1 plasmid. Recombinant protein was expressed with Hi5 cells at 27°C. Cells were harvested at 48 hpi(hour post infection)and resuspended in 25 mM HEPES pH 7.4, 1 M NaCl, 1 mM MgCl₂ and 2mM TCEP. An equal volume of the same buffer supplemented with 0.2% (v/v) Igepal CA-630 (Anatrace) was added and incubated at 4°C for 10 min. Cells were lysed by sonication and the lysate was clarified by ultracentrifugation. Cleared lysates were passed through a 0.22-μm filter film before further purification. The protein was purified by tandem affinity chromatography and SEC.

Surface plasmon resonance (SPR) assay

The affinities between nsp12 and drugs were measured at room temperature (r.t.) using a Biacore 8K system with CM5 chips (GE Healthcare). The nsp12 protein was immobilized on the chip with a concentration of 100 μg/mL (diluted by 0.1mM NaAc, PH 4.0).

Drug samples were prepared according to procedure 29264621AA of GE Healthcare Life Sciences. 1×PBS solution plus 5% DMSO and 0.005% p20 was used for running and diluting drugs. A blank channel of the chip was used as the negative control. Serial diluted drugs were then flowed through the chip surface. The LMW multi-cycle kinetics was analyzed with the Biacore 8K Evaluation Software (version1.1.1.7442) and fitted with a 1:1 binding model.

Statistical analysis

Data are presented as the mean ± SD (Standard Deviation). All analyses were performed using GraphPad Prism version 7.0 for Windows (GraphPad Software, San Diego California, USA). Data were subjected to statistical analysis by two-way ANOVA or two-tailed Student’s t-test. The P values less than 0.05 were considered statistically significant.