Computational screening of camostat and related compounds against human TMPRSS2: A potential treatment of COVID-19.

The global coronavirus pandemic has burdened the human population with mass fatalities and disastrous socio-economic consequences. The frequent occurrence of these new variants has fueled the already prevailing challenge. There is still a necessity for highly effective small molecular agents to prevent severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection. Here, we targeted the human transmembrane surface protease TMPRSS2, which is essential for proteolytic activation of SARS-CoV-2. Camostat is a well-known inhibitor of serine proteases and an effective TMPRSS2 inhibitor. A virtual library of camostat-like compounds was computationally screened against the catalytic site of TMPRSS2. Following a sequential in-depth molecular docking and dynamics simulation, we report the compounds that exhibited promising efficacy against TMPRSS2. The molecular docking and MM/PBSA free energy calculation study indicates these compounds carry excellent binding affinity against TMPRSS2 and found them more effective than camostat. The study will open doors for the effective treatment of coronavirus disease 2019.

Introduction

Coronavirus disease (COVID-19) has become a life-threatening pandemic. The lack of effective and adequate care and high-mortality rate have motivated researchers to design effective COVID-19 prevention strategies and vaccines. COVID-19 incidents are rising rapidly, with over 266 million positive cases and over 5.26 million deaths worldwide by the second week of May 2021 (Organization, 2021). The invasion of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) into the host lung epithelial cells is mediated by the binding of its transmembrane spike glycoprotein with the host angiotensin-converting enzyme-2 (ACE-2) receptors (Ni et al., 2020, Zamorano Cuervo and Grandvaux, 2020). SARS-CoV-2 then uses the human transmembrane surface protease TMPRSS2 to cleave and trigger the spike protein, allowing the virus to participate in host membrane fusion (Glowacka et al., 2011, Huang et al., 2020, Tang et al., 2020).

Usually expressed in epithelial cells, TMPRSS2 is a transmembrane serine protease (Chen et al., 2010). The TMPRSS2 extracellular protease domain can cleave a spike protein domain to instigate membrane fusion. This protease promotes the entrance of various viruses into cells, including influenza, SARS, and the Middle East respiratory syndrome (Glowacka et al., 2011, Hoffmann et al., 2020). Both TMPRSS2 and furin are proteases at the cell surface and are crucial for the proteolytic activation of SARS-CoV-2 in human airway cells. While furin cleaves the spike protein in the S1/S2 region, TMPRSS2 cleaves it at the 2′ site, triggering spike protein membrane fusion activity (Bestle et al., 2020). The crucial role of TMPRSS2 in the viral life cycle has attracted it to be considered a potential target, restricting the viral-host cell entry. Several TMPRSS2 inhibitors have demonstrated successful (in vitro) SARS-CoV-2 infection suppression, suggesting that TMPRSS2 is a mediating factor of viral entry (Padmanabhan et al., 2020, Shang et al., 2020). TMPRSS2 further weakens the detection of viruses by neutralizing antibodies from the host, thereby facilitating viral pathogenesis (Glowacka et al., 2011). The occurrence of new variants with critical mutations that increase their resistance toward antiviral and neutralizing antibodies is a growing concern (Garcia-Beltran et al., 2021, Resende et al., 2021). The frequent occurrence of these new variants has challenged the antiviral drug discovery process (Hoffmann et al., 2021a, Hoffmann et al., 2021b, Shen et al., 2021, Wang et al., 2021). With its vital role in SARS-CoV-2 pathogenesis, several studies recommend TMPRSS2 as an effective target to suppress SARS-CoV-2 infection (Padmanabhan et al., 2020, Ragia and Manolopoulos, 2020).

Camostat is a well-known serine protease inhibitor and an effective TMPRSS2 inhibitor, which may be a plausible antiviral against SARS-CoV-2 (Fig. 1) (Breining et al., 2021a, Breining et al., 2021b). Using camostat to block virus-membrane fusion can reduce viral infection by two-thirds (Breining et al., 2021a, Breining et al., 2021b). In Japan, camostat mesylate is prescribed to treat chronic pancreatitis and drug-induced lung injury (Zhou et al., 2015). Camostat is currently in clinical trials as a COVID-19 treatment option (CTID: NCT04662086, NCT04455815, NCT04662073, and NCT04657497) (Breining et al., 2021a, Breining et al., 2021b, Uno, 2020).

Fig. 1

Camostat structure.

