Virtual Screening for Inhibitors of SARS-CoV-2 Entry by Simultaneous Targeting of TMPRSS2 and CTSL

The global impact of the COVID-19 pandemic continues, driven by the emergence of various SARS-CoV-2 variants and the limited availability of effective therapeutic options. The viral entry process requires the SARS-CoV-2 spike protein to be activated by host cell proteases, particularly TMPRSS2 and CTSL. These proteases facilitate viral membrane fusion and support the endocytic uptake of the virus. In this study, a computational approach using virtual screening was applied to discover compounds that could simultaneously inhibit both TMPRSS2 and CTSL. Two pharmacophore models were constructed from the binding pockets of these proteases, each in complex with its respective ligand. These models were employed to screen a library containing 41,775 compounds, including 10,849 drugs from the ChEMBL database and 30,926 natural products from the NPASS database. This process identified 115 compounds—54 pharmaceutical drugs and 61 natural products—that matched both the TMPRSS2 and CTSL pharmacophore models. The identified compounds were then docked into the protease structures to further refine the list. Molecular docking simulations identified 17 top candidates (5 drugs and 12 natural products) with stronger binding energies than the reference ligands and known inhibitors of these proteases. These candidates were then evaluated using various filters, including ADMET predictions, drug-likeness assessments, and synthetic accessibility. Among the drugs, silibinin emerged as the leading repurposed candidate, showing promise as a dual inhibitor for SARS-CoV-2 treatment. Additionally, the natural product barettin was highlighted as a strong contender for development into a novel dual-target inhibitor of TMPRSS2 and CTSL.

Coronaviruses (CoVs) have been behind several major global health crises due to their ability to cause severe illness and death in humans. This includes the 2002-2004 outbreak of severe acute respiratory syndrome (SARS) and the 2012 outbreak of Middle East respiratory syndrome (MERS). In March 2020, the World Health Organization (WHO) declared COVID-19, caused by SARS-CoV-2, a global pandemic due to its rapid spread. By October 2022, COVID-19 had affected more than 617 million individuals worldwide, resulting in over 6 million deaths [1]. While mass vaccination and convalescent plasma transfusions have played key roles, significant focus has been on drug repurposing to find effective treatments for COVID-19. This approach has led to the investigation of existing and experimental drugs, such as remdesivir, nirmatrelvir/ritonavir, molnupiravir, baricitinib [2], and ivermectin [3], as possible therapies for those infected. Despite the large number of ongoing studies, the hunt for potent and reliable therapeutic agents against SARS-CoV-2 remains a critical priority.

CoVs belong to a diverse group of viruses characterized by a single-stranded RNA genome and an envelope, and are responsible for causing respiratory diseases in both mammals and birds. A hallmark feature of these viruses is the spike (S) protein, which is essential for their ability to bind to and enter host cells [4]. For SARS-CoV-1 and SARS-CoV-2, the S protein interacts with the angiotensin-converting enzyme 2 (ACE2) receptor on type II pneumocytes to facilitate viral entry. The S protein consists of two functional subunits: the S1 subunit, containing the receptor-binding domain (RBD) that binds ACE2, and the S2 subunit, which is responsible for mediating the fusion of the viral and host cell membranes [5]. Once the virus binds to ACE2, it can enter the host cell either through direct membrane fusion or via receptor-mediated endocytosis involving endosomes. Host cell proteases such as transmembrane protease serine 2 (TMPRSS2) and cathepsin L (CTSL) play key roles in activating the S protein, allowing viral entry via membrane fusion and endosomal routes, respectively [6]. Simultaneously blocking both protease-mediated pathways is a promising strategy for preventing SARS-CoV-2 entry and could serve as an effective early intervention or prophylactic measure.

