Protein (Receptor) structures. The crystal structures of NSP3 (PDB ID- 6W02), NSP9(PDB ID- 6W4B), NSP15 (PDB ID- 6VWW) from SARS CoV-2 were obtained from Protein Data Bank133. The crystal structures of COVID-19 main protease (NSP5) in complex with Z44592329 (ID:5r83) were obtained from the Protein Structure Database of Europe134. QHD434158: the I-TASSER models were obtained for NSP3(QHD43415_3), NSP5 (QHD43415_5), NSP11(QHD43415_11), NSP12 (QHD43415_12), NSP13(QHD43415_13), and NSP14 (QHD43415_14) 8.
Target Proteins (Receptors) preparation. Protein preparation was performed using the Protein Preparation Wizard in Maestro9. All co-crystallized atoms including calcium and chlorine were deleted. Co-crystalized ligands within active site were not deleted because they were used for grid generation. All water with less than 3 hydrogen bonds (Sample water orientation) with ‘non waters’ i.e. receptor protein or ligand were deleted. Each of the seven models developed was subjected to preparation by the Wizard in Maestro9. Each of the complexes was optimized with minimized hydrogens of altered species and then minimized using the OPLS3e force field9,135. Further, the models were subjected to MD simulation (MDS) after solvation in water and 0.15M NaCl (physiological solution). A simulation box covering the entire enzyme system was introduced with a 10 Å buffer space. The simulation was run for 100 ns at 300 K and standard pressure (1.01325 bar). The target complex trajectory analysis for interacting residues and was conducted using the Desmond software9.
Target library preparation. A virtual library of n=5903 drug compounds was downloaded from the latest Zinc database10 under the ‘world’ subset (approved drugs in major jurisdictions, including the FDA, i.e DrugBank approved). The database was then checked for redundancy, and duplicates were removed. Ligand preparation was performed using LigPrep, which generated variations of the ligands, eliminated reactive species, and optimized the ligands. Optimization was performed under the OPLS3e force field. EPIK minimization was performed on possible states at pH 7.0 ±2.0. Tautomers were generated for each ligand retaining specific chirality combinations with a maximum of 32 structures per ligand9.
Receptor grid generation. COACH analysis was performed to determine the location and size of the active site136. The shape and properties of the receptor are represented on a grid to ensure that possible active compounds are not missed. The centroid of the COACH predicted binding pocket residues and was used to generate grids using default values of protein atom scaling (1.0 Å) within a cubic box and ligand docking set to length of 20Å. The force field employed for grid generation was OPLS3e135. Rotatable groups were automatically determined by Receptor Grid Generation wizard of Glide module9. COACH predicted binding sites for each target were as follows; Exoribonuclease (ExoN); Based on the active site similarity with PDB ID : 5C8T, S-Adenosyl-L-Methionine (SAM) binding site was selected 292, 333, 335, 352, 353, 366, 367, 368, 385, 386, 387 & 426 PLpro (Papain-like proteinase); Based on PDB similarity with SARS-CoV Plpro model PDB ID: 3mj5B and catalytic triad of Cys–His–Asp the amino acid grid site selected was 209, 223, 230, 233, 234, 235, 236, 237, 238, 239, 293, 318, 321, 324, 325, 327, 349, 350, 352, 353, 354, 356 & 763. 3-Chymotrypsin-Like Protease (3CLpro); Based on PDB id: 2q6gB the substrate binding site was selected comprising of amino acid numbers 25, 26, 27, 41, 49, 140, 141, 142, 143, 144, 145, 163, 164, 165, 166, 168, 172, 187, 189, 190 & 192. 2’-O-Methyltransferase (2'-O-MT); Based on similarity with PDB id. 4N48 the SAM binding site was predicted to be 46, 74, 133, 134, 170 & 203 and was used for grid selection RNA-directed RNA polymerase (RdRP); Based on PDB hit; 3ol9I with similar active site the nucleic acid binding site was selected comprising of amino acid side chain residues 500, 501, 507, 512, 543, 545, 557, 559, 560, 569, 580, 589, 590, 591, 592, 682, 683, 684, 685, 686, 758, 857, 860, 864 & 914 Nonstructural Uridylate-specific endoribonuclease (NendoU): The COACH hit was PDB id:2C1W and the active site selected was 234, 249, 289, 291, 292, 293, 332, 334, 339, 340, 341 & 342 & Helicase (Hel): based on the similarity of active site with hit PDB id: 4N0O nucleic acid binding site was selected comprising of amino acid side chain residues 177, 178, 179, 180, 181, 197, 214, 309, 310, 311, 334, 335, 336, 337, 338, 339, 408, 409, 410, 411, 412, 413, 414, 485, 486, 516, 534, 554 & 560
High-Throughput Virtual Screening (HTVS). The screening was performed with default parameters including ionization states, Epik state penalties within the Glide (Grid-based ligand docking from energetics) module of Schrödinger suite137. The scaling factor was maintained at a default of 0.8 and a partial charge cut-off was limited to 0.15. The OPLS3e force field was used during the docking process. The HTVS ligand docking was the first to be performed, followed by SP and XP docking on the top 10% of scoring hits from each previous step. The XP docking aids in removing false positives and the scoring function is much stricter than the HTVS. The greater the XP Glide score, the better the calculated affinity of the hit in binding to the protein target. Further, the estimation of free binding energies for the best hit-docked complexes using MM force fields and implicit solvation was performed using the molecular mechanics/generalized Born surface area (MM-GBSA) method within the virtual screening workflow of Schrödinger suite9. The binding energy was calculated based on the following equation.
ΔG= E_complex(minimized) -(E_ligand(minimized)+ E_receptor(minimized))
The leads were ranked based on the binding free energy calculation for their respective protein-ligand complexes.
Molecular Dynamics Simulation (MDS). To determine the stability of the ligand receptor complex especially during flexiblity of the binding site an to determine the energetic strain on the docked ligand that could result in immediate flyoff (<10ns). The top hit for each target was subjected to a 20 ns MD simulation (100 ns for XAV-939 + 3CLpro, BIM IX + 3CLpro, & BIM IX + ExoN) as described above in the Target/receptor Preparation section to validate the interaction and thereby the HTVS pipeline and MM-GBSA performance. If the top hit flew off imediately due to major conformation change in the receptor site than the receptor was simulated alone for 20ns ro reach a stable state and screnning was repeated with the resulting minimized structure.
In-vitromicroneutralization assay.BIM IX was shortlisted for its predicted affinity for both 3CLpro (Main protease) and ExoN (Exonuclease) active site. Ivermectin was included as a positive control. Other shortlisted molecules included RdRP hits (paromomycin, amikacin, and lactulose), 3CLpro hits (iopromide & troxerutin), and one NendoU hit (haloperidol).
In vitro studies on SARS-CoV-2 were performed at the United States Army Medical Research Institute for Infectious Diseases (USAMRIID). Both entry and spread assays were performed on ATCC Vero E6 cells infected with the SARS-CoV-2 virus. Entry Assay: Vero E6 cells, in 96 well plates, were pre-treated with 14 different test compounds starting at a concentration of 50 uM with 3 fold dilutions down to 0.76 nM for approximately 1hr at 37°C. SARS-CoV-2/MT020880 was then added to the test compound treated cells for 1 hr. at 37°C at an MOI of 0.4. After 1 hr., the cells were washed with PBS before adding additional compounds back in fresh culture media to the cells for 24 hr. at 37°C. Cell media was removed and washed with PBS and plates were submerged in the formal fixing solution for 24 hrs., then permeabilized with 0.2% Triton-X for 10 min at RT and treated with blocking solution (3% BSA/PBS). Infected cells were detected using a primary detection antibody recognizing the SARS-CoV-2 nucleocapsid protein (Sino Biological, 40143-R001) and Alexa Fluor 488 conjugated antibody (goat α rabbit), then counterstained with DAPI for visualization of the cell nuclei. Infected cells were enumerated using Operetta high content imaging instrument and data analysis was performed using the Harmony software (Perkin Elmer). Spread Assay: A similar protocol was utilized to the entry assay described above with the following modifications. The virus was used at an MOI of 0.02 and the assay was incubated for 48 hrs. rather than 24 hr. after infection of the cells.
Compound screening to identify 3CLpro inhibitors.
Assays were performed in duplicate at room temperature in 96-well black plates at 25 °C. Reactions containing varying concentrations of inhibitor (1-50 µM) and 3CLpro enzyme (0.4 µM) in Tris-HCl pH 7.3, 1 mM EDTA were incubated for approximately five minutes. Reactions were then initiated with TVLQ-AMC probe substrate (40 µM), shaken linearly for 5 s, and then measured for fluorescence emission intensity (excitation λ: 364 nm; emission λ: 440 nm) over time (1 min-3 h) on a Synergy Neo2 Hybrid. Each assay contained 6-8 positive control wells (DMSO replacing inhibitor) and 2 negative control wells (assay components without 3CLpro).
3CLpro enzyme inhibition assay.
Fluorescence-based biochemical assay for inhibitors of 3CLpro. Inhibition assays were performed in 96-well black plate in triplicate at 25 °C. Reactions containing varying concentrations of inhibitor (0-1000 µM) and 3CLpro enzyme (0.4 µM) in Tris-HCl pH 7.3, 1 mM EDTA were incubated for approximately five minutes. Reactions were then initiated with TVLQ-AMC probe substrate (40 µM), shaken linearly for 5 s, and then measured continuously for fluorescence emission intensity (excitation λ: 364 nm; emission λ: 440 nm). Data were fit using nonlinear regression (dose-response inhibition, variable slope) analysis in GraphPad Prism 7.0.