This study intended to identify potent inhibitors targeting the TMPRSS2 catalytic site. In this regard, compounds analogous to camostat were retrieved from the Pub-Chem database. Several in silico evaluations (Fig. 2), comprising virtual screening, molecular dynamics simulation, and free energy calculation studies, have identified four compounds with high binding efficacy against TMPRSS2. The selected compounds demonstrated better and stable binding affinity, significantly improving over the known inhibitor, camostat.

Fig. 2

Binding of SARS-CoV-2 protein to human ACE-2 and its priming by TMPRSS2 for fusion and internalization with ACE-2 spike complex. Inhibition of TMPRSS2 by camostat and workflow for identifying camostat analogs with improved inhibitory potentials against TMPRSS2.

Materials and methods

Protein structure preparation

The protein ID: 7MEQ was used to collect the structure information of the TMPRSS2 protein from well-known structure database at RCSB (). The heteroatoms were removed before optimization for energy utilizing the steepest descent algorithm for 5000 steps in UCSF’s chimera software (Pettersen et al., 2004a, Pettersen et al., 2004b). Docking studies were carried on the optimized and validated structure.

Ligand similarity search and library preparation

The structures of camostat (CID: 2536) and related analogous compounds were extracted in the SDF format from an open chemistry PubChem database (Kim et al., 2016). A total of 223 compounds were retrieved, and those with camostat were energy-minimized. Energy minimization was performed using the UCSF Chimera molecular modeling package (Pettersen et al., 2004a, Pettersen et al., 2004b) for 5,000 steps with the steepest-descent method. The minimized compound structures were then employed for virtual screening.

Virtual screening

CCDC GOLD was utilized to screen all the compounds based on their binding affinity within the active site of TMPRSS2 (Jones et al., 1997). The crystal ligand was selected as a reference for assigning the active binding site. An area of 8 Å was used for generating the grid surface. For each compound, a total of 100 possible conformations were generated, and the best pose was selected using the ChemPLP score. The top-scoring compounds were further screened for their ADMET properties using the ADMETlab 2.0 webserver (Xiong et al., 2021). The compounds passing the ADMET test were finally selected.

Molecular dynamics simulation

The selected complexes underwent dynamic simulations to investigate the stability of camostat and other top-scoring compounds in complex with TMPRSS2. The complexes of camostat and other selected compounds were prepared using molecular docking. The molecular dynamics simulations were performed using the GROMACS 2020 package with Charmm27 force field (Bjelkmar et al., 2010, Hess et al., 2008, Pronk et al., 2013). GROMACS is widely used in molecular dynamics and protein–ligand simulation studies (Baig et al., 2016, Baig et al., 2019, Baig et al., 2014, Bao et al., 2018, Liu et al., 2018). The selected complexes were solvated within the dodecahedron box water model with a box wall and solute margin set to 0.1 nm. The system was neutralized by adding Na+/Cl− counterions (Mark and Nilsson, 2001, Mark and Nilsson, 2002). The long-range coulombic interactions (Darden et al., 1993) were estimated by particle-mesh Ewald method and for van der Waals interactions we used Lennard–Jones method at a cutoff distance of 0.1 nm. The Linear Constraint Solver method was utilized to constrain the bond lengths. The time step was set to 0.002 ps (Hess, 2008, Hess et al., 1997). The system build was steepest-descent energy minimized for 10,000 steps and subjected to equilibration for 1 ns. Berendsen weak coupling systems was set on to maintain biological simutaion framework at 300 K temperature and 1 bar pressure (Berendsen et al., 1984, Izaguirre et al., 2001). To the stable equilibrated system generated production run was performed for 100-ns. PyMol and xmgrace were used for graphical inspections, analysis and graph generation ().

Molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) free energy calculations

The MM-PBSA approach plays a more efficient role in drug discovery than the traditional free energy calculations (Kollman et al., 2000). The binding free energy was estimated by considering the vacuum potential energy and solvation free energy (polar and nonpolar). Poisson–Boltzmann equation and solvent accessible surface area methods were used to calculate polar and nonpolar solvation energies (Rizzo et al., 2006, Still et al., 1990). The Poisson–Boltzmann equation approximates the electrostatic component of biological macromolecules and helps study the ligand-binding affinity to the protein. The SASA method helps identify the surface surrounding the protein with van der Waals contact probed by the solvent sphere. The MMPBSA.py module was used to for MM-PBSA calculations at AMBER interface (Miller et al., 2012).

Binding free energy (ΔGbinding) were calculated as per the following equations:

