Relaxed complex scheme and molecular dynamics simulation suggests small molecule inhibitor of human TMPRSS2 for combating COVID-19

Abstract As the coronavirus disease 19 (COVID-19) pandemic continues to pose a health and economic crisis worldwide, the quest for drugs and/or vaccines against the virus continues. The human transmembrane protease serine 2 (TMPRSS2) has attracted attention as a target for drug discovery, as inhibition of its catalytic reaction would result in the inactivation of the proteolytic cleavage of the SARS-CoV-2 S protein. As a result, the inactivation prevents viral cell entry to the host’s cell. In this work, we screened and identified two potent molecules that interact and inhibit the catalytic reaction by using computational approaches. Two docking screening experiments were performed utilizing the crystal structure and holo ensemble structure obtained from molecular dynamics in bound form. There is enhancement and sensitivity of docking results to the holo ensemble as compared to the crystal structure. Compound 1 demonstrated a similar inhibition value to nafamostat by interacting with catalytic triad residues His296 and Ser441, thereby disrupting the already established hydrogen bond interaction. The stability of the ligand–TMPRSS2 complexes was studied by molecular dynamics simulation, and the binding energy was re-scored by using molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) binding free energy. The obtained compounds may serve as an initial point toward the discovery of potent TMPRSS2 inhibitors upon further in vivo validation. Communicated by Ramaswamy H. Sarma


Introduction
The outbreak of coronavirus disease-19 (COVID-19) caused by a severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has resulted in global health and economic crisis. SARS-COV-2 is a single-stranded RNA virus belonging to the same family with SARS-CoV and the Middle East respiratory syndrome (MERS) . These viruses (SARS-COV-2, SARS-CoV and MERS-CoV) have been identified as highly pathogenic human coronaviruses, with an ability to infect human hosts . The spread rate of SARS-CoV-2 has been reported to exceed those of MERS-CoV and SARS-CoV . As of 8 June, 2021 approximately 3,764,513 deaths with over 174,822,669 cases were confirmed (https://www. worldometers.info/coronavirus/). Until now, there is no effective treatment, although some drugs have been reported and vaccines developed (Shadrack et al., 2021;Zhu et al., 2021). The long-term effects and the associated side effects are yet to be established, thus, scientists across the globe are working to develop effective drugs and vaccines to combat the disease. Among the approaches, drug repurposing has emerged as a promising way to identify drugs, whereby some drugs such as hesperidin and chloroquine have been repurposed for COVID-19 treatment and are in clinical trials (Haggag et al., 2020;Hu et al., 2021). Many other researchers have devoted efforts in searching for molecules that could block the SARS-COV-2 entry to human cells.
The SARS-CoV-2 entry into the human cell is mediated by a transmembrane spike (S) glycoprotein (Hu et al., 2021). Generally, the release of the viral RNA genome into host cells occurs after the S protein binds to the host receptor angiotensin-converting enzyme 2 (ACE2), and upon endocytosis the viral membrane fuses into the host via endosomal membrane (Hu et al., 2021). The binding and fusion of the viral S protein is mediated and activated by the human proteases (Hu et al., 2021). The cleavage sites S1/S2 and S 0 found at the boundary between S1 and S2 subunits are proteolytically cleaved by the proteases [transmembrane protease serine 2 (TMPRSS2), furin and cathepsin L and B] . Although both proteases have shown a potential in priming and activating the S protein, the TMPRSS2 plays an important role in activating the membrane functioning of the SARS-COV-2 S protein (Hu et al., 2021). Therefore, small molecules can be targeted to inhibit the proteolytic activity of TMPRSS2, hence blocking the viral cell membrane entrance to the host cell. Targeting the human TMPRSS2 has an advantage of overcoming the drug resistance, unlike targeting the viral proteins. Thus, in the effort to develop new anti-SARS-COV-2 drugs, the TMPRSS2 is considered as a promising target.
Several research groups have considered TMPRSS2 as the molecular target for identifying small molecules that can bind and inhibit the proteolytic process Hu et al., 2021;Huggins 2020;Hoffmann et al., 2020Hoffmann et al., , 2021Peiffer et al., 2021). For instance, camostat was first suggested by Hoffmann et al. (2020) as a potent inhibitor of TMPRSS2 which can block the entry of SARS-CoV-2. Later, the same author demonstrated that camostat mesylate and its active hydrolyzed metabolite 4-(4-guanidinobenzoyloxy)phenylacetic acid (GBPA) exerts antiviral activity by blocking the TMPRSS2 through covalent binding to the enzyme , the metabolite was further observed to be slowly metabolized, hence, exerting a prolonged antiviral activity of the administered single dose . Several research groups have considered other protease inhibitors as potential candidates for COVID-19, and some including nafamostat and the previously reported TMPRSS2 inhibitors camostat, successfully chanced into in clinical trials as inhibitors of TMPRSS2 . Although both nafamostat and camostat have demonstrated in vitro and in vivo inhibition to TMPRSS2 , the former has shown more potence nearly two-fold higher than the latter . Efforts to screen libraries containing small molecules for inhibiting TMPRSS2 are using in silico methods has been reported since the emergency of COVID-19. Hu and co-workers, using in silico approaches screened a library and identified potent small molecules (Hu et al., 2021). In addition, Huggins (2020) reported approved small molecules targeting inhibition of TMPRSS2 in silico (Huggins, 2020). These reports suggest that targeting TMPRSS2 has a great potential to fight against this pandemic.
Indeed, atomistic simulation methods have shown a great promise in drug discovery. Methods such as molecular docking are now widely used in screening large libraries of small molecules. However, docking algorithms are limited due to global protein flexibility. Incorporation of full flexible receptor in virtual screening has been suggested and showed improvement in docking score. Several groups have considered the relaxed complex scheme (RCS) docking based approach. In this approach, the apo or holo protein is subjected to molecular dynamics (MD) simulation and ensemble structures are extracted from the stable equilibrated structure and used for docking. In an effort to discover a more potent inhibitor of TMPRSS2, our group carried out a combination of in silico methods to screen potential candidates from NANPDB database using different screening methods based on the apo and holo ensembles. In this work, we found that holo ensemble structure improved docking score as compared to the crystal structure, the stability of the compound was established by using MD and binding free energy calculations based on molecular mechanics Poisson-Boltzmann surface area (MM-PBSA).

