key: cord-222664-4qyrtzhu
authors: Coban, Mathew; Morrison, Juliet; Freeman, William D.; Radisky, Evette; Roch, Karine G. Le; Caulfield, Thomas R.
title: Attacking COVID-19 Progression using Multi-Drug Therapy for Synergetic Target Engagement
date: 2020-07-06
journal: nan
DOI: nan
sha: 
doc_id: 222664
cord_uid: 4qyrtzhu

COVID-19 is a devastating respiratory and inflammatory illness caused by a new coronavirus that is rapidly spreading throughout the human population. Over the past 6 months, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the virus responsible for COVID-19, has already infected over 11.6 million (25% located in United States) and killed more than 540K people around the world. As we face one of the most challenging times in our recent history, there is an urgent need to identify drug candidates that can attack SARS-CoV-2 on multiple fronts. We have therefore initiated a computational dynamics drug pipeline using molecular modeling, structure simulation, docking and machine learning models to predict the inhibitory activity of several million compounds against two essential SARS-CoV-2 viral proteins and their host protein interactors; S/Ace2, Tmprss2, Cathepsins L and K, and Mpro to prevent binding, membrane fusion and replication of the virus, respectively. All together we generated an ensemble of structural conformations that increase high quality docking outcomes to screen over>6 million compounds including all FDA-approved drugs, drugs under clinical trial (>3000) and an additional>30 million selected chemotypes from fragment libraries. Our results yielded an initial set of 350 high value compounds from both new and FDA-approved compounds that can now be tested experimentally in appropriate biological model systems. We anticipate that our results will initiate screening campaigns and accelerate the discovery of COVID-19 treatments.