While the precise function of TMPRSS2 remains unclear, its overexpression has been associated with prostate cancer, and it is known to be essential for the priming of SARS-CoV-1 and MERS-CoV S proteins [7]. CTSL, on the other hand, is a cysteine protease that helps maintain cellular balance and is implicated in various cancers. It is overexpressed during COVID-19 infection and correlates with disease severity [8]. Given that both TMPRSS2 and CTSL contribute to viral entry in cell-type-specific ways, there is growing interest in identifying drugs that target both host proteases. Drugs that target two critical cellular enzymes are more likely to be potent and less prone to developing resistance. Moreover, targeting host rather than viral proteins reduces the risk of resistance due to viral mutations [9].

In contemporary drug discovery, computational techniques such as pharmacophore modeling and molecular docking are widely used to efficiently identify potential drug leads. A pharmacophore represents a set of spatially arranged chemical features necessary for ligand-receptor interaction [10]. However, while pharmacophore-based screening identifies potential hits based on structural alignment, it does not always predict binding affinity. Therefore, molecular docking is employed to predict the strength of interaction between the ligand and receptor, thereby helping refine and validate potential drug candidates. These methods are often used in combination as a part of virtual screening approaches, providing a cost-effective and efficient means of identifying promising compounds for drug development [11]. Our prior studies have demonstrated the successful application of pharmacophore modeling and docking to discover potential drug leads against M. tuberculosis targets [12-16].

This study aims to identify existing therapeutic drugs and natural compounds that can simultaneously target both TMPRSS2 and CTSL. To do so, two distinct pharmacophore models were developed using the crystal structures of these proteases bound to their respective ligands. These models were then used to screen compounds from the ChEMBL and NPASS databases. Common hits from both screenings were further evaluated through molecular docking to identify the best candidates. The top-performing compounds were analyzed using ADMET, drug-likeness, and synthetic accessibility filters to identify drugs and natural products with potential for repurposing to treat SARS-CoV-2 infections.

Computational analyses were carried out on a Windows 11 machine equipped with an Intel® Core™ i5-10300H processor (2.50 GHz, 8 cores), 8 GB of RAM, and an NVIDIA GeForce GTX 1650 Ti GPU.

Development of pharmacophore models

To design structure-based pharmacophore models, crystal structures of the human TMPRSS2 enzyme in complex with nafamostat (PDB ID: 7MEQ; resolution 1.95 Å) and CTSL bound to the dipeptidyl nitrile inhibitor NOW (PDB ID: 3HHA; resolution 1.27 Å) were employed. For the TMPRSS2 structure (7MEQ), the ligand GBS was restored to its original nafamostat form, as seen in the complex with TMPRSS2. The CTSL structure (3HHA) was directly imported into the LigandScout 4.4 software using the PDB search tool [17]. The active binding sites for both receptors were automatically detected through the co-crystallized ligands, and optimization was performed using the MMFF94 force field. Subsequently, pharmacophore models were constructed based on the chemical features of the binding pockets that interact with these ligands.

Pharmacophore model evaluation