Screening of natural products
TMPRSS2 is a serine protease responsible for activating and mediating the S-protein of SARS-CoV-2 proteolytic cleavage for viral cell entry through ACE2. SARS-CoV-2 proteolytic cleavage processes result in the formation of S1/S2 and S 0 fragments which attach to the TMPRSS2 (Hu et al., 2021;Zhu et al., 2021). The catalytic triad site residues; His296, Asp345 and Ser441 and substrate binding sites residues Asp435, Ser460 and Gly462 form the human TMPRSS2 active site and are important for the SARS-COV-2 proteolytic cleavage and attachment Sgrignani & Cavalli 2021;Sonawane et al., 2021). Inactivation of TMPRSS2 enzymatic activities requires small molecules that bind strongly to these key residues. To identify novel small molecule inhibitors of TMPRSS2, we screened a database of African medicinal compounds against the catalytic and substrate binding site using the crystal structure (PDB ID: 7MEQ). Screened compounds showed the distribution binding free energy range of À9.6 to À3.2 kcal/mol (see Figure 1(a)). Nafamostat, an experimentally known inhibitor of TMPRSS2, was used as a standard and had a binding free energy of À8.1 kcal/mol. It is important to mention that, the control drug, nafamostat, which is a wide spectrum protease inhibitors was chosen based on its potent inhibition to TMPRSS2 with IC 50 of 55 nM as compared to other known inhibitors such as camostat whose IC 50 is 142 nM . Additionally, although both camostat and nafamostat have similar structures, nafamostat is slightly shorter and possess many aromatic rings making it more rigid unlike camostat, as a result, camostat been long and more flexible may have several off-target binding hotspots, hence making it not a suitable control drug in docking calculations Zhu et al., 2021). Thus, compounds with binding free energies À8.1 could be more potent and were selected for further studies. It is worth mentioning that the observed binding free energy distribution using the crystal structure suggested potential compounds that could strongly inhibit the catalytic proteolytic cleavage of TMPRSS2 to SARS-COV-2.
Although molecular docking is an indispensable tool for screening large libraries and predicting the binding modes, the docking algorithms have limitations on accommodating full (global) receptor flexibility, hence may result in generation of false positive lead molecules. To overcome the challenge of global protein flexibility in virtual screening, accommodating receptor flexibility is necessary in drug discovery. A relaxed complex scheme (RCS) employing apo or holo receptor structure has shown a great promise. Our previous experience (Shadrack et al., 2020) shows that holo ensemble improves the binding free energy as compared to the apo receptor ensembles. A similar observation was also reported by Xu and Lill (2011). To perform virtual screening with a flexible receptor, we performed a MD simulation of the TMPRSS2 bound with compound 1, which ranked higher with binding free energy of À9.6 kcal/mol compared to À8.1 kcal/mol of nafamostat (Figure 1(b)). In this work, ensemble structures used for RCS were obtained from the equilibrium trajectory with minimum free energy computed using RMSD as a reaction coordinate, and subjected to docking calculation ( Figure S1).
It should be noted that, during free energy calculations, the first 10 ns were regarded as equilibration period and hence excluded from the analysis. Interestingly, during the holo ensemble screening, we noted two observations; first, the binding free energy generally improved in the holo ensemble as compared to the crystal structure. Figure 1(a) shows a shift in normal distribution curve for the holo structure with an improved docking score to À10.4 kcal/mol. Another feature observed was that, some compounds which ranked higher in the crystal structure ranked least during the ensemble screening. This suggests that accommodating receptor flexibility is necessary in virtual screening which can ultimately help in eliminating false positive binders. For example, compound 1 ranked higher in the crystal structure screening with binding free energy of À9.6 kcal/mol. However, during the RSC, the binding free energy decreased to À7.6 kcal/mol. These observations suggest sensitivity of docking results to the crystal and holo ensemble structure.
To aid our understanding on the inhibition mechanism of nafamostat and compound 1 which ranked higher during the docking calculations, in Figure 1(c), we show the binding modes of the two compounds inside the TMPRSS2. It is interesting to observe that the guanidinium groups of nafamostat, and the sugar moiety of compound 1 (Figure 1(c)) all bound to the substrate site by interacting with Gly462 and Trp461 as shown in Figure 1(c). Although the two molecules possess different chemical scaffolds, they showed similar binding modes to TMPRSS2, by interacting with residues at the catalytic triad and substrate binding site by forming hydrogen bonds with His296, Gly462 and Trp461, a similar interaction observed to nafamostat was also reported by Zhu et al. (2021) and Hempel et al. (2021). Although the RCS suggested an improvement in docking score to some molecules, the thermodynamic stability and dynamics of the complexes were investigated by means of unbiased classical MD and end-point free energy based on MM-PBSA to ascertain the observed binding free energy.