where

and

Results and discussion

TMPRSS2, a 70-kDa protein, is a serine protease mediating the entry of SARS-CoV-2 via the ACE-2 enzyme (Hoffmann et al., 2020). The binding of the SARS-CoV-2 S protein to ACE-2 is an essential step required for cellular entry. TMPRSS2 primes this binding, thereby promoting the endocytic entry of the virus (Hoffmann et al., 2020). The emergence of novel SARS-CoV-2 variants viz. delta, kappa, and epsilon has threatened the effectiveness of vaccines (Tian et al., 2021). The mutations carried by these novel variants increase viral transmission and immune escape (Kannan et al., 2021, Raheem et al., 2021). Studies have suggested that targeting this serine protease (TMPRSS2) may be a crucial checkpoint for controlling the viral entry of SARS-CoV-2 within human cells (Glowacka et al., 2011, Hoffmann et al., 2020, Iwata-Yoshikawa et al., 2019, Kawase et al., 2012, Zhou et al., 2015). Thus, the selection of TMPRSS2 as a therapeutic target holds a significant scope for successfully treating SARS-CoV-2 infection. Camostat mesylate, a protease inhibitor, is widely reported to be capable of blocking the virus-activating cellular protease TMPRSS2 and thereby inhibiting the SARS-CoV-2 infection (Hoffmann et al., 2020). This clinically proven protease inhibitor approved for human use in Japan as for the treatment of pancreatitis (Abe, 1980, Ohshio et al., 1989), is being widely investigated as a COVID-19 treatment option (Breining et al., 2021a, Breining et al., 2021b, Hoffmann et al., 2021a, Hoffmann et al., 2021b, Hofmann-Winkler et al., 2020). The search for other such molecules capable of blocking TMPRSS2 may open new therapeutic gateway for the treatment of COVID-19.

Molecular docking-based virtual screening and ADMET analysis

Recently reported structure of TMPRSS2 with nafamostat was used for docking studies. The binding affinity of the studied molecules was evaluated based on their ChemPLP scores (Table 1). Several molecules exhibited a higher binding affinity against TMPRSS2 than camostat (Table s1). The top 20 molecules were subjected to ADMET screening, and it was found that 7 molecules demonstrated outstanding ADMET properties (Table s2). The binding affinity of these seven compounds is shown in Table 1. The selected compounds displayed ChemPLP fitness scores of 73.12, 73.36, 71.71, 73.95, 76.99, 75.45, and 73. 50, which is considerably higher than the camostat (58.08) (Table 1).

Table 1

The binding details of camostat and other compounds against TMPRSS2.

Compounds	ChemPLP	Residues
		Hydrogen bonds	Alkyl/Amide-Pi Staked	Pi-Pi Staked	Attractive Charge	Van der Waals
Camostat	−58.08	C297, E299, K300, D435, S436, S441, W461		H296		V298, P301, L302, K342, C437, Q438, D440, T459, S460, G462, G464, A466, P471, G472
Compound 1 (20155148)	73.12	K342, D435, S436, S441, G464	C437	H296, W461	D435	Q438, G462
Compound 2 (53682039)	73.36	K342, D435, S436, S441, G464	K449, S436		K342, D435, W461	S339, K340, T341, L419, C437, Q438, G439, D440, T459, S460, S463, C465, A466, R470, P471, G472, Y474
Compound 3 (53793692)	71.71	K340, K342, D435, S436, S441, G464	K342, C437	H296, W461	D435	E299, Y337, S339, T341, L419, M424, C437, G439, D440, T459, S463, C465, A466, P471, G472
Compound 4 (53964549)	73.95	K342, D435, S436, S441, G464	K342, C437	W461	D435	H296, E299, S339, K340, T341, L419, Q438, G439, T459, G462, C465, A466, G472
Compound 5 (134379672)	76.99	C297, E299, D435, S436, S441, W461, G464, C465	C437		D435	V280, H296, V298, K300, P301, L302, K342, Y416, Q438, G439, T459, S460, G462, S463, A466, R470, P471, G472
Compound 6 (134379673)	75.45	E299, K300, D435, S436, Q438, S441, G464, C465	C437	W461		H296, C297, P301, L302, K342, E389, K390, Q438, T459, S460, G462, S463, C465, A466, G472
Compound 7 (139645059)	73.86	E389, D435, S436, G464	C437, C465	W461	D435	V280, H296, Q438, G439, T459, S460, G462, S463, A466, R470, P471, G472