To assess the ability of the pharmacophore models to distinguish active compounds from decoys, experimentally validated active ligands for TMPRSS2 and CTSL were identified from scientific literature and uploaded in SMILES format to the DUDE•E repository (http://dude.docking.org/). These ligands, along with decoy compounds, were converted into database files (.ldb format) using LigandScout 4.4. The validation process involved screening both active and decoy libraries using the pharmacophore models. The performance of these models was evaluated through receiver operating characteristic (ROC) curve analysis, calculating the area under the curve (AUC) and enrichment factor (EF) to gauge their predictive accuracy.

Database compilation

A compound library comprising 41,775 structures was assembled, with 10,849 small-molecule drugs sourced from ChEMBL (https://www.ebi.ac.uk/chembl/) and 30,926 natural products obtained from NPASS (http://bidd.group/NPASS/). These molecules were processed and optimized in OpenBabel 2.3.1 [18], including conversion to 3D structures, addition of polar hydrogens, assignment of MMFF94 partial charges, and removal of duplicates. After this optimization, the library consisted of 40,039 unique compounds, including 9,945 drugs and 30,094 natural products. These structures were imported into LigandScout 4.4, where they were subjected to MMFF94 minimization and organized into an .ldb database file. To prepare for further screening, each compound was processed to generate up to 25 conformations per molecule.

Screening of compound collections

The compound library was evaluated using the two pharmacophore models under optimized conditions derived from the validation phase. The selected scoring approach was pharmacophore fit, and the screening mode was set to ensure that all essential features of the query were matched, with a tolerance of at most 1 feature omission, while incorporating exclusion volumes for accuracy. The first stage of screening targeted the TMPRSS2 pharmacophore, followed by a secondary screening against the CTSL model. The compounds that passed these screens were organized by average pharmacophore-fit score and subsequently exported to Excel for further examination. Overlapping compounds from both models were manually curated and saved in the mol2 file format for subsequent use.

Molecular docking process

For the preparation of the TMPRSS2 (PDB ID: 7MEQ) and CTSL (PDB ID: 3HHA) structures, missing residues were supplemented using homology modeling through SWISS-MODEL. The structures were then refined using UCSF Chimera 1.17 and the Dock Prep tool, which involved removing heteroatoms and solvent molecules, adding polar hydrogens, and assigning Gasteiger charges with Amber’s Antechamber module. Protein minimization was carried out using the MMFF94 force field in OpenBabel 2.3.1. Before docking the hits from virtual screening, the AutoDock Vina protocol [19] was validated by re-docking the co-crystallized ligands, GBS for TMPRSS2 and NOW for CTSL, into their binding pockets after optimization. The grid box for TMPRSS2 was set to 25 Å × 25 Å × 25 Å, centered at coordinates (-8.55, -2.05, 16.07), and for CTSL, it was set to 25 Å × 25 Å × 25 Å, centered at (8.65, 8.06, -2.91). Following validation, the compounds from the virtual screening were added to UCSF Chimera 1.17 as mol2 files, optimized, and docked into the active sites of TMPRSS2 and CTSL. The best docking poses for each compound were evaluated, and their binding energies were recorded. To further validate the results, known inhibitors such as camostat (for TMPRSS2) and CLIK-148 (for CTSL) were docked along with the co-crystallized ligands. The resulting receptor-ligand complexes were saved in pdb format, and interaction diagrams were generated for analysis in Discovery Studio Client 2021.

ADMET profiling

The ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) properties of the compounds were predicted using computational tools such as SwissADME (http://www.swissadme.ch/) and OSIRIS Property Explorer (https://www.organic-chemistry.org/prog/peo/). The evaluation focused on key properties such as water solubility, gastrointestinal absorption, blood-brain barrier penetration, P-glycoprotein (P-gp) binding, and inhibition of cytochrome P450 2D6. Toxicological risks, including mutagenic, tumorigenic, irritant, and reproductive effects, were also assessed. The drug-likeness and lead-likeness of the compounds were evaluated using Lipinski’s rule of five, along with synthetic accessibility (SAscore) calculated in SwissADME.

In the search for new scaffolds that could effectively interact with both TMPRSS2 and CTSL, we developed two distinct pharmacophore models using the available crystallographic data of these enzymes bound to their respective ligands (Figure 1). During the creation of the pharmacophore models, exclusion volumes were included to preserve the steric features and overall shape of the binding pockets for both targets. For the TMPRSS2 model, the key features consisted of 2 hydrogen bond donors, 1 hydrogen bond acceptor, 2 hydrophobic interactions, 1 positively charged ionizable site, and 18 exclusion volumes. In contrast, the CTSL model incorporated three hydrogen bond donors, three hydrogen bond acceptors, three hydrophobic interactions, one residue-binding site, and 20 exclusion volumes.

Before screening large compound databases, it’s essential to validate pharmacophore models to ensure their quality and reliability, particularly after selecting features and adjusting tolerances. In this study, the structure-based pharmacophore models for TMPRSS2 and CTSL were assessed for their ability to differentiate active from decoy compounds. Each model’s test set included five active compounds identified in the literature and 250 decoy compounds sourced from the DUDE•E directory. These decoys were designed with similar physical properties to the active compounds (e.g., molecular weight, logP, HBD, HBA, and the number of rotatable bonds) but had different structural topologies to minimize binding similarities [20].

To assess the models’ performance, the receiver operating characteristic (ROC) curve was used. This curve compares the true positive rate (sensitivity) to the false positive rate (1-specificity), providing a visual representation of the model’s accuracy. An ideal model would produce a curve closer to the upper-left corner, while poor performance results in a diagonal line [21]. The performance of the models was summarized using two metrics: the Area Under the Curve (AUC) and the Enrichment Factor (EF). The AUC, ranging from 0 to 1, represents the model’s ability to distinguish between active and inactive compounds, with higher values indicating better performance. The EF shows how many more active compounds are expected to appear in the top-ranked portion of the screening compared to a random model. EF values can range from 1 (random screening) to values greater than 100, indicating the model’s effectiveness in identifying true positives [22].

Both models successfully identified all five active compounds in the test set. Among the 250 decoys, the TMPRSS2 model had only three false positives, while the CTSL model had 15. At the 1% threshold, both pharmacophore models exhibited AUCs of 1.00 and EFs of 25.5, demonstrating their excellent ability to discriminate active compounds. As shown in Figure 2, an EF value of 25.5 indicates that the top 1% of screening results would likely contain 25 times as many active compounds as a random selection from the database.

Pharmacophore-based screening of the ligand library resulted in the identification of 115 compounds (comprising 44 drugs and 61 natural products) that were found in both the TMPRSS2 and CTSL pharmacophore screens. The fit scores for these overlapping compounds ranged from 42.08 to 47.08 for TMPRSS2 and from 71.31 to 72.04 for CTSL. These scores assess the degree to which ligand substructures align with the pharmacophore model features, with higher values indicating a stronger fit with the target proteins.

To assess ligand binding affinity and pose, molecular docking was performed, focusing on the receptor-binding pockets. Before docking the overlapping compounds from the pharmacophore screening, the docking protocol’s performance was validated. Ligands from the crystal structures of TMPRSS2 and CTSL were redocked into their respective active sites, and AutoDock Vina returned poses that aligned closely with their original co-crystallized structures. The RMSD metric, used to quantify structural similarity between superimposed molecules, is generally considered accurate when it falls below 1.5 Å [23]. As shown in Figure 3, the RMSD values were 0.367 Å for the GBS ligand on TMPRSS2 and 0.850 Å for the NOW ligand on CTSL, both below the 1.5 Å threshold, indicating high structural similarity and reliable binding pose predictions for both proteins.

The re-docking of the GBS ligand to TMPRSS2 preserved several significant interactions, including hydrogen bonds with key residues such as ASP435 (substrate-binding) and SER441 (catalytic) (Figure 1). GBS, which represents the phenylguanidino acyl group of nafamostat, slowly hydrolyzes upon contact with the enzyme, resulting in the absence of the other catalytic triad residues, HIS296 and ASP345, in the re-docked structure. Despite this, GBS still interacted with GLY462 and additional residues, such as SER436, GLY439, and GLY464, in TMPRSS2.

Meanwhile, the NOW ligand, upon re-docking, retained hydrogen bonds with the catalytic dyad CYS25 and HIS163, as observed in the experimental structure. Moreover, it maintained interactions with several hydrophobic and electrostatically interacting residues within the substrate-binding pocket, including ASP162, MET161, ALA135, ALA214, LEU69, MET70, and GLY68.

After fine-tuning the docking grid parameters and confirming the protein preparation and AutoDock Vina protocols, the 115 pharmacophore hits from the virtual screening were docked into the active sites of TMPRSS2 and CTSL. The known inhibitors camostat (TMPRSS2) and CLIK-148 (CTSL) were also included as controls, along with the co-crystallized ligands. Hits were identified based on their binding energies, where compounds with more negative values than the controls were considered favorable, indicating a more stable receptor-ligand interaction, with the negative sign denoting a spontaneous binding process [24]. For TMPRSS2, the binding energies of the re-docked GBS and camostat were -7.6 and -7.5 kcal/mol, respectively. For CTSL, the binding energies for NOW and CLIK-148 were -7.0 and -6.7 kcal/mol, respectively.

Out of the 54 drugs and 61 natural products tested, 5 drugs and 12 natural products exhibited binding energies equal to or more negative than the known inhibitors and co-crystallized ligands. Consequently, 17 compounds from the pharmacophore-based hits passed both screening tests and showed potential as dual inhibitors of TMPRSS2 and CTSL, essential targets for hindering SARS-CoV-2 cell entry. The binding energies of these potential inhibitors ranged from -7.6 to -9.1 kcal/mol for TMPRSS2 and -7.0 to -8.7 kcal/mol for CTSL (Table 1).

In silico evaluation of the pharmacokinetic properties, medicinal chemistry profiles, and toxicity of the top compounds is a critical step in identifying promising candidates for further optimization or progression in the drug development process. The compounds identified in the docking analysis were assessed using the SwissADME server to evaluate their potential as drug-like molecules, and the OSIRIS Property Explorer was used to assess their mutagenic, tumorigenic, irritant, and reproductive toxicity risks (Table 2). Additionally, some compounds that had already been tested in clinical trials were included in the screening to verify their drug-likeness and safety profiles.

^aSynthetic Accessibility score of 6.0 was set as a threshold between easy- and hard-to-synthesize compounds [25].

Among the five ChEMBL drugs selected, nearly all demonstrated favorable pharmacokinetic profiles, although they exhibited limited gastrointestinal absorption, with only three passing the Lipinski drug-likeness filter. Pilaralisib, Cloguanamil, and Carubicin were identified as having moderate to high toxicity risks related to mutagenic, irritant, or reproductive effects. In contrast, PCI-27483 and Carubicin failed to meet the Lipinski criteria, which diminishes their potential for further pharmaceutical development. Silibinin, however, emerged as the most promising candidate. It displayed oral bioavailability, despite lower gastrointestinal absorption, and did not cross the blood-brain barrier nor interact with P-gp or CYP2D6. Additionally, it adhered to the Lipinski rules for drug-likeness, maintained an acceptable SAscore, and showed no signs of toxicity, making it a strong candidate for repurposing as a COVID-19 treatment, particularly by inhibiting SARS-CoV-2 cell entry through dual targeting of TMPRSS2 and CTSL.

Interaction diagrams for Silibinin (Figure 4) revealed key interactions within the TMPRSS2 active site, including hydrogen bonds with SER441 and van der Waals contacts with HIS296, both of which are part of the catalytic triad. Furthermore, silibinin’s ether and hydroxyl groups formed conventional hydrogen bonds with residues LYS342, CYS437, and GLY462. However, interactions with ASP440, characterized as donor-donor bonds, introduced repulsive forces that could weaken the complex’s stability. In the CTSL binding site, silibinin established interactions with the catalytic dyad, including pi-sulfur bonding with CYS25 and weak C-H bonding with HIS163. Additional hydrogen bonds were formed with ASN66, GLN21, TRP189, and GLY164, complemented by hydrophobic interactions with MET70 in the enzyme’s binding pocket.

Silibinin is the principal active component in silymarin, an extract derived from the seeds of Silybum marianum (milk thistle). It is approved in Europe for treating conditions such as hepatotoxicity, chronic hepatitis, and cirrhosis. Additionally, silibinin is being investigated for its potential as an anticancer agent due to its inhibitory effects on STAT3 signaling, protection against environmental toxins and UV-induced damage, and anti-inflammatory and immunomodulatory properties [26]. Recently, silibinin has been proposed as a therapeutic option for COVID-19, owing to its dual function: inhibiting the host cytokine storm via STAT3 suppression and blocking the RdRp enzyme required for viral replication. Clinical trials are currently underway in Spain to assess silibinin’s efficacy in patients with oncohematological conditions and COVID-19 [27]. Previous docking studies have shown that silibinin forms stable complexes with the SARS-CoV-2 spike protein RBD and Mpro [28], suggesting its potential to disrupt multiple stages of SARS-CoV-2 pathogenesis, including viral entry, replication, post-translational modifications, and immune regulation.

In this study, a natural product database was also explored to identify potential agents for early-stage or prophylactic COVID-19 treatments. However, all 12 selected natural compounds exhibited poor drug-likeness or undesirable ADME profiles, with four showing significant toxicity risks, such as mutagenic, irritant, or reproductive effects. Nevertheless, NPC304187 (Barettin) demonstrated the most favorable pharmacokinetic and lead-likeness characteristics. While Barettin is a P-gp substrate, it does not pose any toxicity concerns, making it the most promising natural product candidate for further development as a dual inhibitor targeting TMPRSS2 and CTSL.

Based on its interaction analysis (Figure 5), Barettin’s guanidine group formed multiple hydrogen bonds and electrostatic interactions with key residues in the TMPRSS2 binding site, including ASP435, SER436, GLY464, and GLY462, in addition to van der Waals interactions with several hydrophobic residues. The pyrazinol moiety interacted with the catalytic residue HIS296 through hydrogen bonding, while the bromoindole moiety formed van der Waals interactions with SER441. In the CTSL binding pocket, Barettin’s pyrazinol moiety exhibited hydrogen bonding with CYS25 and MET161, as well as van der Waals interactions with HIS163. The guanidine moiety also interacted with ASP114 and ASP71, resulting in attractive electrostatic interactions. Additionally, the bromoindole moiety showed amide-pi stacking with GLY67 and hydrogen bonding with ASN66, with additional van der Waals interactions occurring in the enzyme’s hydrophobic pocket.

Barettin, the primary bioactive compound from the deep-sea sponge Geodia barretti in the North Atlantic, has demonstrated antioxidant and anti-inflammatory properties by inhibiting protein kinases, including RIPK2 and CAMK1α [29]. However, its antiviral properties had not been previously identified. This compound warrants further investigation using in silico approaches, followed by in vitro and in vivo testing to evaluate its potential effectiveness against SARS-CoV-2.

In this investigation, we used an innovative in silico method to identify compounds targeting both TMPRSS2 and CTSL, which are pivotal in the SARS-CoV-2 entry process. Through a combination of structure-based pharmacophore modeling, molecular docking, virtual screening, and ADMET predictions, we identified silibinin, a Phase 3 clinical candidate, and Barettin, a natural product, as promising dual inhibitors of these enzymes. The results of this study pave the way for future experimental and clinical research aimed at developing novel therapeutic strategies against COVID-19, leveraging the molecular features of the top compounds identified.

The authors express gratitude to Inte: Ligand for providing access to the LigandScout software and to The Scripps Research Institute for providing AutoDock Vina. ILSV also thanks the Department of Science and Technology – Science Education Institute for their generous scholarship.

None

None

None

Sanjay Kulkarni, Meenal Joshi & Aniket Deshmukh contributed to this work.

Authors and affiliations

Department of Pharmacology and Drug Design, Faculty of Pharmacy, Savitribai Phule Pune University, Pune, India
Sanjay Kulkarni & Meenal Joshi

Department of Toxicology and Safety Assessment, Faculty of Medicine, IIT Bombay, Mumbai, India
Aniket Deshmukh

Corresponding author

Correspondence to Meenal Joshi

Open Access The author(s) retain copyright. This article is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. It may be shared and adapted for non-commercial purposes with appropriate attribution, an indication of changes, and distribution of adaptations under the same license. Third-party material may be subject to separate terms identified in its credit line. View the license at https://creativecommons.org/licenses/by-nc-sa/4.0/.