MD simulation
MD simulation is an important tool in structure-based drug design, it provides dynamical and related stability of a complex. It should be noted that lower binding free energy obtained from docking calculations does not necessarily imply stability of the complex. In addition, end-point free energy methods are very useful in re-scoring poses obtained from docking calculations. In this work, nafamostat, a known TMPRSS2 inhibitor and compound 1 which ranked the least in holo and higher in crystal structure docking calculations (Figure 1(b)) were selected for MD calculations to establish the time-lapse stability.
The general stability of the complexes was attained by measuring the Ca root mean square deviation (RMSD) for TMPRSS2 bound with compound 1 and nafamostat. The timedependent RMSD value shows that TMPRSS2-compound 1 complex quickly equilibrated within 10 ns and remained stable with few noticeable fluctuations ( Figure S2). However, TMPRSS2-nafamostat complex equilibrated during the 40 ns and remained stable with few oscillations.
The complex formed by TMPRSS2-compound 1 showed an average Ca RMSD value of %0.23 nm and TMPRSS2-nafamostat showed a slight increase in Ca RMSD value of %0.3 nm, the observed values of Ca RMSD i.e.
3 nm suggests stability of the complexes (Figure 2(a)) (Wu et al., 2007). We further noted a longer equilibration phase in TMPRSS2-nafamostat complex, we hypothesized that the longer equilibration could be attributed by orientation and pose changes of nafamostat which occurred during the first 20 ns as confirmed by the ligand pose RMSD (Figure 2(b)). Interestingly, changes in binding orientation of nafamostat at the catalytic site have recently been noted . However, such orientations are not observed in TMPRSS2-compound 1 complex. As shown in Figure S3(a,b), the observed stability in compound 1 was enhanced by several hydrogen bonds formed between the residues mediated with interfacial water molecules at the vicinity of the catalytic triad. The number of hydrogen bonds formed between compound 1 are higher than that of nafamostat.
It is important to note that although complex RMSD is useful in assessing the general stability, it does not provide an in-depth understanding on ligand dynamics (pose changes) and stability. To understand whether compound 1 and nafamostat changed from their initial docking pose (association/dissociation), we computed the pose ligand RMSD which predicts drug efficacy in vitro or in vivo. The time-dependent pose RMSD suggests that compound 1 remained in its initial docking pose with few noticeable unbinding (dissociation/association) processes observed from 65 ns (Figure 2(b)). It is important to mention that even when compound 1 dissociated and diffused into solvents with the RMSD distance !8 nm, it eventually returned to the catalytic triad site (Figure 2(b)). Contrary to compound 1, nafamostat displayed a different binding behavior, (Figure  2(b)), in which, during the first 5 ns, it remained in its docking pose then moved to occupy a different pose ( Figure  2(d)). The changes in nafamostat binding site suggest an existence of transient binding pockets in TMPRSS2 which would result in an allosteric inhibition, however, this calls for further exploration. Although nafamostat is known to inhibit TMPRSS2 by binding to the catalytic and substrate site, the existence of a putative allosteric binding pocket in TMPRSS2 has been suggested recently (Sgrignani & Cavalli, 2021) however, as mentioned before, the allosteric mechanism of nafamostat still need further exploration. To gain a better understanding on the dynamical stability of compound 1 and nafamostat, the free energy as the function of ligand RMSD and ligand pose RMSD was calculated (Figure 2(c,d)). It is observed that, compound 1 although exhibited some dissociation, remained stable as shown in Figure 2(c). However, the binding of nafamostat exhibited two RMSD values with the most stable conformation at pose 2 ( Figure  2(d)). The observed change in binding pose is further confirmed by the distance between His296-nafamostat which formed hydrogen bond in pose 1 (Figure 2(e)). Our results are somehow different from the finding by Zhu et al. (2021), who reported the RMSD value for nafamostat to be 0.13 nm. However, in our work, we observed two conformational states with RMSD at 0.1 nm and 0.17 nm which still suggests stability. The RMSD value for compound 1 was found to be 0.22 nm, suggesting that all the drugs bind with stability to the TMPRSS2 region.
To understand the backbone stabilities and fluctuations of the two complexes, root mean square fluctuations (RMSF) was computed (Figure 2(f)). Two regions with residue 213-224 and 385-395 exhibited high fluctuations in the nafamostat complex. Such fluctuations are not pronounced in the complex formed with compound 1. This suggests that the binding of compound 1 to TMPRSS2 resulted in structural stabilities unlike nafamostat. The region with residue 460-466 also showed higher fluctuation when bound with nafamostat. Such fluctuations at this region are not observed in the apo and TMPRSS2-compound 1 complex. The observed fluctuations correspond to loops. Generally, the backbone of TMPRSS2-compound 1 exhibited lower fluctuations throughout the simulation as compared to the complex formed by TMPRSS2-nafamostat. The compactness of the complexes was assessed by measuring the radius of gyration (Rg) ( Figure S4). It is observed that all the protein showed preservation of the compactness with an average Rg value below 2.2 nm. Some fluctuations near 65 ns are observed for TMPRSS2-compound 1 which are explained by the dissociation of compound 1.
The observed dissociation and association of compound 1 could be attributed to; first, the thermal motion of the solvents, when compound 1 dissociated, interacted with solvent molecules until it was attracted by the residues at the TMPRSS2 surface. Second, the surface and the inner region of the catalytic region are made up of several negatively charged residues Glu/Asp, the positively charged groups of compound 1 were electrostatically attracted to the vicinity of TMPRSS2 binding substrate and eventually to the catalytic region where stabilizes by forming several hydrogen bonds. The observed changes in pose for nafamostat are explained by different conformational and positional changes in the catalytic site before adapting to the final pose. The hydrogen bonds formed at the region stabilized the ligands, hence causing critical disturbances in the already established hydrogen bonds between the catalytic residues His296 and Ser441; a similar observation for nafamostat was reported by Zhu et al. (2021).
The two residues His296 and Ser441 play an important role in the catalytic reaction for proteolytic cleavage of the S-Protein (Figure 3(a)). Any disturbances of these residues would result in the inactivation of the TMPRSS2 catalytic process, hence stopping the virus cell entry. As shown in Figure 3(a), the crystal structure shows that the catalytic triad residues are in the order of hydrogen bond distances. The binding of compound 1 has the role in disturbing the already established His296-Ser441 hydrogen bonding by sitting on the top of Ser441 and forming a p-stacking interaction, thereby inducing a steric hindrance resulting to an increased distance (Figure 3(b,d)). However, nafamostat binds differently; it interacts with the loop and b-sheet forming the His296 thereby inducing disturbances in the recognizing Ser441 (Figure 3(c,e-f)). Figure 3(b,c) shows that when compound 1 and nafamostat bound to TMPRSS2 resulted in an increased distances between His296 and Ser441 but with different magnitudes. As observed in the snapshot, the two molecules did not induce appreciable changes in Asp345-His296 distances ( Figure S5), suggesting that these compounds work by disturbing the His296-Ser441 interaction. The structural re-orientation of nafamostat resulted in the interaction with His279 and Val280 stabilized by the formation of hydrogen bonds; such interactions were previously noted by Zhu et al. (2021).