The study also reveals the critical residues involved in accommodating camostat and other molecules within the binding site of TMPRSS2. D435, S436, S441, and G464 were engaged in hydrogen bond formation with the compounds (Table 1 and Fig. 3). Several other residues played a prominent role in accommodating all the compounds via Amide-Pi Staked, Pi-Alkyl, Alkyl, and Van der Waals interactions. Most notably, H296 and W461 were involved in pi-pi stacking interactions in most of the selected compounds. There were several other residues, namely K342, C437, Q438, D440, T459, S460, G462, G464, A466, P471, and G472, involved in van der Waals interaction with most of the compounds. Previous studies have well documented the contributory role of these active site residues (Rolta et al., 2021, Tateyama-Makino et al., 2021).

Fig. 3

The binding of Camostat and top selected compounds within the active site of TMPRSS2.

Molecular dynamics simulation studies

The top molecules in complex with TMPRSS2 were further subjected to molecular dynamics simulations to study the extent of the interactions of these compounds with TMPRSS2. Root mean square deviation (RMSD) analysis is one of the most significant approaches to investigating protein dynamics. We explored the TMPRSS2 protein backbone dynamics in complex with camostat and other selected inhibitors (Fig. 4a). The RMSDs of these inhibitors during the simulation period were determined (Fig. 4b). Ligands binding to their respective target proteins results in conformational structural changes within the resulting complex (Frimurer et al., 2003, Seeliger et al., 2010). The protein backbone RMSD was analyzed, and it was observed that protein backbone remains stable during the 100 ns time frame. Fig. 4a illustrates the backbone RMSDs of all the selected complexes. It was found that the backbone RMSD of all the complexes was stable throughout, with slight fluctuation observed (less than 2.5 Å). Overall, the backbone RMSD for the camostat and other compounds bound TMPRSS2 was stable throughout the simulation time period. The analysis of the ligand RMSD illustrates that the maximum fluctuation was observed in Compounds 1, 2, and 3 (Fig. 4b). Constant pose variation in the molecular conformation of the ligand was observed in the case of these compounds. For compound 4, these fluctuations were stabilized after the initial 45 ns, with minor deviations observed in the remaining period. For compounds 5 and 6, these fluctuations were observed to be stable with low structural deviations. For these compounds, the ligand RMSD was observed to be stable during the 100 ns MD run with deviations of less than 2 Å. Compound 7 was also observed to show some fluctuation level, but was stabilized after the initial 40 ns. Compared with the camostat, the protein backbone and ligand RMSD plots clearly show that the compounds 5, 6, and 7 bound complex of TMPRSS2 were stable, while compound 4 acquires stability after the initial deviations. The intermolecular hydrogen bond analysis (Fig. 4c) shows that all the compounds made an average of 2–4 hydrogen bonds throughout the 100 ns. For compounds 5,6, and 7, the highest hydrogen bonds were observed (upto 8). For compound 7, constant stable 6 hydrogen bonds were observed. As can be seen in the RMSF plot, a large degree of fluctuation was noticed in the amino acid residues involved in the loop region of TMPRSS2, with the maximum degree of fluctuation up to 4 Å (Fig. 4d). Careful analysis indicates that the residues involved in the binding of the ligand were observed to be having low RMSF values.

Fig. 4

Molecular dynamics results of the camostat and selected compounds bound complexes of the TMPRSS2. (a) The backbone RMSD of the TMPRSS2 in complex with the camostat and selected compounds (b) The ligand RMSD of the TMPRSS2 bound compounds and camostat (c) The intermolecular hydrogen bond formations of the TMPRSS2 in complex with the camostat and selected compounds (d) The RMSF of TMPRSS2 residues during the 100 ns. (e) The image indicates the SASA values of camostat and ligand molecules.

Further, the solvent-accessible surface area (SASA) was studied to evaluate the selected complexes' structural folding–unfolding dynamic under the solvent environment (Fig. 4e). For camostat, the average ligand SASA was 7.4 nm/NS2, while for other compounds, namely 1, 2, 3, 4, 5, 6, and 7, were 8.8 nm/NS2, 8.5 nm/NS2, 8.5 nm/NS2, 8.1 nm/NS2, 8.0 nm/NS2, 8.1 nm/NS2 and 7.4 nm/NS2, respectively. All the compounds were observed to have more solvent-accessible surface area than the camostat, suggesting that they should have more possibility of interaction with solvents. Moreover, constant flat SASA values for all compounds (except for compounds 1, 2, and 3) during the MD indicate the possibility of high stability.