COVID-19 is a disease cause by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It was identified in Wuhan city, in the Hubei province of China in December 2019 (Chen et al., 2020; Huang et al., 2020; Zhu et al., 2020) . The virus is spread between people via small droplets produce by talking, sneezing and coughing. The disease was declared a global pandemic by the World health organization (WHO) on March 11th, 2020. While a large proportion of the cases results in mild symptoms such as fever, cough, fatigues, loss of smell and taste, as well as shortness of breath, some cases progress into more acute respiratory symptoms such as pneumonia, multiple-organ failure, septic shock and blood clots. These more severe symptoms can lead to death and are likely to be precipitated by a cytokine storm after infection and multiplication of the virus in humans. Indeed, recent data indicate that the levels of IL-6 correlate with respiratory and organ failures (Gubernatorova et al., 2020) . So far, the estimated death rate of SARS-CoV-2 is above 1.3%, which is more than 10 times higher than the death rate of seasonal influenza (Abdollahi et al., 2020) . Older patients and patients who have serious underlying medical conditions such as hypertension, diabetes, and asthma are at higher risk for severe disease outcomes . A clear understanding of the genetics and molecular mechanisms controlling severe illness remains to be determined. SARS-CoV-2 is a positive-sense, single-stranded RNA betacoronavirus, closely related to SARS-CoV-1, which caused severe acute respiratory syndrome (SARS) in 2003, and Middle East respiratory syndrome coronavirus (MERS-CoV), which caused MERS in 2012. Positive-strand RNA viruses are a large fraction of known viruses including common pathogens such as rhinoviruses that cause common colds, as well as dengue virus, hepatitis C virus (HCV), West Nile virus. The first genome sequence of SARS-CoV-2 was released in early January on the open access virological website (http://virological.org/) (Zhou et al., 2020) . Its genome is ~29.8 kb and possesses 14 open reading frames (ORFs), encoding 27 proteins (Wu et al., 2020a) . The genome contains four structural proteins: spike (S) glycoprotein, envelope (E) protein, membrane (M) protein, and nucleocapsid (N) protein. The E and M proteins form the viral envelope, while the N protein binds to the virus's RNA genome. The spike glycoprotein is a key surface protein that interacts with cell surface receptor, angiotensinconverting enzyme 2 (Ace2) mediating entrance of the virus into host cells (Zhu et al., 2018) . In addition to its dependence on the binding of S to Ace2, cell entry also requires priming of S by the host serine protease, transmembrane serine protease 2 (Tmprss2). Tmprss2 proteolytically processes S, promoting membrane fusion, cell invasion and viral uptake (Heurich et al., 2014; Hoffmann et al., 2020) . Blocking viral entry by targeting S/Ace2 interaction or Tmprss2-mediated priming may constitute an effective treatment strategy for COVID-19. The non-structural proteins, which include the main viral protease (nsp5 or M pro ) and RNA polymerase (nsp12), regulate virus replication and assembly. They are expressed as two long polypeptides, pp1a and pp1ab, which are proteolytically processed by M pro .

The key role of M pro in viral replication makes it a good therapeutic target as well. A third group of proteins are described as accessory proteins. This group is the least understood, but its members are thought to counteract host innate immunity (Kim et al., 2020, Cell 181, 914-921) (Fig. 1A) .

There is currently no treatment or vaccine available to prevent or treat COVID-19 (Baden and Rubin, 2020; Lurie et al., 2020) (https://www.fda.gov/news-events/press-announcements/coronavirus-covid-19update-daily-roundup-june- . While the FDA has granted emergency use authorization (EUA) for the 65-year-old antimalarial drug, hydroxychloroquine, COVID-19 treatment based on early results from clinical trial in China and France Gautret et al., 2020a; Gautret et al., 2020b; Million et al., 2020) , more recent results reported that hydroxychloroquine does not decrease viral replication, pneumonia or hospital mortality, and may in fact increase cardiac arrest in patients infected with COVID-19 Rosenberg et al., 2020) . The accuracy of the statistical analyses in these studies raised serious concerns in the scientific community. More accurate data are needed to reach a conclusion about the effect of hydroxychloroquine in COVID-19 patients. In another recent study published in the New England Journal of Medicine, the antiviral remdesivir, an unapproved drug that was originally developed to fight Ebola, seemed to improve patients with severe breathing problems (Beigel et al., 2020) and has also recently been granted EUA by the FDA. Repurposing drugs that are designed to treat other diseases is one of the quickest ways to find therapeutics to control the current pandemic. Such drugs have already been tested for toxicity issues and can be granted EUA by the FDA to help doctors to treat COVID-19 patients.

Another efficient way to attack the virus is to use drug cocktails to target multiple enzymes/pathways used by the virus. Combination therapy has the advantage of being less likely to select for treatmentresistant viral mutants. Such a strategy has been successfully used to treat hepatitis C virus (HCV) and human-immunodeficiency virus (HIV) infections. In the case of HCV, the treatment, Enpclusa, combines sofosbuvir, which inhibits the viral RNA-dependent RNA polymerase (NS5B), and velpatasvir, a defective substrate that inhibits NS5A. Antiretroviral therapy (ART) against HIV combines drugs from different drug classes to target disparate aspects of the HIV replication cycle. These drug classes include nucleoside reverse transcriptase inhibitors, non-nucleoside reverse transcriptase inhibitors, protease inhibitors, fusion inhibitors, CCR5 antagonists, post-attachment inhibitors, and integrase inhibitors. One example from HIV-AIDS literature is the randomized comparison of 4 groups of patients comparing monotherapy to combination therapies: zidovudine (ZDV) monotherapy; ZDV zidovudine and didanosine; ZDV plus zalcitabine; or didanosine monotherapy. This randomized trial showed positive results when ZDT was combined with didanosine or zalcitabine, and for didanosine compared to ZDT monotherapy in raising CD4 counts greater than 50% (Hammer et al., 1996) . Combination therapy has become standard of care initial treatment in other infectious diseases such as Mycobacterium tuberculosis and failure to cure with monotherapy and requires multidrug therapy (MDT) (Collaborative Group for the Meta-Analysis of Individual Patient Data in et al., 2018) . Similar MDT is also found effective in hepatitis C virus infection using glecaprevir and pibrentasivr combination therapies which lead to sustained virological response rates as far out as 12 weeks' post-treatment (Wang et al., 2019) .

We propose an effective combination therapy for COVID-19 could target the SARS-CoV-2 replication cycle at multiple levels to synergistically inhibit viral spread and dissemination. Using a computational pipeline that aimed to expeditiously identify lead compounds against COVID-19, we combined compound library preparation, molecular modeling, and structure simulations to generate an ensemble of conformations and increase high quality docking outcomes against two essential SARS-CoV-2 viral proteins and their host protein interactions; S/Ace2, Tmprss2, Cathepsin L and K, and M pro that are known to control both viral binding, entry and virus replication (Fig. 1A) . Our in silico approach (Fig.  1B) , which will most likely lead into experimental virus screening, structural characterization of binding interactions by X-ray crystallography, and compound safety profiling. Virtual screening (VS) is a rational driven controller for identification of new hits from compound libraries (Willett, 2006) using either ligandbased (LBvs) or structure-based (SBvs) virtual screening (Dror et al., 2004) . LBvs tactics use structural and biological data of known active compounds to select favorable candidates with biological activity from experiments (Jahn et al., 2009; Maldonado et al., 2006) . SBvs approaches, on the other hand, examine quantitative structure-activity relationships (QSAR), clustering, pharmacophore and 3D shape matching (Villoutreix et al., 2007) . The utility of VS is evident in the growth of our knowledge base of new compounds and existing drugs as well as the expansion of our structural databases. SBvs is generally the preferred approach when access to the target 3D-information derived from NMR, X-ray crystallography or homology models (Jahn et al., 2009; Maldonado et al., 2006) is possible. Molecular docking (docking) is the most common SBvs approach used today (Bottegoni et al., 2009; Corbeil et al., 2012; Fernandez-Recio et al., 2005; Friesner et al., 2006; McGann, 2012; Morris et al., 2009 ) and searches for the ideal position and orientation (called "pose") of the small molecule within a target's binding site, which gives a score for the pose. When including knowledge of experimentally known compounds ("actives") from a 3D target, LBvs and SBvs can be combined to increase likelihood of obtaining new actives from searches (Kruger and Evers, 2010) .

Hit identification in VS also requires careful selection of the methods used based on the goal of the project (e.g. compound databases and libraries can be either proprietary, commercial or public) (Bender, 2010) . ZINC is one such large public database often used in VS (Irwin and Shoichet, 2005) , which contains millions of compounds. By contrast, other libraries have structure-activity relationships (SAR) databases (Scior et al., 2007) that integrate information about compound interactions with their known targets. DrugBank, Chem-Space are other attractive sources of compounds for drug repurposing (or repositioning) (Ashburn and Thor, 2004; Duenas-Gonzalez et al., 2008; O'Connor and Roth, 2005) (Wishart et al., 2008) , and maintain drug diversity that is useful for scaffold development (Gozalbes et al., 2008; Schreiber, 2000) .

Advances in computing power have increased utility of in silico screening capabilities and balanced the need for accuracy with virtual high-throughput screening approximations and assumptions (Anthony, 2009; Lee et al., 2008; McGaughey et al., 2007; Plewczynski et al., 2009) , while recent techniques have improved accuracy without sacrificing CPU time (Caulfield and Devkota, 2012; Caulfield et al., 2011; Jiang et al., 2014; MacKerell et al., 1998; Phillips et al., 2005) (Fig. 1B) . Further innovations in docking methods have improved the exactness of empirical docking equations (Corbeil et al., 2012; Fernandez-Recio et al., 2005; Friesner et al., 2006; Kalid et al., 2012; Kruger and Evers, 2010; McGann, 2012) . Accuracy is improved by incorporating molecular flexibility with simulations (Caulfield, 2012; Caulfield et al., 2019; Caulfield and Medina-Franco, 2011; Caulfield et al., 2011; Caulfield et al., 2014; Kayode et al., 2016) , thus capturing conformational information on structural changes that directly impact compound docking results.

To target the COVID19 problem on multiple fronts (e.g. Ace2:S protein, Tmprss2, M pro , and Cathepsin L and K), as well as improve our screening accuracy using our selected repurposing libraries and new chemical entity libraries (ZINC database), we implemented a novel method that integrates protein flexibility/shape, adaptive biasing algorithms, machine learning from drug data, and final Z-score matrix weighting to our drug modeling. We matched all FDA compounds with our realistic (X-ray derived) protein structures over a dynamic range of protein conformations with accelerated dynamics using our algorithms, such as Maxwell's demon molecular dynamics (MdMD); this approach combines docking with simulations for exploration of both ligand and protein flexibility (Caulfield, 2012; Caulfield et al., 2019; Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield et al., 2014; Kayode et al., 2016; von Roemeling et al., 2018) . We then refined the drug-target interface our specific leaderlike hit compounds using the quantum mechanics (QM)-based scoring within our MdMD matrix (Caulfield, 2012) to make our go/no-go assessment, which is particularly useful with NCEs and de novo compound design (DCDs). The protocol for library, structural modeling, dynamics, refinement, and hit identification as part of a pipeline is given (Fig. 1B) .

To improve our docking outcome, we constructed x-ray structure-based models of Ace2 bound to Sprotein, M pro , and Tmprss2 in our molecular dynamics simulations (MDS) and virtual screening (Fig.  1B,S1 ). As S-protein interfaces with Ace2 at a distinct region from the active site ( Fig. S1A-D) , inhibition of the binding site by ligands may disrupt the Ace2/S-protein interaction. Canonical inhibitors of Ace2 bind at the active site where angiotensin interacts, whereas drugs directed at the structural region for S-protein binding are not overlapping with the binding site. The modulation of Ace2/S-protein interaction by canonical Ace2 inhibitors is likely allosteric and suboptimal. Therefore, directly targeting the interface of the interaction should increase efficacy of the approach and block COVID viral binding, precluding entry (Fig. S1) . Additional investigation into the glycosylation sites of the S-protein demonstrated that the Ace2 binding site is mostly unaffected by these additions (Fig. S2) .

A. S-protein:Ace2 interaction (protein-protein inhibitor, PPI) requires dynamics to reveal binding site To get the optimal interface for drug screening, we used our grid searching algorithms, as well as site mapping and protein-protein docking, to examine the protein-protein interactions surface using MDS ( Fig. 2-3 ,S1) (Bhachoo and Beuming, 2017; Caulfield and Devkota, 2012; Caulfield et al., 2011; Caulfield and Harvey, 2007; Fernandez-Recio et al., 2005; Kozakov et al., 2006) . The protein-protein inhibitor (PPI) interaction complex did not identify any immediate binding site on the surface of the PPI interfaces. Nevertheless, a small pore around one single beta-sheet in the center of the PPI interaction area could be exploited as a weak point that may perturb the interface equilibrium. Using UniProt, which contains information about a number of confirmed mutations, we determined the relative potencies of PPI binding residues, identifying those that would likely affect the integrity of the complex (Fig. 2) . Residues K353 and Y41, which interact with D155 at the center of the PPI, are likely stabilizing its surface, potentially forming a useful "hot spot" for targeted druggability ( Fig. 2-3,S2) .

To check whether this is true and to understand how Ace2:S-protein cooperation functions, we performed two MD simulations, one with and one without the mutation of Y41A. This mutation causes strict inability to form the S-protein:Ace2 complex. Analysis of the trajectory of the wild-type protein, which possessed an intact complex, revealed the three most stable conformations of the "hot spot" region with expanded pores inside the triangle of residues K353, D155, Y41. Since it is impossible to determine which of these three conformations is the most stable, we ran three high-throughput screenings based on the donor-acceptor atoms and hydrophobic areas of the region. We then performed three MD simulations with top pose ligands. As demonstrated in Figure 3K , ligands failed binding within 10 ns, while docked ligands became leaders, as determined by energetic stability, during MD and interaction energy values (electrostatic -red, Van der Waals -blue) ( Fig. 3J/L) .

To identify inhibitors of the S-protein:Ace2 interaction via docking, we used the best scoring compounds obtained after combination of molecular docking and molecular dynamics simulations, which feeds into the pipeline for constraint-based screening. The high-throughput screening (HTS) of a PPI library did not produce any results, since the PPI binding sites were weakly identified shallow regions (Fig. 2,4A -D,S1). Compounds that made good insertion into the sites situated between Ace2 and S-protein were able to perturb the association of S-protein with Ace2 via steric hindrance of S-protein association (Fig.  3) . From the MDS, we detected compounds that decreased energy of stability between the Ace2:Sprotein complex, which is desired in an inhibitor of protein-protein interaction. As a whole, this approach identified a deep and narrow binding site to disturb the S-protein interaction with Ace2 (Fig. 3,4A-D) .

C. Tmprss2 and M pro modeling requires dynamics to reveal optimal inhibitor binding To optimize the binding site of our inhibitors, we constructed a full-length (zymogen) model of Tmprss2 (epitheliasinogen), as well as a mature version of the protease (epitheliasin), as described in our method section ( Fig. 4E-G) . The mature protease model was used for MDS studies to generate a reference dynamical profile that can be used to assist in silico screening of Tmprss2 inhibitors. A control experiment was also completed with the uncleaved (non-catalytic) form of Tmprss2 to demonstrate the pocket's instability and poor ligand binding capacity ( Fig. S3 ) (Ko et al., 2015; Lucas et al., 2014; Wilson et al., 2005) . A full-length model of monomeric M pro was also constructed, as well as a homodimer ( Fig. 4H-K,S1 ). The structure derived from PDB code 6Y2F with its ligand was used for a consensus virtual screen . In addition, we used the dimer to generate a reference dynamical profile to assist with in silico screening and study its interdomain behavior.

We acquired the dimer protein sequence from the UniProt database. BLAST search showed the highest identification values against factor XI, prothrombin, kallikrein proteases (~41-42%). However, we focused on ligands that could be active against active form of Tmprss2 protein. Thus, we found the ligand:

(2s)-1-[(2r)-2-(Benzylsulfonylamino)-5-Guanidino-Pentanoyl]-N-[(4-Carbamimidoylphenyl)methyl]pyrrolidine-2-Carboxamide, contained within the ChemblDB repository (CHEMBL1229259) and active against Tmprss2, prothrombin, and Factor XI. Likewise, another docked model was recovered with macrocyclic ligand (CHEMBL3699198), called: Ethyl14-[[(E)-3- [5-chloro-2-(tetrazol-1-yl) phenyl]prop-2-enoyl]amino] -5-(methoxycarbonylamino)-17-oxo-8,16 diazatricyclo[13.3.1.02,7] nonadeca-1(18),2(7), 3,5,15(19) -pentaene-9-carboxylate. We launched several molecular dynamics simulations (up to 75 ns of duration) to understand the interaction with the target protein-binding site. Figure S3 shows the initial and stable/final states of our various models ( Fig. 4E-G) . The MD analysis provided useful results for selecting the appropriate model. After 15 ns MD, the putative binding site collapsed (Fig. S3,4E-G) . Although the active form of thrombin was used for Tmprss2 modeling, as a negative control we also examined the region with prothrombin-based binding site for completeness of the docking study (Fig. S3) . The overlay of the average homology model structure from MD and structure 3F68 (PDB code) was used as a template to compare protein-ligand interaction map and assign docking constraints (Baum et al., 2009) . Two optimal inhibitors for Tmprss2 were selected for demonstration purpose in Figure 5 . We also modeled Cathepsins L and K for preliminary work, since these can be implicated in late-endosomal entry of the virus (Fig. S4 ).

For the viral main proteinase, M pro , a key enzyme for coronavirus replication (SARS-CoV-2), and a potential target for anti-SARS drug development, several peptidomimetics synthetized in early 2012 against SARS-CoV-1 proteases were identified as selective. There is a high degree of sequence identity between the SARS-CoV-1 and SARS-CoV-2 M pro . This means that SARS-focused ligands could form similar interaction map with M pro protein and offers good launching points for 3D-QSAR/Machine Learning-drive based drug design for future iterations. To perform the virtual screening, protein structure was taken from the PDB code 7BQY complex and significant attention was paid to the interaction between the crystallized ligand from the complex and protein-binding site (Jin et al., 2020) (Fig. 6) . As the binding site is quite large (Fig. 6A) we used a set of additional crystal structures (PDB code 6Y2F and fragment-like compounds from https://www.diamond.ac.uk/) to narrow the source of possible conformations. The binding of the compounds inserted into this region demonstrated a very canonic and recurring interacting motif, represented with α-Keto amide group flanked with aliphatic or saturated rings. We then performed molecular dynamics of 75 ns for the ligand-free dimer structure of the M pro to evaluate and "catch" the most flexible elements of the binding site. Our simulation revealed that the extended binding pocket was not very stable, unlike its individual sub-pocket, which contains active cysteine (C145) residue (Fig. 6B,6C) . We began our molecular docking after assigning several combinations of constraints that should define specific interactions with the protein-binding site. We performed several high-throughput screening procedures using the same set of features in different combinations of constraints by partial matching algorithm ( Fig. 6D-E) . We then ranged docking scores and compared obtained conformations inside the binding site with the co-crystalized ligands from 7BQY, 6Y2F structures to select the most potent compounds.

By disrupting the SARS-CoV-2 viral process in three different critical routes: Binding, Entry, and Replication with our virtual screening approaches against dynamic structures, we were able to identify 350 compounds (Dataset S1) and compile data reflecting physiochemical and chemoinformatic properties. An exemplar top hit from each target is summarized for docking score in Table 1 . To classify the compounds and their chemical space, we completed various regression, K-means analyses and fingerprint measurements, and provide further details about their structures and properties, including commonly evaluated traits: MW, HBA, HBD, docking score, Rule of Three (Jorgensen), Rule of Five (Lipinksi), logP o/w , and logS (Dataset S2). We focused on new compound searches. The MW for these initial screening compounds ranges from large fragment (~250 Da) to mature drug sized molecules (~500 Da) with only 10 of the 310 top scoring compounds being over 500 Da in size and the smallest fragment-based compound measured 178 Da. Overall the docking scores were very good with median around -7 kcal/mol using the Glide XP calculations. We also generated a list of most commonly related drugs and discuss some of our best hits to known and clinical trial drugs (Dataset S3). The general process for pruning the >30 million total chemical fragments and compounds from commercially available compounds for the initial round of virtual screening is described (Fig. 1B) , which reduces the primary large set to 3 million per conformation of target. As an example, when examining some prototype compounds from our selected dataset of >300 NCEs screened from >10 million total compounds, we find the predicted interactions between drug and protein ( Table S2 ) have some common binding modalities. When looking at the dynamical data for the drugs binding to the protein-protein site on Ace2, we find the RMSD, RMSF, and H-bond occupancy evidence strong binding capability, as calculated from three separate simulations of Ace2 with different ligands, referred to as 300, 392, and 488 ( Fig. 4,S1) . These observations can be applied to generate constraints for additional virtual screening to improve the performance at higher throughput. Based on these results, ligand 392 reduced the overall RMSD and per residue RMSF, while maintaining strong hydrogen bonds, as demonstrated by its greater occupancy during the simulation (Table S1 ). This information, particularly H-bond occupancy and modulation of interface residue RMSFs, can be used in conjunction with docking and other data to profile the compounds more thoroughly (Fig. 4) . In some cases, where constraints were utilized, the docking score underrepresents the compound and testing is needed to get important single-point data to clarify actives from non-actives, as well as determine the real IC50s for the selected active compounds. We will enrich our dataset with the top compounds for future rounds of parallel chemical screening and eventual de novo chemical design for novel chemical entities. Current results of our approach are presented on all three targets (Ace2, Tmprss2, M Pro ).

For each of our targets, we screened for hits from a library of FDA-approved compounds alongside the more extensive library of NCEs. Our final result across all three targets identified a total of 350 specific compounds, with 167 against Ace2, 40 against Tmprss2, and 103 against M pro . Among these are FDAapproved drugs that could be repurposed: 21 against Ace2, 11 against Tmprss2, and 8 against M pro (Supplemental Dataset TableS1_topNCE-FDA-hits.xlsx).

Isoprenaline hydrochloride (isoprotenerol) is an adrenoreceptor agonist that can be repurposed as a vasopressor to augment cardiovascular function with a beta-receptor side benefit of bronchodilation to improve breathing function. Metaraminol bitartrate, a stereoisomer of meta-hydroxynorephedrine, is a potent sympathomimetic amine to raise blood pressure. Atenolol and nadolol are beta-receptor blocking agents used in chronic hypertension, a comorbid risk factor in COVID-19 patients. Propafenone is an anti-arrhythmic agent approved for patients with life-threatening ventricular tachycardia. Levosulpiride is an atypical antipsychotic medication with prokinetic function that can be used in patients with agitated delirium, and gut immotility. Valganciclovir hydrochloride is an antiviral agent used for cytomegalovirus (CMV), varicella zoster virus (VZV), and preventative medication in HIV patients (Wu et al., 2020b) . Recent data shows COVID-19 deplete CD8 T helper cells similar to HIV (Zheng et al., 2020) . Amikacin sulfate and cephalexin are antibiotic anti-bacterial drugs that can treat bacterial super-infection. Prochlorperazine dimaleate is a phenothiazine derivative prescribed in medicine for nausea. Isoetharine mesylate is a selective adrenergic beta-2 agonist and fast-acting aerosolized bronchodilator for COVID-19 respiratory distress. Benserazide hydrochloride is an aromatic L-amino acid decarboxylase (DOPA decarboxylase inhibitor) used with levodopa for the treatment of Parkinsonism. Glucosamine hydrochloride is constituent found in cartilage and used for osteoarthritis joint pains. S4701 or 2-Deoxy-D-glucose (2D-DG) compound can induce ketogenic state, a powerful pathway involved in reducing systemic inflammation. Inulin is a natural prebiotic agent that enhances GI function and digestion by increasing prebiotic GI homeostasis critical to stabilize downstream anti-inflammatory effects and prevent overgrowth of harmful bacteria. Metaproterenol is a bronchodilator (beta-2 receptor agonist) that is commonly used to treat a variety of respiratory disorders including asthma, COPD, bronchitis and wheezing associated with viral pneumonias in clinical practice. The novelty of this drug is that is aerosolized and can be given as a breathing treatment and similar reach the lungs, which have a tremendous surface area and enter the blood rapidly. By inhalation this drug acts rapidly and potentially with or in combination with other aerosolized drugs or oral or IV combination drugs. Its inhalational route of delivery also can reach alveolar type II cells which express Ace2 for dual synergism. Metaraminol bitartrate, a stereoisomer of meta-hydroxynorephedrine, is a potent sympathomimetic amine. This drug is used in patients with hypotension or low blood pressure. COVID-19 hospitalized patients in the intensive care unit (ICU) setting often need vasopressor agents to raise blood pressure in a condition called shock (dangerously low blood pressure) from COVID-19 disease or sepsis. Therefore, metaraminol has dual purpose of antiviral function at Ace2 docking site /entry as well as helping with systemic blood pressure in those acutely ill COVID-19 patients. This drug has immediate repurposing use in this patient population.

Atorvastatin is a statin drug with anti-inflammatory, immunomodulatory (Diamantis et al., 2017) and endothelial benefits (Ackermann et al., Varga et al., 2020) . Carbenicillin disodium is a penicillin derivative antibacterial antimicrobial agent. Catechins are derived from plants with many beneficial properties in human health including anticancer, anti-obesity, antidiabetic, anti-cardiovascular, antiinfectious, hepatoprotective, and neuroprotective effects (Isemura, 2019) . These substances fall outside FDA purview since supplements and generally have a wide safety margin that will be tested on the multidrug platform. Epicatechine S5105 is a naturally occurring flavonoid found in chocolate with anti-sarcopenic effects on skeletal muscle (Gutierrez-Salmean et al., 2014) . Ivosidenib is an experimental drug for treatment of several forms of cancer. Bezafibrate is a fibrate lipid-lowering drug, which creates a favorable anti-inflammatory ratio against cardiovascular diseases. PF299804 or dacomitinib is an EGFR inhibitor used in cancer therapeutics. Metaproterenol is a bronchodilator (beta-2 receptor agonist) that is commonly used to treat a variety of respiratory disorders with viral pneumonias in clinical practice. Carbenicillin disodium is a penicillin derivative antibacterial antimicrobial agent that as mentioned above can be used in conjunction with other anti-SARS-Cov-2 agents to shut down antiviral effects and used in combination with those COVID-19 patients with secondary super-infection with bacterial infection of lung, blood, or skin.

Bumetanide is a loop-diuretic used to remove extra fluid in the body (edema) such as pulmonary edema. Aloin is an anthraquinone glycoside found naturally in aloe vera plants, a natural cathartic, and decreases 16s rRNA sequencing of dysbiosis-producing butyrate producing bacterial species via an emodin breakdown product (Gokulan et al., 2019) . Emodin blocks Ace2 and viral docking (Ho et al., 2007) . Salbutamol sulfate (albuterol) is a bronchodilator used in various breathing disorders. S4953 usnic acid is a naturally occurring dibenzofuran derivative found in lichen plant species, in some kombucha teas, with adrenergic function to raise blood pressure and potential bronchodilator. Usnic acid is an active ingredient in some and a preservative in others and has a wide array of antimicrobial action against human and plant pathogens with antiviral, antiprotozoal, antiproliferative, antiinflammatory, and analgesic activity (Ingolfsdottir, 2002) . Avanafil is a class of medications called phosphodiesterase (PDE) inhibitors, which are pulmonary artery and circulation dilators. S3612 Rosmarinic acid is a naturally occurring compound found in plants (rosemary and sage), which has broad range of antimicrobial activity including antiviral activity including HIV (Shekarchi et al., 2012) . Ractopamine is a beta-agonist function used for bronchodilatation. Neohesperidin dihydrochalcone (NHDC) is a naturally derived plant sweetener (bitter orange) with anti-Tmprss2 effects. Cidofovir and zidovudine (ZDV) are both antiviral drugs used in HIV patients.

There is a critical unmet patient need for therapeutics to treat the acute phase of COVID-19 disease now and for the future. Efforts to create and trial a vaccine are underway, but 11.6 million patients are confirmed infected globally (>540K deaths) with 25% infected within the United States and we are just at the midpoint of 2020. Therefore, there is an urgent need to rapidly speed drug discovery from the bench to the bedside. In order to accelerate drug discovery, translation and human application, a design funnel using high-powered artificial intelligence is needed to screen millions of compounds against macromolecular mechanistic targets against the virus. At the back end of this funnel 40 drug candidates emerged, many of which may represent repurposing candidates for use in humans due to known safety and tolerability profiles. However, the approach with the highest probability of overall clinical therapeutic success may be not a single drug therapy for this viral RNA disease but rather a multi-pronged drug approach gleaned from decades of HIV-AIDS epidemic research. A multidrug approach for HIV has improved survival, markedly reduced viral loads, and vastly improved management of the disease by preventing AIDS end-stage fatal complications. We therefore suggest that a multifaceted drug approach for SARS-Cov-2 may prove superior by attacking 3 viral entry and replication cycle sites simultaneously: Ace2 receptor docking site and entry, Tmprss2 endosomal packaging, and M Pro viral replication. Multiple drug targets for each of the 3 sites also allow permutations and optimization for combinatorial success.

A recent study that screened commercially available >10,000 clinical-staged and FDA-approved small molecules against SARS-CoV-2 in a cell-based assay (Riva et al., 2020) identified interesting compounds for alternative targets that complement our results. These FDA approved compounds included MDL-28170, a selective Cathepsin B inhibitor; VBY-825, a non-specific Cathepsin B, L, S, V inhibitor; Apilimod, an inhibitor of production of the interleukins IL-12 and IL-23; Z-LVG-CHN2, a tripeptide derivative inhibitor for cysteine proteinases; ONO 5334, a selective Cathepsin K inhibitor; and SL-11128, a polyamine analogs designed against E. cuniculi, a antimicrobial agents used as an adjuvant treatment for opportunistic AIDS-associated infections. Overall these compounds are Cathepsin-centric or antibiotic in nature, with little to no effect on our intended targets (Tmprss2, Ace2, M Pro ). Additional top hits identified by Riva et al. include: AMG-2674, an AMGEN compound inhibitor of TRPV-1 (Vanilloid Receptor); SB-616234-A that possesses high affinity for human 5-HT1B receptors; SDZ 62-434 that strongly inhibited various inflammatory responses induced by lipopolysaccharide (LPS) or function-activating antibody to CD29; Hafangchin A (also called "Tetrandrine"), a bisbenzylisoquinoline alkaloid, which acts as a calcium channel blocker; Elopiprazole an antipsychotic drug of the phenylpiperazine class (antagonist at dopamine D2 and D3 receptors and an agonist at serotonin1A receptors) that was never marketed; YH-1238, which inhibits dipeptidyl peptidase IV (DPP-IV) enzyme prolonging the action of the incretin hormones, glucagon-like peptide-1 (GLP-1) and glucose-dependent insulinotropic polypeptide (GIP); KW-8232, an anti-osteoporotic agent that can reduce the biosynthesis of PGE2; Astemizole, an antihistamine; N-tert-butyl Isoquine (also called "GSK369796"), an antimalarial drug candidate; and Remdesivir, a broad-spectrum antiviral medication developed by the biopharmaceutical company Gilead Sciences. Again, none of these compounds were geared toward targeting Tmprss2 or Mpro, and are also not specific to Ace2. While the lack of overlap may be surprising, results generated by Riva and colleagues are not in opposition to our findings and both approaches can complement each other. Most importantly, these approved FDA compounds can be combined with our set of identified NCE (310 compounds) that have been demonstrated to have low toxicity issues based on our chemoinformatics filtering (Fig. 1B) . All NCE compounds identified were chemical moieties that do not overlap any FDA drugs. Altogether, the data presented here complements previously generated data and should help prioritize and rapidly identify safe treatments for COVID-19. Future work will rely on advanced 3D-QSAR, fragment-based drug design principles for de novo drug optimization.

Among millions of potential COVID-19 drugs screened the majority of the final 40 drug candidates have known medical use and/or FDA approval for a primary indication (e.g., hypertension, cardiac indication, hyperlipidemia) with well-established patient safety and tolerability profiles from large phase III human trials and post-market (Phase IV) analyses. These large human data provide both a clinically significant and scientifically innovative window of opportunity to test 40 compounds on the multidrug platform, and, in conjunction, observe longitudinal human survival outcomes of COVID-19 patients on these drugs for comparative effectiveness within established and ongoing patient registries. An emerging example of this important parallel is Ace2 pathway drugs (Ace inhibitors [AceI] and angiotensin receptor blocking drugs [ARB]), which are increasingly observed in humans with COVID-19 to be associated with improved survival advantage (Jarcho et al., 2020; Mancia et al., 2020; Mehta et al., 2020; Patel and Verma, 2020; Vaduganathan et al., 2020) . However, there is a scientific knowledge gap within human registries data regarding a scientifically robust and testable translational platform to test mechanistic effects of these different molecular compounds. Therefore, creation of a "pandemic platform" using newer technology of compound AI drug throughput screening combined with animal multi-drug screening models creates an early Phase I/II safety, tolerability and early efficacy platform which is rapidly needed to expedite bedside human use for COVID-19 pandemic, and as a platform that can be used in future pandemics.

A flurry of activity to identify compounds for SARS CoV-2 targets has been underway by academic labs globally. Here in our approach we introduce our novel Maxwell's demon molecular dynamics method for screening flexibility required to get rare and essential conformational transitions and pathways to find the most likely druggable state. We also used our quantum docking technique (QM-driven adaptive molecular dynamics scanning docking) (Caulfield, 2012) to identify compounds effective for targeting Ace2, Tmprss2 and M pro . The compounds identified by our large-scale in silico platform can next be experimentally validated as binders for intended targets and for efficacy in models of the disease, evaluated for EC50/safety-toxicity data, and carried into hit-to-lead and lead optimization in a drug development pipeline. Structural studies such as X-ray crystallography will also be important to generate structural SAR data for these efforts.

In sum, our leading edge in silico methods incorporating structural dynamics have produced a set of 350 candidate compounds suitable for screening in biological disease models. Among these, 40 FDAapproved compounds are eligible for rapid clinical trial testing. Additionally, our results bring forward 310 NCEs predicted to possess potency and specificity for viral or human accessory target proteins to lower the viral load. Moreover, this resource offers the community a set of chemical tools to probe the behavior of these enzymes essential for SARS-CoV-2 progression, namely, binding, entry and replication. As SARS-CoV-2 is already endemic, the rapid identification of effective antivirals remains a paramount focus until we have an efficient vaccine to provide long-lasting protection.

In general, COOT was used for building in missing residues and regularizing geometry (Emsley and Cowtan, 2004; Emsley et al., 2010) . More details for the preparation of each model are given in the respective subsections. Since these structures were all used in downstream computational studies, a uniform structural preparation was implemented. The full-length structures are comprised of all residues and side chains. We added missing atoms in rotamers and de-clashed atoms, added missing residues for chain continuity, and removed extraneous molecules/atoms (e.g. artifacts of crystallography or alternative conformations of residues were removed (keeping the highest occupancy)), and the Bfactors were set to isotropic. The PDBePISA server was used to data mine the interface between Ace2 and S-protein (Krissinel and Henrick, 2007) . Surface interactions data is provided (Supplemental). Calculations on molecular dynamics trajectories including RMSD, RMSF, and H-bonds were performed using VMD and internal tools thereof (RMSD trajectory tool and Tk Console). Prior to calculations, the backbone (CONCα) atoms of each frame of the trajectories were aligned to the first frame as a reference, to remove the effect of random rotation/translation. After alignment, the per residue average of RMSF or RMSD per frame in Å across the entire MDS trajectory is given. For the Ace2-ligand simulations, the number of hydrogen bonds between the protein and ligand were recorded for each frame, and the occupancy of each specific H-bond is defined as the percentage of frames the bond is present. RMSD, RMSF, and H-bond data were plotted in 2D format in Excel. The RMSF was also appended to the beta column of the PDB and heat-mapped to the structure using a custom Tcl/Tk script and PyMOL. All molecular graphics were generated in PyMOL (Mooers, 2016) .

Molecular Dynamics and Monte Carlo simulations were performed on the protein to allow local regional changes for full-length structure for all acids of each structure. The X-ray refinement for Monte Carlo was built using YASARA SSP/PSSM Method (Altschul et al., 1997; Hooft et al., 1996a; Hooft et al., 1996b; King and Sternberg, 1996; Krieger et al., 2009; Qiu and Elber, 2006) . The structure was relaxed to the YASARA/Amber force field using knowledge-based potentials within YASARA. The side chains and rotamers were adjusted with knowledge-based potentials, simulated annealing with explicit solvent, and small equilibration simulations using YASARA's refinement protocol (Laskowski RA, 1993) . The entire full-length structure was modeled, filling in any gaps or unresolved portions from the X-ray.

Refinement of the finalized models was completed using either Schrodinger's LC-MOD Monte Carlo-based module or NAMD2 protocols. These refinements started with YASARA generated initial refinement of Tmprss2 (Altschul et al., 1997; Hooft et al., 1996a; Hooft et al., 1996b; Krieger et al., 2009) . The superposition and subsequent refinement of each protein regions yields a complete model. The final structures were subjected to energy optimization with PR conjugate gradient with an Rdependent dielectric.

Atom consistency was checked for all amino acids of the full-length wild-type structure, verifying correctness of chain name, dihedrals, angles, torsions, non-bonds, electrostatics, atom typing, and parameters. Model was exported to the following formats: Maestro (MAE), YASARA (PDB). Model manipulation was done with Maestro (Macromodel, version 9.8, Schrodinger, LLC, New York, NY, 2010), or Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) .

MDS and MC searching were completed on each model for conformational sampling, using methods previously described in the literature (Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield, 2011; Caulfield et al., 2011) . Briefly, each protein system was minimized with relaxed restraints using either Steepest Descent or Conjugate Gradient PR, then allowed to undergo the MC search criteria, as shown in the literature (Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield, 2011; Caulfield et al., 2011) . The primary purpose of MC, in this scenario, is examining any conformational variability that may occur with each protein.

For Ace2/S-protein, PDB code 6VW1 was used to construct the model (Shang et al., 2020) . While the structure was mostly complete, chain F (S-protein) was missing more residues, though it had residue Ala522. Chain E (S-protein) was only missing residue 522. Residue Ala522 was built into chain E using COOT and where the extraneous molecules (solvent/cryoprotectant) and chains were deleted to leave only the heterodimer Ace2/S-protein, which was processed to be used for computational studies, not to generate a de novo model or complete structure with missing atoms and sections. All information about the protein was found on the corresponding Uniprot page. After identifying the hot spot residues using SiteMap or protein-protein interfaces, we used MD to find out how the Y41A mutation can affect of PPI inhibition. We performed MD for wild type and mutated protein. Residual mutation was also performed using PyMol's built-in tools. Gromacs 2018 and amber99 force field were used to conduct MD and further analysis of the results (Baugh et al., 2011; Dilip et al., 2016; Janson et al., 2017; Kazmierkiewicz, 2013, 2016; Mooers, 2016) . Visual inspection of every 10 frames allowed us to determine some tendency of structural deformation in a certain place on the protein surface. According to the literature data and our finding, we focused on the predicted binding site. Then, each trajectory was analyzed via the built-in clustering tool based on the RMSD distribution. Three the most stable conformations of the binding site were chosen for the docking studies. All received docking poses from each docking study were evaluated based on the docking scores, interaction diagrams and solvent exposure. To make some prediction regarding the binding method, we carried out another molecular dynamics simulation for the upper poses of each docking. After such a confirmation of our assumptions, we selected the most powerful and accurate compounds from the results of docking.

A homology model was constructed on the basis of prothrombin crystal structure in complex with the ligand analog (PDB code 3F68) (Baum et al., 2009) . We modeled the 492 amino acid Tmprss2 protein two different ways: YASARA based and SwissModel server based (Krieger et al., 2002; Waterhouse et al., 2018; Zoete et al., 2011) . First, the YASARA based model begins with the FASTA sequence: MALNSGSPPAIGPYYENHGYQPENPYPAQPTVVPTVYEVHPAQYYPSPVPQYAPRVLTQASNPVVCT QPKSPSGTVCTSKTKKALCITLTLGTFLVGAALAAGLLWKFMGSKCSNSGIECDSSGTCINPSNWCDG VSHCPGGEDENRCVRLYGPNFILQVYSSQRKSWHPVCQDDWNENYGRAACRDMGYKNNFYSSQGI VDDSGSTSFMKLNTSAGNVDIYKKLYHSDACSSKAVVSLRCIACGVNLNSSRQSRIVGGESALPGAWP WQVSLHVQNVHVCGGSIITPEWIVTAAHCVEKPLNNPWHWTAFAGILRQSFMFYGAGYQVEKVISHPN YDSKTKNNDIALMKLQKPLTFNDLVKPVCLPNPGMMLQPEQLCWISGWGATEEKGKTSEVLNAAKVLL IETQRCNSRYVYDNLITPAMICAGFLQGNVDSCQGDSGGPLVTSKNNIWWLIGDTSWGSGCAKAYRP GVYGNVMVFTDWIYRQMRADG. Topological domains have the following characteristics: residues 1 -84 forms the cytoplasmic sequence; residues 85 -105 form the transmembrane domain region (helical 21 aa); and residues 106 -492 form the Signal-anchor for type II membrane protein (extracellular), where the protein as two main chains: non-catalytic chain (Met1-Arg225) and catalytic chain (Ile256-Gly492), where each domain modeled as a separate unit built together in composite. Disulfide bonds exist between several residues (113 ↔ 126), (120 ↔ 139), (133 ↔ 148), (172 ↔ 231), (185 ↔ 241), (244 ↔ 365), (281 ↔ 297), (410 ↔ 426), (437 ↔ 465), which can be informative for building the structure. Glycosylation sites are also possible at residues N213 and N249. Cleavage site (active) exists between Arg255 and Ile256 (see refinement section).

The second method, homological modeling was performed using the SwissModel server after performing a BLAST search on available protein structures in the RCSB database. Molecular dynamics simulations of 100 ns of both, suggested and re-modeled protein structures, was performed with GROMACS 2018 Kazmierkiewicz, 2013, 2016) . Based on the structural analysis and the generated Connolly surfaces, we identified critical changes in the binding site of the proposed model and began creating a mesh for the binding site of the new homology model. Since our model was based on the structure of thrombin, we used its co-crystallized ligand as a template for assigning constraints and ensured we built the catalytically active state.

For M pro (PDB 6Y2F) co-crystallization with tert-butyl (1-((S)-1-(((S)-4-(benzylamino)-3,4-dioxo-1-((S)-2oxopyrrolidin-3-yl)butan-2-yl)amino)-3-cyclopropyl-1-oxopropan-2-yl)-2-oxo-1,2-dihydropyridin-3yl)carbamate (also referred to as alpha-ketoamide 13b) was used; the structure was also mostly complete. Residues E47 and D48 were built in using COOT, where the other preparations previously described were also performed. To build the missing residues, the coordinates and structure factors were downloaded, generated 2mFo-DFc and FEM maps, and real space refine zone/regularize zone were used to fit to electron density and optimize local geometry. The ligand (alpha-ketoamide 13b) was left for usage as a cognate ligand for virtual screening.

The protein structure was initially studied using MD to find out if the binding site is cruel enough or can break down without a ligand molecule during the simulation. Simulation of the dimeric complex for 100 ns was sufficient to compare conformational changes from different MD states. A set of positional and hydrogen bonds were assigned based on the available peptidomimetic structure. Thus, two screenings were conducted with an emphasis on positional constraints or interactions of hydrogen bonds.

Using MDS and MC refinement with Schrodinger and/or YASARA SSP/PSSM methods (Altschul et al., 1997; Hooft et al., 1996a; Hooft et al., 1996b; King and Sternberg, 1996; Krieger et al., 2009; Qiu and Elber, 2006) , each structure was relaxed to the YASARA/Amber force field using knowledge-based potentials within YASARA. The side chains and rotamers were adjusted with knowledge-based potentials, simulated annealing with explicit solvent, and small equilibration simulations using YASARA's refinement protocol (Laskowski RA, 1993) . The entire full-length structure was modeled, filling in any gaps or unresolved portions from the X-ray structure.

Refinement of the finalized models was completed using either Schrodinger's Monte Carlobased module or in-house protocols. These refinements started with generated initial refinement for each independent structure (Altschul et al., 1997; Hooft et al., 1996a; Hooft et al., 1996b; Krieger et al., 2009) . The superposition and subsequent refinement of the overlapping regions yields a complete model for all four proteins. The final structures were subjected to energy optimization with PR conjugate gradient with an R-dependent dielectric.

Atom consistency was checked for all amino acids (and atoms) of the full-length wild-type model, verifying correctness of chain name, dihedrals, angles, torsions, non-bonds, electrostatics, atom typing, and parameters. A multimeric-complex model is predicted, including cofactors and ions. All of the models were exported in the following formats Maestro (MAE), YASARA (PDB). Model manipulation was done with Maestro (Macromodel, version 9.8, Schrodinger, LLC, New York, NY, 2010), or Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) . Analyses were emphasized on the protein-protein interaction regions containing.

Monte Carlo dynamics searching (MC-search) was completed on each model for additional conformational sampling, using methods previously described in the literature (Caulfield and Devkota, 2012; Caulfield, 2011; Caulfield et al., 2011) . Briefly, each protein system was minimized with relaxed restraints using either Steepest Descent or Conjugate Gradient PR, then allowed to undergo the MC search criteria, as shown in the literature (Caulfield and Devkota, 2012; Caulfield, 2011; Caulfield et al., 2011) . The primary purpose of MC, in this scenario, is examining any conformational variability that may occur with different orientations in the region near to protein-protein interfaces.

The total atomic force field was used to minimize the energy of the system, namely, the descent algorithm for 20,000 steps with an iteration interval of 2 fs. The equilibrium of the solvent was carried out using positional restrictions imposed on the atoms of protein structures, while the solvent molecules remained mobile for all 100 ps. Each system was placed in a box in which the layer of the TIP3P water molecule was 10 Å. The final systems were neutralized by the addition of Na + and Cl-ions to a concentration of 150 mM. All simulations were performed under periodic boundary conditions using the V-Rescale Thermostat algorithm to maintain temperature (310 K) and the Parrinello-Rahman Barostat algorithm for constant pressure (1 bar) (Bussi et al., 2007; Parrinello and Rahman, 1981) . Long-range unrelated interactions were calculated using the Particle-Mesh-Ewald (PME) method (Abraham and Gready, 2011). All molecules were relaxed with a molecular dynamics simulation of 100 ns. Ligand topologies were created using the antechamber module from the AmberTools18 package (Case et al., 2005) .

We used SiteMapper (Bhachoo and Beuming, 2017) to identify possible binding sites for docking affinity with the proteins Ace2 (allosteric site), Tmprss2, and M pro . We also used our novel MDS biasing technique algorithm, Maxwell's demon MD, for searching within these sites for potential flexible zones that would have beneficial peptide interactions, which served as a reductive filter limiting the total number of possible sites screened on the proteins to those with adequately deep binding grooves (Caulfield, 2011; Kayode et al., 2016) or interesting insertion sites (Ace2).

Prior to the docking with the Ace2 (allosteric site), Tmprss2, and M pro , we had completed rigorous molecular dynamics simulations (MDS) and Monte Carlo (MC) conformational searching for each model for additional conformational sampling, using methods previously described in the literature (Caulfield and Devkota, 2012; Caulfield, 2011; Caulfield et al., 2011) . The primary purpose of MC, in this scenario, is examining any conformational variability that may occur with different orientations in the region near to protein-protein interfaces.

Over three million compounds were docked to each site using the Glide XP docking program (Bhachoo and Beuming, 2017). All compounds were accounted for using OPLS3 within Maestro program (Maestro-9.4, 2014). Using our published docking protocols on each identified site, we reductively scanned from 100s to the top 10 poses from each docking and then did cross-comparisons of docking scores to retain only the top binding pose of each compound from each site in a winnertakes-all strategy.

Each compound has been converted into a set of energy minimized three-dimensional shapes with the Ligprep module. Without protein preparation, it was used for the correct distribution of protonation and post-minimization in the OPLS3 force field. In the case of assigning restrictions based on ligands (M pro , Tmprss2), we tried to cover the most important and strong interactions. In the case of Ace2, a set of constraints was generated in sufficient quantities to generate combinations of possible interactions. Positional constrains (1.8 A radius) and h-bond constraints were generated in the Schrodinger Glide module, namely in the mesh generation tool. Aromatic and hydrophobic features were represented with short SMARTS. A partial matching protocol for applying constraints has also been used to improve process accuracy. A high throughput screening protocol with regulated ligand flexibility was applied.

Each compound has been converted into a set of energy minimized three-dimensional shapes with the Ligprep module. Without protein preparation, it was used for the correct distribution of protonation and post-minimization in the OPLS3 force field. In the case of assigning restrictions based on ligands (M pro , Tmprss2), we tried to cover the most important and strong interactions. In the case of Ace2, a set of constraints was generated in sufficient quantities to generate combinations of possible interactions. Positional constrains (1.8 A radius) and h-bond constraints were generated in the Schrodinger Glide module, namely in the mesh generation tool. Aromatic and hydrophobic features were represented with short SMARTS. A partial matching protocol for applying constraints has also been used to improve process accuracy. A high throughput screening protocol with regulated ligand flexibility was applied.

Conformations of compound orientations were generated using our standard protocols (Bhachoo and Beuming, 2017; Kalid et al., 2012; Unger et al., 2015) . The starting conformation of relaxed protein structures was first obtained by the method of Polak-Ribière conjugate gradient (PRCG) energy minimization with the Optimized Potentials for Liquid Simulations (OPLS) 2005 force field (Jorgensen, 2004; Jorgensen and Tiradorives, 1988) for 5000 steps, or until the energy difference between subsequent structures was less than 0.001 kJ/mol-Å units. Our docking methodology has been described previously (Caulfield and Devkota, 2012; Friesner et al., 2006; Loving et al., 2009; Vivoli et al., 2012) .

Briefly, compounds were docked within the Schrödinger software suite (Mohamadi et al., 1990 ) using a virtual screening workflow (VSW) (Bhachoo and Beuming, 2017; Friesner et al., 2006; Jacobson et al., 2002; Kalid et al., 2012; Kozakov et al., 2006) . Alternative docking methods were also employed, including in-house software techniques for top leads for SAR elucidation. The top seeded poses were ranked and unfavorable scoring poses were discarded. Top favorable scores from initial dockings yielded hundreds of poses with the top five poses retained. Molecular interactions of the ligand-protein interfaces were used to help determine the optimal binding set, which included descriptors were used to obtain atomic energy terms like hydrogen bond interaction, electrostatic interaction, hydrophobic enclosure and π-π stacking interaction that result during the docking run. Molecular modeling for importing and refining the proteins was completed (Maestro-9.4, 2014).

Examinations of structure stability were examined for all proteins investigated, S-protein:Ace2, Tmprss2, and M pro , respectively (Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield, 2011; Reumers et al., 2005; Schymkowitz et al., 2005; Zhang et al., 2013) . Object stability was used to determine if any changes in structure that were deleterious to function from immediate inspection, which the FoldX algorithm can provide, prior to docking studies. Thus, we examined the local residues around the docking site and determined an electrostatic calculation may be useful to explain the change in function. The molecular model for the full structure and its truncated form are given (Fig. S1 ) using our state of the art methods, which have been established (Abdul-Hay et al., 2013; Ando et al., 2017; Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield, 2011; Caulfield et al., 2011; Caulfield et al., 2014; Caulfield et al., 2015; Fiesel et al., 2015a; Fiesel et al., 2015b; Puschmann et al., 2017; Zhang et al., 2013) .

Local residues within the 12Å cutoff near docking sites were analyzed ( Fig. S1-S2) . Any interactions requiring inducible fit, or Threonine/Serine hydroxyl rotation or other docking parameter (πstacking/halogen-directionality) were also included. Mapping electrostatics was accomplished using the Poisson-Boltzmann calculation for solvation on all amino acids for each docked structure (Caulfield and Devkota, 2012; Caulfield and Medina-Franco, 2011; Caulfield, 2011; Reumers et al., 2005; Schymkowitz et al., 2005; Zhang et al., 2013) E. Libraries used Compounds were derived from either a set of all FDA approved and clinical tested compounds, bioactive set of compounds, or a large multi-million compound set from ZINC database. In the all cases the libraries were prepared using LigPrep described above. The ZINC database was pruned using parameters for better drug-like profile and removal of reactive functional groups and poor chemoinformatics properties delivering a large set suitable for screening on all targets across dynamic time points from MDS.

Supplemental Information can be found online at: XXX ACKNOWLEDGEMENTS AUTHOR CONTRIBUTIONS T.C. designed and conducted most experiments, analyzed data, and wrote the manuscript with inputs from J.M., K.L., M.C. and E.R.; T.C. and M.C. performed MDS, docking, and generated analyses; T.C. and M.C. performed post simulation analyses; C.B. and M.C. performed bioinformatics analysis; T.C. supervised M.C.; T.C. provided expertise on data analysis; E.R. provided expertise and insight interpreting experimental structures and homology models; and T.C. proposed the project to K.L., whom helped with formatting, detailing analysis and edited the manuscript.

The authors declare no competing interest. Zhu, Z., Zhang, Z., Chen, W., Cai, Z., Ge, X., Zhu, H., Jiang, T., Tan, W., and Peng, Y. (2018) . Predicting the receptor-binding domain usage of the coronavirus based on kmer frequency on spike protein. Infect Genet Evol 61, 183-184. Zoete, V., Cuendet Ma Fau -Grosdidier, A., Grosdidier A Fau -Michielin, O., and Michielin, O. (2011) . SwissParam: a fast force field generation tool for small organic molecules. (G) Ligand Interaction Diagram rendered with Maestro for ACE2 with ligand 300 at the allosteric site impacting S-protein binding from SAR-CoV2. This 2D "flat" representation shows the interactions at this particular compounds interface on Ace2 that would interfere with S protein binding. In particular, extending from deeply inserted to superficial, the interactions are described in the subsequent sentences. D382 and D350 are hydrogen bond acceptors (side chains) from the opposite NH+ on the piperazine-like ring deeply inserted into the binding pocket. R393 is a hydrogen bond donor (side chain) to the alcohol group connecting the piperazine-like ring to the fused ring. E37 is a hydrogen bond acceptor (side chain) to one of the NH on the fused ring. The fluorocyclohexane group is entirely solvent-exposed at the mouth of the binding pocket.

Although glycosylation sites at residues N165, N234, N343 from S-protein (PDB code 6VSB), are nearby the ACE2:S-protein binding interface, they do not overlap and interfere with the protein-protein interface, offering an adjacent site is readily available for PPI docking (S-prot glycosylation analysis: DOI: 10.1126/science.abb9983; 10.1101/2020.04.29.069054}. The majority of glycosylation sites are not on the RBD (Fig. S2) , the glycosylation site that is actually present on the RBD, N343, is not in 3D proximity to the binding interface.

Recently, a variant of the S-protein, D614G, was identified to possess enhanced transmissibility and resistance to contemporary interventions and this site is not present on the RBD. Neither the glycosylation sites, nor the enhanced transmissibility variant D614G, are within the 3D proximity to the drug binding site for our targeted protein-protein interface disrupting therapeutics for Ace2. 3-methyl-7-{2-methyl-3-[(1-phenyl-1H-1,2,3,4tetrazol-5-yl) sulfanyl]propyl}-8-(4-methylpiperazin-1-yl)-2,3,6,7-tetrahydro-1H-purine-2,6-dione ACE2 -7.596339 [ (2-chloro-1-benzofuran-3-yl) rac-(3R,4S)-1-[ (3-cyclopentyl-1,2,4-oxadiazol-5yl) methyl]-4- (1-methyl-1H-imidazol-5-yl) 2-(decahydroisoquinolin-2-yl)-N- (2-oxo-2,3-dihydro-1H-1,3-benzodiazol-5-yl) 

Predicting molecular interactions in silico: I. A guide to pharmacophore identification and its applications to drug design

The prince and the pauper. A tale of anticancer targeted agents

A common feature pharmacophore for FDAapproved drugs inhibiting the Ebola virus

Coot: model-building tools for molecular graphics

Features and development of Coot

Optimal docking area: a new method for predicting protein-protein interaction sites

Patho-)physiological relevance of PINK1-dependent ubiquitin phosphorylation

Structural and Functional Impact of Parkinson Disease-Associated Mutations in the E3 Ubiquitin Ligase Parkin

Extra precision glide: docking and scoring incorporating a model of hydrophobic enclosure for protein-ligand complexes

Breakthrough: Chloroquine phosphate has shown apparent efficacy in treatment of COVID-19 associated pneumonia in clinical studies

Clinical and microbiological effect of a combination of hydroxychloroquine and azithromycin in 80 COVID-19 patients with at least a six-day follow up: A pilot observational study

Dose-Dependent Effects of Aloin on the Intestinal Bacterial Community Structure

Development and Experimental Validation of a Docking Strategy for the Generation of Kinase-Targeted Libraries

IL-6: Relevance for immunopathology of SARS-CoV-2

Effects of (-)-epicatechin on molecular modulators of skeletal muscle growth and differentiation

A Trial Comparing Nucleoside Monotherapy with Combination Therapy in HIV-Infected Adults with CD4 Cell Counts from 200 to 500 per Cubic Millimeter

TMPRSS2 and ADAM17 cleave ACE2 differentially and only proteolysis by TMPRSS2 augments entry driven by the severe acute respiratory syndrome coronavirus spike protein

Emodin blocks the SARS coronavirus spike protein and angiotensin-converting enzyme 2 interaction

SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor

The PDBFINDER database: a summary of PDB, DSSP and HSSP information with added value

Errors in protein structures

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

VMD: visual molecular dynamics

Usnic acid

ZINC--a free database of commercially available compounds for virtual screening

Catechin in Human Health and Disease

On the role of the crystal environment in determining protein side-chain conformations

Optimal assignment methods for ligandbased virtual screening

The ReFRAME library as a comprehensive drug repurposing library and its application to the treatment of cryptosporidiosis

PyMod 2.0: improvements in protein sequence-structure analysis and homology modeling within PyMOL

Inhibitors of the Renin-Angiotensin-Aldosterone System and Covid-19

Generalized Scalable Multiple Copy Algorithms for Molecular Dynamics Simulations in NAMD

Structure of M(pro) from SARS-CoV-2 and discovery of its inhibitors

The many roles of computation in drug discovery

The Opls Potential Functions for Proteins -Energy Minimizations for Crystals of Cyclic-Peptides and Crambin

Consensus Induced Fit Docking (cIFD): methodology, validation, and application to the discovery of novel Crm1 inhibitors

An Acrobatic Substrate Metamorphosis Reveals a Requirement for Substrate Conformational Dynamics in Trypsin Proteolysis

Identification and application of the concepts important for accurate and reliable protein secondary structure prediction

Androgen-Induced TMPRSS2 Activates Matriptase and Promotes Extracellular Matrix Degradation, Prostate Cancer Cell Invasion, Tumor Growth, and Metastasis

PIPER: an FFT-based protein docking program with pairwise potentials

Improving physical realism, stereochemistry, and side-chain accuracy in homology modeling: Four approaches that performed well in CASP8

Increasing the precision of comparative models with YASARA NOVA--a self-parameterizing force field

Inference of macromolecular assemblies from crystalline state

Comparison of structure-and ligand-based virtual screening protocols considering hit list complementarity and enrichment factors

Online structure-based screening of purchasable approved drugs and natural compounds: retrospective examples of drug repositioning on cancer targets

Procheck -a Program to Check the Stereochemical Quality of Protein Structures

Optimization of high throughput virtual screening by combining shape-matching and docking methods

Energetic analysis of fragment docking and application to structure-based pharmacophore hypothesis generation

The androgen-regulated protease TMPRSS2 activates a proteolytic cascade involving components of the tumor microenvironment and promotes prostate cancer metastasis

Developing Covid-19 Vaccines at Pandemic Speed

All-atom empirical potential for molecular modeling and dynamics studies of proteins

Molecular dynamics simulation by GROMACS using GUI plugin for PyMOL

Improvements in GROMACS plugin for PyMOL including implicit solvent simulations and displaying results of PCA analysis

Molecular similarity and diversity in chemoinformatics: from theory to applications

Renin-Angiotensin-Aldosterone System Blockers and the Risk of Covid-19

FRED and HYBRID docking performance on standardized datasets

Comparison of topological, shape, and docking methods in virtual screening

Advances in the computational development of DNA methyltransferase inhibitors

Hydroxychloroquine or chloroquine with or without a macrolide for treatment of COVID-19: a multinational registry analysis

Association of Use of Angiotensin-Converting Enzyme Inhibitors and Angiotensin II Receptor Blockers With Testing Positive for Coronavirus Disease

Early treatment of COVID-19 patients with hydroxychloroquine and azithromycin: A retrospective analysis of 1061 cases in Marseille

Macromodel-an integrated software system for modeling organic and bioorganic molecules using molecular mechanics

Simplifying and enhancing the use of PyMOL with horizontal scripts

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

Finding new tricks for old drugs: an efficient route for publicsector drug discovery

Polymorphic transitions in single crystals: A new molecular dynamics method

COVID-19 and Angiotensin-Converting Enzyme Inhibitors and Angiotensin Receptor Blockers: What Is the Evidence?

Coinhibition of the deubiquitinating enzymes, USP14 and UCHL5, with VLX1570 is lethal to ibrutinib-or bortezomib-resistant Waldenstrom macroglobulinemia tumor cells

Scalable molecular dynamics with NAMD

A medicinal chemistry perspective of drug repositioning: Recent advances and challenges in drug discovery

Virtual high throughput screening using combined random forest and flexible docking

Heterozygous PINK1 p.G411S increases risk of Parkinson's disease via a dominant-negative mechanism

SSALN: an alignment algorithm using structure-dependent substitution matrices and gap penalties learned from structurally aligned protein pairs

Homology model-based virtual screening for GPCR ligands using docking and target-biased scoring

SNPeffect: a database mapping molecular phenotypic effects of human non-synonymous coding SNPs

A Large-scale Drug Repositioning Survey for SARS-CoV-2 Antivirals. bioRxiv

Association of Treatment With Hydroxychloroquine or Azithromycin With In-Hospital Mortality in Patients With COVID-19 in

Target-Oriented and Diversity-Oriented Organic Synthesis in Drug Discovery

Prediction of water and metal binding sites and their affinities by using the Fold-X force field

Large compound databases for structure-activity relationships studies in drug discovery

Structural basis of receptor recognition by SARS-CoV-2

Comparative study of rosmarinic acid content in some plants of Labiatae family

Predictors of mortality in hospitalized COVID-19 patients: A systematic review and meta-analysis

Selection of nanobodies that block the enzymatic and cytotoxic activities of the binary Clostridium difficile toxin CDT. Sci Rep 5, 7850. Vaduganathan

Aldosterone System Inhibitors in Patients with Covid-19

Endothelial cell infection and endotheliitis in COVID-19

Free resources to assist structure-based virtual ligand screening experiments

Inhibition of Prohormone Convertases PC1/3 and PC2 by 2,5-Dideoxystreptamine Derivatives

Accelerated bottom-up drug design platform enables the discovery of novel stearoyl-CoA desaturase 1 inhibitors for cancer therapy

Synthesis and evaluation of derivatives of the proteasome deubiquitinase inhibitor b-AP15

Efficacy and safety of glecaprevir/pibrentasvir for chronic hepatitis C virus genotypes 1-6 infection: A systematic review and meta-analysis

SWISS-MODEL: homology modelling of protein structures and complexes

Similarity-based virtual screening using 2D fingerprints

The membrane-anchored serine protease, TMPRSS2, activates PAR-2 in prostate cancer cells

DrugBank: a knowledgebase for drugs, drug actions and drug targets

Genome Composition and Divergence of the Novel Coronavirus (2019-nCoV) Originating in China

Analysis of therapeutic targets for SARS-CoV-2 and discovery of potential drugs by computational methods

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved alpha-ketoamide inhibitors

The dual functions of the extreme N-terminus of TDP-43 in regulating its biological activity and inclusion formation

Functional exhaustion of antiviral lymphocytes in COVID-19 patients

A pneumonia outbreak associated with a new coronavirus of probable bat origin

A Novel Coronavirus from Patients with Pneumonia in China

Fc1c(F)c(F)c(CC(=O)Nc2cccc(c2)C(=O) NC2CC2)c(F)

COc1ccc(CCC(=O)Nc2ccccc2C(=O)OCC (=O)NC(C)c2ccc(F)cc2F)

-methoxyphenyl)methyl]-4H-1,2,4-triazol-3-yl}sulfanyl)-N-(3,5-dimethylphenyl

=O)Nc3cc(C)cc( C)c3)c3ccccc3)n2N