MM-PBSA free energy calculation
The main driving forces responsible for ligand-TMPRSS2 stable interaction were further analyzed by MM-PBSA free energy calculations. Table 1 shows that compound 1 binds more strongly than nafamostat; this is in agreement with the observed complex RMSD values (Figure 2(a)). As shown in Table 1, van der Waals forces are the main driving forces followed by non-polar energy terms. Electrostatic forces contributed least to the interaction while the polar forces resulted in an opposing effect on the ligand-TMPRSS2 interaction. The observed van der Waals forces contribution could be attributed to the group closeness between compound 1/nafamostat and TMPRSS2 in the binding site.
A closer observation shows that compound 1 had higher electrostatic contribution as compared to nafamostat. This could be attributed to the presence of many hydroxyl groups in compound 1 which interacts with the catalytic site of TMPRSS2. The free energy decomposition and per residue contribution for the two ligands are shown in Figure 4(a,c). The catalytic residue His296, and substrate residues Val280 and Trp461 highly contributed to the interaction with compound 1. It is observed that, the catalytic residue Ser441 contributed the least in the interaction as compared to His296, and, there is no any interaction contribution observed from Asp345 (Figure 4(b)). The distances measured between His296-Asp345 ( Figure S5) also support the per residue decomposition and contribution to the binding free energy analysis. It seems that, compound 1 would work by   Compound 1 À310.7 ± 23 À24.6 ± 10 158.5 ± 18 À76.2 ± 5 À253 ± 21 Nafamostat À160.2 ± 16 À3.8 ± 3 48.9 ± 13 À31.2 ± 4 À146.4 ± 14 strongly interacting with His296 and Ser441 thereby blocking their recognition. As mentioned earlier, nafamostat changed its binding pose, and still induced an appreciable distance between His296-Ser441. Per residue contribution and decomposition (Tables 2 and 3) suggests that, nafamostat would bind strongly with the residues Val280, Val275, Leu302 and His307, thereby indirectly inducing disturbances in the His296-Ser441 interaction as shown in Figures 3(d) and 4(d). Furthermore, as indicated in the two tables, the molecular mechanics energy term highly contributed to the stabilization of the complexes. Water is known to play an important role in protein-ligand interaction and stabilization. The role of water on the binding of compound 1 and nafamostat was also investigated to provide details on the observed binding stability and pose changes. During the first 5 ns of the simulation, it is observed that the head of nafamostat is directly bound at the interface of His296-Ser441. The two residues are observed to form a hydrogen bond which is surrounded by several water molecules; in which two vicinal water molecules are observed to mediate the interaction between nafamostat and protein by forming hydrogen bonds. The hydrogen bond between nafamostat with residues Lys342 and Gln438 are stabilized and mediated by water molecules. However, it is observed that water molecules formed stable interaction with the residues, thereby pushing nafamostat to another position, which is observed during the 40 ns ( Figure  5(a,b)). The head of nafamostat has slightly moved from the catalytic site and interacts strongly with Val275 and with His307 as confirmed with MM-PBSA free energy decomposition. However, such movement is not observed for compound 1. This is explained by the presence of several hydroxyl groups responsible for the formation of hydrogen bonds at the catalytic site mediated by water molecules (Figure S3(a,b)).