Free energy calculation study

Accurately predicting the binding affinity of the compounds against the TMPRSS2 by measuring the binding free energy allows us to select the optimal compounds for TMPRSS2 inhibition.

The multi-layer screening of the compounds followed by precise prediction of their binding affinity against TMPRSS2 by measuring the binding free energy allows us to identify the best TMPRSS2 inhibitors. Table 2 displays the effect of each compound against TMPRSS2 by comparing the binding free energy values during the simulations. We evaluated the binding free energy of the complexes, estimated by using the MM-PBSA tool by AMBER, which is based on Poisson Boltzmann calculations performed using an internal PBSA solver in sander. Compounds 4, 5, 6 and 7 displayed the highest binding free energy of 50.793 +/-4.70 Kcal/mol, −64.121 +/-4.49 Kcal/mol, −60.052 +/-3.47 Kcal/mol and −51.318 +/-3.44 Kcal/mol respectively. These compounds demonstrated a much higher binding affinity than camostat (–33.352 +/- 3.41). Only Compound 1 showed a binding free energy value less than camostat. Detailed free energy values of the different compounds have been summarized in Table 2. The multi-layer screen of compounds increased the probability of getting an effective lead. The analysis of the selected compounds on different parameters suggests that these molecules are likely to be good hits in discovering TMPRSS2 inhibitors.

Table 2

Computed (MM-GBSA) binding free energies of the top selected compounds against TMPRSS2.

Compounds	ΔG bind (KCal/mol)	van der Waal energy (KCal/mol)	Electrostatic energy (KCal/mol)	Electrostatic solvation energy (KCal/mol)	Non-polar solvation energy (KCal/mol)	ΔG solv (KCal/mol)	ΔG gas (KCal/mol)
Camostat	–33.352 +/- 3.41	−43.649 +/- 3.19	−50.605 +/- 10.20	66.600 +/- 7.65	−5.699 +/- 0.27	60.902 +/- 7.62	−94.253 +/- 9.55
Compound 1	−30.582 +/-3.69	−25.573 +/-4.59	−152.802 +/-11.06	151.567 +/- 10.3	−3.773 +/-0.61	147.793+/-10.00	−178.375 +/- 12.32
Compound 2	−37.747 +/-5.15	−37.507 +/- 4.29	−148.037 +/-13.47	153.326 +/- 9.82	−5.529 +/-0.44	147.797 +/-9.76	−185.544 +/- 12.93
Compound 3	−41.194 +/-4.80	−34.012 +/-4.36	−151.259 +/-11.83	149.042 +/- 9.36	−4.966 +/-0.53	144.076+/-9.26	−185.271 +/- 12.28
Compound 4	−50.793 +/-4.70	−39.623 +/-4.65	−136.831 +/-12.14	131.811 +/-10.21	−6.150 +/-0.36	125.661 +/-10.14	−176.454 +/- 12.56
Compound 5	−64.121 +/-4.49	−44.930 +/-4.35	−194.950 +/-11.45	182.203 +/- 9.02	−6.444 +/-0.33	175.759 +/-8.91	−239.880 +/- 11.53
Compound 6	−60.052 +/-3.47	−46.597 +/-3.85	−174.761 +/-9.98	167.585 +/- 7.61	−6.279 +/-0.20	161.306 +/- 7.59	−221.359 +/- 9.01
Compound 7	−51.318 +/-3.44	−28.209 +/-3.62	−192.874 +/-12.68	173.990 +/- 9.98	−4.225 +/-0.30	169.765 +/- 9.88	−221.083 +/- 11.84

Conclusions

The critical role of TMPRSS2 in the viral entry within the host cell and replication makes it an attractive therapeutic target for inhibiting SARS-CoV-2 infection. We utilized state-of-the-art in silico approaches, including structure-based virtual screening, molecular docking, ADMET analysis, molecular dynamics simulation, and free energy calculations to identify TMPRSS2 inhibitors. Analysis of MD simulations along with the MM-PBSA calculations led to the selection of the four compounds (CID 53964549, 134379672, 134379673, and 139645059) as promising inhibitors for TMPRSS2. The selected compounds showed steady and stable binding to TMPRSS2 and better binding potential than camostat. Our findings clearly indicate the high TMPRSS2 inhibitory potential of the identified compounds and may serve as a therapeutic strategy to combat the SARS-CoV-2 infections efficiently.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.