CCOC(=O)Nc1ccc2c(COC(=O)CNS(=O)( =O)c3cccc(Br)c3)cc(=O)

Cc1ccc(c(c1)C(=O)Nc1cccc(OCC(=O)Nc 2ccccc2)c1)-n1cnnc1

O=C(COC(=O)Cc1ccsc1)Nc1cccc2C(=O) c3ccccc3C(=O)

CN(c1ccccc1)S(=O)(=O)c1cccc(c1)C(=O )OC(C(=O)Nc1cccc(c1)C(C)=O)

=O)c(c1)C(F)(F)F

NC(=O)CNC(=O)Nc1cccc(NC(=O)

CCOC(=O)c1c(NC(=O)C(C)Sc2nc3ccc(O CC)cc3[nH]2)sc2CC(C)

-fluorophenyl)-3-methyl-1H-pyrazol-5-yl

Cc1cc(NC(=O)CSc2nnc3ccccn23)n(n1)-c1ccc(F)

FC(F)(F)c1cccc

Cn1cc(nc1C1CCCN(C1)C(=O)Cc1ccc(O)

C(F)(F)F

4-tetrahydropyrimidin-1-yl)-Nmethyl-N-{[(2R,3S)-1-methyl-2-(1-methyl-1H-pyrazol-5-yl)

C(=O)Cn1ccc(=O)[nH]c1=O

3-carbamoylphenyl)-4-(2-fluorophenyl)-1,4-diazepane-1-carboxamide MPRO -7.652268 NC(=O)c1cccc(NC(=O)N2CCCN(CC2)c2 ccccc2F

-(difluoromethyl)-5-methyl-1H-pyrazol-4-yl

)F)c1NC(=O)c1cnc2[nH] c(=O)