Discussion
TMPRSS2 is an important target toward developing an effective drug for COVID-19 treatment. Inhibitions of its proteolytic catalytic process will result in inactivation of the S protein from binding to the ACE2. Targeting inhibition by using small molecules has shown great promise in treatment of various viral diseases (Sarker et al., 2021;Wang et al., 2019) and has been suggested for COVID-19 in several reviews and research articles (AlAjmi et al., 2021;Al-Shari, 2020;Chikhale et al., 2021;Panoutsopoulos, 2020;Tian et al., 2021;Xiang et al., 2021). The use of natural products has shown a great potential in the pharmaceutical industry and a number of natural products approved for other diseases has increasingly reported (Newman & Cragg, 2016. This work has screened small molecules using computational approaches viz; molecular docking, MD and MM-PBSA binding free energy to identify small molecules from African medicinal plant databases that could be used to combat the disease. First, we demonstrate that the binding free energy of small molecules can be improved by holo ensemble structure generated from molecular simulation which overcomes the limitations of protein flexibility. Then, we have found that, the top ranked molecule in virtual screening also demonstrated an appreciable stability which is comparable to nafamostat, a known TMPRSS2 inhibitor. Although docking study showed that both molecules can bind to the catalytic triad site, however, the molecules exhibited some fluctuations. The RMSD and RMSF suggested that both compound 1 and nafamostat, an experimentally known TMPRSS2 inhibitor had stability with RMSD and RSMF value 0.3 nm, however, some RMSF oscillations above 0.3 nm were observed in the loop regions, which is expected. Vicinal water molecules in the catalytic triad site also played an important role in mediating the interaction between protein and the molecules in a different manner. Having several hydroxyl groups, compound 1 strongly interacted with catalytic site residues mediated with water molecules. Binding free energy based on MM-PBSA further suggested strong affinity of compound 1 with binding free energy of À253 kJ/mol over nafamostat with À146 kJ/ mol in the TMPRSS2 pocket, however, the difference in binding free energy is within the error. van der Waal forces were found to play big role in the stabilization of the complex. We have also observed an existence of a transient-like pocket in TMPRSS2 as shown by the binding of nafamostat, such pocket is also suggested by a recent work by Sgrignani and Cavalli (2021), however, further investigations are still needed to validate this observation. The methodology applied here suggests that compound 1 could strongly bind and inhibit the catalytic process of TMPRSS2 upon further in vitro and in vivo validation.

Relaxed complex scheme-based docking screening
Molecular docking using crystal structure was performed by using Autodock Vina (Trott & Olson, 2010). The structure was obtained from Northern and South African natural compounds (NANPDB) (https://african-compounds.org/nanpdb/) in their SMILE file format, energy minimized and finally converted to .pdbqt file format using OpenBabel tool (O'Boyle et al., 2011). The protein structure was obtained from the RCSB Protein Data Bank (PDB) with the following PDB ID: 7MEQ. Missing hydrogen at pH 7.4 and gasteiger charges Table 2 The per residue decomposition contribution to binding free energy (DG (kJ/mol)) for compound 1. À10.5 ± 0.16 6.8 ± 0.09 À1.1 ± 0.02 À4.8 ± 0.17 Table 3 The per residue decomposition contribution to binding free energy (DG (kJ/mol)) for nafamostat. were added to the protein structure and then converted to .pdbqt format. Docking screening experiments against TMPRSS2 were performed using our customized shell script. A restricted docking calculation targeting the catalytic and substrate region was performed. Achieving full protein flexibility is a major problem in docking calculations. Some docking tools are able to allow side chain residue flexibility. In this study, relaxed complex scheme based docking screening was performed as described in our previous works (Shadrack et al., 2020). Briefly, the complex with compound 1 which ranked higher in crystal structure docking was selected and subjected to MD simulation of 100 ns. RMSD was used to assess the stability of the complex, the system RMSD attained equilibrium state during the first 10 ns. To obtain the configurations for docking calculation, we computed the free energy (F) as the function of the RMSD as shown in Equation (1).

Residues
where, k B is the Boltzmann constant, T is the temperature and P is the probability distribution of the RMSD. Structures with RMSD value corresponding to minimum free energy were sampled for docking calculation. As for the crystal structure, ensemble structures were prepared in a similar way and docking screening calculations were performed using AutoDock Vina (Trott & Olson, 2010). It is important to highlight that our docking calculations were validated using a set of experimental data reported in our recent developed protocol (Shadrack et al., 2020) as well as comparison to RMSD values of the redocked molecule which gave an acceptable value of 2.5 Å.

MD simulation
Classical MD simulation of the complexes was performed as described in our previous works (Shadrack et al., 2020(Shadrack et al., , 2021. Here, we describe the methodology briefly, nafamostat and compound 1 were selected for classical MD simulation, the structures obtained from docking calculations were used as the starting coordinates for classical MD simulation. MD for holo proteins was performed using Gromacs 2018 employing Amber03 force field. The systems were energy minimized using the steepest descent algorithm, followed by equilibration at NVT (constant Number of particles, Volume and Temperature) ensemble for 1000 ps and NPT (constant Pressure) ensemble for 500 ps, with position restraints. The equilibrated system was then subjected to MD production run for 100 ns. During the equilibration stage, both temperature and pressure were maintained using Berendsen method (Nos e, 1986). For the production stage, Parinello-Rahman and v-rescale were used for pressure and temperature coupling at 300 K, respectively (Bussi et al., 2007;Parrinello & Rahman, 1981). Periodic boundary conditions (PBCs) were applied in all directions. Particle Mesh Ewald (PME) was used to treat long-range electrostatic interactions with a cutoff distance on 11 Å for both electrostatic and van der Waals interactions, while covalent bonds were constrained using the LINCS algorithm (Cheatham et al., 1995;Hess et al., 1997), a time step of 2 fs was used for all calculations. The MD configuration was used for analysis, gromacs built-in tools and our inhouse scripts were used for analysis.

MM-PBSA binding free energy calculations
Recently, it has been suggested that binding free energy and per-residue energy decomposition are very useful in establishing the underlying inhibition mechanism and biological activities target compound (Gupta et al., 2018). Herein, free energy calculations and the per-residue decomposition were performed to gain further understanding on the inhibition mechanisms of the studied molecules. MM-PBSA was calculated using the g mmpbsa tool (Kumari et al., 2014). The binding free energy (DG bind ) for protein-ligand interaction is expressed as where, G complex represents the total free energy of the protein-ligand complex while G protein and G ligand are total free energies of the individual protein and ligand in solvent, respectively. The free energy for each individual entity was computed as where, x stands for the protein or ligand or protein-ligand complex, T and S denote the temperature and entropy, respectively, and TDS refers to the entropic contribution to the free energy in a vacuum. The average molecular mechanics (MM) energy in vacuum is denoted as hE MM i, this term includes bonded and non-bonded interactions, hE MM i is calculated based on MM force-field parameters as follows where E bonded is bonded interactions consisting of bond, angle, dihedral and improper interactions. The nonbonded interactions (E nonbonded ) include both electrostatic (E elec ) and van der Waals (E vdW ) interactions and are modeled using a Coulomb and Lennard-Jones (LJ) potential function, respectively. The free energy of solvation hG sol i includes G polar and G nonÀpolar and can be calculated as where c is a coefficient related to surface tension, SASA stands for solvent accessible surface area and b is the fitting parameter, TDS is the entropic contributions to free energy. In g mmpbsa , an entropic contribution is not considered (Kumari et al., 2014). The binding free energy was calculated using a single trajectory, where a total of 200 snapshots were evenly extracted at a predetermined time from the first 40 ns and last 40 ns of the trajectory. The solvent and solute dielectric constants were 80 and 2, respectively, and c was 0.0227, the PB equation was solved by using the linear PBsolver.