COC(=O)CC(NC(=O)C1CC(=O)Nc2cc(F) ccc12)

C=CC(=O)c2ccc3[nH]c(=O)[nH] c3c2)cc1

6-dichlorophenyl)prop-2-enoyl

Clc1cccc(Cl)c1C=CC(=O)c1ccc2[nH]c(= O)

2-yl)methyl]sulfamoyl}phenyl)-3-(2-oxo-1,2,3,4-tetrahydroquinolin-3-yl)propanamide

S(=O)(=O)

-fluorophenyl)-1-(2-methoxyacetyl)-4,5-dihydro-1H-pyrazol-3

Fc1ccc(cc1)S(=O)(=O)Nc1ccccc1C(=O) Nc1cccc(c1)

COc1ccc(cc1)N(C)S(=O)(=O)c1cccc(c1)

C(=O)Nc1cccc(c1)C#N 4-dihydrophthalazine-1-carboxamide MPRO -7.71278 CC(NC(=O)c1nn(Cc2ccccc2)c(=O)c2ccc cc12)c1ccc(cc1)S(N)(=O)=O

-methyl-2,3-dihydro-1H-indol-1-yl

CC1Cc2ccccc2N1C(=O)CSc1cc(C(N)=O)

2-dimethyl-3-oxo-1,2,3,4-tetrahydroquinoxalin-1-yl

COc1ccc(cc1NCC(=O)N1c2ccccc2NC(= O)C1(C)C)S(=O)(=O)

C(=O)NCc2ccccc2)c1

Cn1cnnc1SCC(=O)Nc1cccc(c1)C(=O)N

CN(C)c1ccc(CN(C)C(=O)COc2ccc3NC(= O)CCc3c2)

benzothiazol-2-yl)piperidine-1-carbonyl

Cc1ccccc1NC(=O)Cc1nc(COC(=O)

)=O)c2ccccc2)cs1

NC(=O)NC(C)c1cccc(NC(C)=O)c1)c1 cccc(F)

COc1ccc(OC)c(CNC(=O)NC(C)c2cccc(N C(C)=O)c2)

-ethoxypyridin-3-yl)methyl

CCOc1ncccc1CNC(=O)Cc1sc(C)nc1-c1ccc(C)

3-dihydro-1-benzofuran-2-yl)methyl]-2-oxo-1,2,3,4-tetrahydroquinoline-6-carboxamide MPRO -7.632513 O=C(NCC1Cc2ccccc2O1)c1ccc2NC(=O)

yl)acetyl]pyrrolidin-2-yl}-4-methyl-4,5-dihydro-1H-1,2,4-triazol-5-one

4-dichlorophenyl)methyl]-N-({4-oxo-4H-pyrido[1,2-a]pyrimidin-2-yl}methyl

Clc1ccc(CN2CCC3(C2)CCCN(C3)C(=O)N Cc2cc(=O)n3ccccc3n2)

CCC(=O)Nc1ccc(cc1)C(C)NC(=O)

CCC(=O)Nc1ccc(cc1)C(C)NC(=O)

NC(=O)c1ccn(n1)-c1cccc(NC(=O)C2CCN(CC2)S(=O)(=O)C =Cc2ccccc2)

CC(N1CCN(CC=Cc2ccccc2)CC1)C(=O)N c1cccc(c1)C#N 4,5-tetrahydro-1H-1

CCCn1c(SCC(=O)N2C(C)CC(=O)Nc3ccc cc23

-oxo-1,2-dihydropyridin-1-yl)acetamide MPRO -7.623712 FC(F)(F)c1cc(NC(=O)Cn2ccccc2=O)cc(c

C(F)(F)F

3-dihydro-1,4-benzodioxin-6-yl)ethyl

Fc1cccc(c1)N(CCN1C(=O)c2cccc3cccc( C1=O)c23)C(=O)

10-dihydroanthracen-1-yl)-2-(thiophen-2-yl)quinoline-4-carboxamide MPRO -7.800731 O=C(Nc1cccc2C(=O)c3ccccc3C(=O)c12 )

C(=O)c3ccccc3C(=O)

COc1ccc(NS(=O)(=O)c2cccc(c2)

COc1ccc(cc1)S(=O)(=O)NN=Cc1ccc(OC (=O)Nc2ccccc2)c(OC)

3-dimethylcyclohexyl)carbamoyl]methyl 5-(benzylsulfamoyl

S(=O)(=O)NCc2ccccc2)