Molecular modelling identification of phytocompounds from selected African botanicals as promising therapeutics against druggable human host cell targets of SARS-CoV-2.

The coronavirus disease 2019 (COVID-19), caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), is highly pathogenic and transmissible. It is mediated by the binding of viral spike proteins to human cells via entry and replication processes involving human angiotensin converting enzyme-2 (hACE2), transmembrane serine protease (TMPRSS2) and cathepsin L (Cath L). The identification of novel therapeutics that can modulate viral entry or replication has been of research interest and would be germane in managing COVID-19 subjects. This study investigated the structure-activity relationship inhibitory potential of 99 phytocompounds from selected African botanicals with proven therapeutic benefits against respiratory diseases focusing on SARS-CoV-2's human cell proteins (hACE2, TMPRSS2, and Cathepsin L) as druggable targets using computational methods. Evaluation of the binding energies of the phytocompounds showed that two compounds, Abrusoside A (-63.393 kcal/mol) and Kaempferol-3-O-rutinoside (-58.939 kcal/mol) had stronger affinity for the exopeptidase site of hACE2 compared to the reference drug, MLN-4760 (-54.545 kcal/mol). The study further revealed that Verbascoside (-63.338 kcal/mol), Abrectorin (-37.880 kcal/mol), and Friedelin (-36.989 kcal/mol) are potential inhibitors of TMPRSS2 compared to Nafamostat (-36.186 kcal/mol), while Hemiphloin (-41.425 kcal/mol), Quercetin-3-O-rutinoside (-37.257 kcal/mol), and Myricetin-3-O-galactoside (-36.342 kcal/mol) are potential inhibitors of Cathepsin L relative to Bafilomycin A1 (-38.180 kcal/mol). The structural analysis suggests that these compounds do not compromise the structural integrity of the proteins, but rather stabilized and interacted well with the active site amino acid residues critical to inhibition of the respective proteins. Overall, the findings from this study are suggestive of the structural mechanism of inhibitory action of the identified leads against the proteins critical for SARS-CoV-2 to enter the human host cell. While the study has lent credence to the significant role the compounds could play in developing potent SARS-CoV-2 candidate drugs against COVID-19, further structural refinement, and modifications of the compounds for subsequent in vitro as well as preclinical and clinical evaluations are underway.

Introduction

The coronavirus disease 2019 (COVID-19), caused by the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), is highly pathogenic and transmissible. The first identified and reported case of COVID-19 was in Wuhan, China in December 2019. As of 14 January 2022, WHO has reported 315 345 967 confirmed cases of COVID-19 with 5 510 174 deaths [1]. The most common clinical features of COVID-19 are tiredness, pyrexia, cough, and dyspnea, while some may remain asymptomatic [2].

Although, some vaccines have been approved and are currently in use, the vaccination program itself is riddled with problems including, misinformation and disinformation about the effectiveness of the various vaccines that are available, reduced potency against the variants, and post-vaccination symptoms, such as thromboembolism and thrombocytopenia and other adverse effects [3]. Hence, the development of novel drug candidates that can stand the test of time in offering therapeutics that will be globally acceptable, affordable, and easily accessible, is imperative. Drug development from plants has been demonstrated against several debilitating diseases including viral infections [4,5] and recent reports have lent credence to exploration of plant secondary metabolites as leads against druggable targets in SARS-CoV-2 infection using computational approaches [[6], [7], [8]].

As part of the human host cell proteins, important druggable targets including the human angiotensin converting enzyme 2 (hACE2), transmembrane protease serine 2 (TMPRSS2), and cathepsin L (Cath L), are required to facilitate and aid the viral entry and replication in SARS-CoV-2 infection. Specifically, infection with SARS-CoV-2 occur when the interaction between the viral spike protein and the human host cell is facilitated through processes involving hACE2, TMPRSS2 and Cath L [9]. In fact, the inhibitory effect of the depletion of hACE2 and TMPRSS2 on SARS-CoV-2 replication in different cell lines has been reported [10]. Similarly, the expression of Cath L is up-regulated during chronic inflammation as observed with cytokine storm in COVID-19 infection and is also implicated in extracellular matrix degradation, a critical process for SARS-CoV-2 viral entry into the host cell [11]. Generally, the host cell proteins are genetically more stable than viral structural proteins, which are also typical druggable targets on the viral genome [12,13]. Hence, the development of novel phytotherapeutics that can modulate viral entry or replication through inhibition of hACE2, TMPRSS2, and Cath L would be of immense benefit to mankind as we grapple with the COVID-19 pandemic.

To identify Phytocompounds as potential therapeutic agents against SARS-CoV-2 and its associated target proteins, several in silico and in vitro studies have explored the potentials of African rich plants. Dwarka et al. [14], reported four compounds (Uzarin, Hypoxide, l-canavanine and arabic acid) from South African medicinal plants as potent therapeutics against COVID-19 in an in silico study [14]. Similarly, Amaranthin from Amaranthus tricolor commonly used in Nigeria, Kenya and Tanzania was reported to be a potential inhibitor of the main protease enzyme of SARS-CoV-2 [15]. Compounds such as licoleafol, methyl rosmarinate, A-terpineol, P-cymene, T-anethole, and thyhydromoquinone are some of the other therapeutics that have been identified from African medicinal plants against the main protease, hACEE-2 and other key druggable targets of SARS-CoV-2 [[15], [16], [17], [18], [19], [20], [21]]. As part of the ongoing efforts to develop anti-SARS-CoV-2 drugs, the current study adopted a computational structure-activity relationship approach exploring 99 plant secondary metabolites from eight African medicinal plants (Leonotis leonurus, Ocimum gratissimum, Macaranga barteri, Abrus precatorius, Artemisia afra, Carapa procera, Alepidea amatymbica, and Drosera madagascariensis) with proven therapeutic benefits against respiratory diseases and related infections [[22], [23], [24], [25], [26], [27], [28], [29], [30]], against the SARS-CoV-2 host cell targets.

Materials and methods

Molecular docking and simulation

Protein acquisition, preparation, and docking

The protein structures of the hACE2 (PDB ID: 1R4L), TMPRSS2 (PDB ID: 5CE1) and Cath L (PDB ID: 5MQY) were obtained from the RSCB Protein Data Bank (https://www.rcsb.org). The preparation of all the protein structures were carried out on the UCSF Chimera version 1.14 by eliminating water molecules, nonstandard naming, and protein residue connectivity [31]. Prior to molecular docking, protein structures with missing atoms in their sidechains and protein backbone are rectified by adding such missing atoms. The reference/standard drugs as well as the phytocompounds used in this study were accessed and downloaded from PubChem [32]. Using Avogadro software [33], the 3-D structures of the selected 99 phytochemicals (Supplementary Table S1) and reference drugs [bafilomycin A1 (BFA), nafamostat (NFM), and MLN 4670 (MLN)] were all prepared in readiness for docking. While BFA was used as a known Cath L inhibitor, NFM and MLN were adopted as potent inhibitors of TMPRSS2 and hACE2, respectively [[34], [35], [36]].

For the molecular docking, the Autodock package on Chimera version 1.14 was used with default parameters [8]. Briefly, to the ligands were added Gasteiger charges and all the hydrogen atoms that are non-polar were joined to carbon atoms. Docking of all prepared secondary metabolites (ligands) into the binding site of the respective proteins (hACE2, TMPRSS2, and Cath L) by defining the grid box with a spacing of 1 Å each and size (26 × 24 × 25), (60 × 76 × 68) and (50 × 40 × 46) pointing in x, y and z directions, respectively. The resulting complex with the best geometric pose in each case was thereafter subjected to molecular dynamics (MD) simulation.

Molecular dynamics simulations, post-dynamic analyses and binding free energy calculation

The MD simulations were done as earlier reported [8] with the AMBER 18 suite Leap module [37] where hACE2 protein was numbered from residues 19 to 615, while Cath L, and TMPRSS2 proteins were numbered from residues 1 to 220 and 1 to 370, respectively. ANTECHAMBER was used to generate atomic partial charges for the ligands by utilizing the Restrained Electrostatic Potential (RESP) and the General Amber Force Field (GAFF) procedures. The Leap module of AMBER 18 allowed for the addition of hydrogen atoms, Na+ and Cl− counter ions for the three proteins, respectively, to neutralize all systems. Thereafter, the systems were anchored within an orthorhombic box of TIP3P water molecules within 10 Å [38], allowing a stepwise heating 0–300 K (50 ps) and a pressure of 1 bar [39]. The total time for the MD simulations conducted was 100 ns. In each simulation, the SHAKE algorithm was employed to constrict hydrogen atoms' bonds [40].

For the post-dynamic analysis of root mean square fluctuation (RMSF), root mean square deviation (RMSD), and radius of gyration (RoG), the CPPTRAJ module was adopted [41], and the resulting plots were generated in Origin V1.4 [42]. For the binding affinity, the binding free energy, Molecular Mechanics/GB Surface Area method (MM/GBSA) was estimated over 1 × 105 snapshots drawn from the 100 ns trajectory [43].

Pharmacokinetic properties

An assessment of the chemistry and drug likeness of the lead phytocompounds was done using the SwissADME online software [44].

Results and discussion

Docking scores and thermodynamic binding free energy

The results of the molecular docking of the 99 phytocompounds investigated against hACE2, TMPRSS2, and Cath L in this study are presented in Supplementary Table S1. Molecular docking enables the assessment of the geometric fitness and affinity of a molecule upon its binding at the active site of the receptor and the higher the negative score, the better the pose and interaction of the compound with the protein [45,46]. For hACE2 (Table S1), 16 phytocompounds had docking scores between −10.1 and −9.1 kcal/mol and better affinity for the protein compared to the reference standard, MLN-4760 with a docking score of −7.3 kcal/mol (Table S1). For TMPRSS2 (Table S1), 11 phytocompounds had docking scores between −9.2 and −7.9 kcal/mol compared to −8.1 kcal/mol for the reference standard, NFM while for Cath L (Table S1), 27 phytocompounds had docking scores between −8.7 and −7.0 kcal/mol compared to BFA (−8.1 kcal/mol). Since molecular docking assesses only the pose and affinity of a molecule in the receptor active site or protein, the most promising of these compounds against each target were further taken through MD simulation and the results of those with better or close binding free energy post-MD simulation are presented in Supplementary Tables S2, S3, and S4 while the results of the most promising ones judging by the highest negative values are shown in Table 1, Table 2, Table 3 . Specifically, for hACE2, Abrusoside A (ABA) and Kaempferol-3-O-rutinoside (KOR) had the best binding free energy values of −63.393 kcal/mol and −58.939 kcal/mol, respectively, relative to MLN (−54.545 kcal/mol) (Tables 1 and S2), while Verbascoside (VBS), Abrectorin (ABC), and Friedelin (FDL) had the best binding free energy values of −63.338 kcal/mol, −37.880 kcal/mol, and −36.989 kcal/mol, respectively against TMPRSS2 relative to −36.186 kcal/mol for NFM (Tables 2 and S3). While the observation on hACE2 is indicative of ABA's affinity for the protein and suggestive of its ability to be a better inhibitor of hACE2 than KOR, VBS could be suggested as a strong inhibitor of TMPRSS2 relative to both ABC and the reference standard NFM. Since TMPRSS2 is responsible for cleaving the viral spike glycoprotein, VBS could be a promising therapeutic drug candidate targeting the entry stage of SARS-CoV-2 replication cycle. It was also noteworthy that both ABC and FDL bound strongly with TMPRSS2 judging by their binding energy values that were comparable to that of NFM. On the other hand, Hemiphloin (HPN): -41.425 kcal/mol, Quercetin-3-O-rutinoside (QOR): -37.257 kcal/mol, and Myricetin-3-O-galactoside (MOG): -36.998 kcal/mol were identified as prominent compounds against Cath L (Tables 3 and S4). The results suggest that HPN is the best inhibitor against Cath L as it showed higher binding affinity relative to BFA and other compounds used in this study. The observations noted regarding the binding energy values of the tested compounds against hACE2, TMPRSS2 and Cath L are similar with previous reports [47,48], where compounds with the highest negative binding energy values were proposed as having the best affinity for SARS-CoV-2 host cell proteins and suggested as potential inhibitors of the respective druggable target.

Table 1

Thermodynamic binding free energy values for MLN, ABA, and KOR towards hACE2.

Complex	MLN 4670 (MLN)	ABA	KOR
ΔEvdW	−32.851 ± 4.761	−67.659 ± 3.783	−43.747 ± 6.103
ΔEelec	−942.921 ± 19.627	−23.862 ± 0.163	−99.020 ± 15.204
ΔGgas	−975.773 ± 20.535	−91.521 ± 7.533	−142.768 ± 12.274
ΔGsolv	921.228 ± 16.585	28.127 ± 5.659	83.828 ± 9.960
ΔGbind	−54.545 ± 7.029	−63.393 ± 4.757	−58.939 ± 5.538

ΔEelec electrostatic energy, ΔEvdW van der Waals energy, ΔGbind total binding free energy, ΔGsolv solvation free energy, and ΔGgas gas-phase free energy.

Table 2

Thermodynamic binding free energy values for NFM, VBS, ABC, and FDL towards TMPRSS2.

Complex	NFM	VBS	ABC	FDL
Δ EvdW	−34.007 ± 3.207	−56.875 ± 4.933	−33.932 ± 2.578	−36.624 ± 3.949
ΔEelec	−250.376 ± 20.327	−80.487 ± 16.014	−8.281 ± 0.068	−7.647 ± 0.116
ΔGgas	−284.383 ± 20.514	−137.363 ± 14.311	−42.213 ± 4.033	−47.272 ± 6.592
ΔGsolv	248.196 ± 19.035	74.025 ± 9.886	17.333 ± 2.702	15.235 ± 3.517
ΔGbind	−36.186 ± 4.572	−63.338 ± 7.493	−37.880 ± 2.596	−36.989 ± 4.036

ΔEelec electrostatic energy, ΔEvdW van der Waals energy, ΔGbind total binding free energy, ΔGsolv solvation free energy, and ΔGgas gas-phase free energy.

Table 3

Thermodynamic binding free energy values for BFA, QOR, HPN, and MOG towards Cath L.

Complex	BFA	QOR	HPN	MOG
Δ EvdW	−49.969 ± 5.756	−38.005 ± 3.345	−40.191 ± 3.216	−32.119 ± 5.521
ΔEelec	−13.595 ± 5.542	−38.535 ± 4.666	−26.703 ± 3.067	−42.1420 ± 0.564
ΔGgas	−63.564 ± 9.447	−76.540 ± 14.689	−71.896 ± 7.478	−74.261 ± 9.889
ΔGsolv	25.384 ± 4.956	44.283 ± 2.141	30.470 ± 5.052	42.263 ± 5.297
ΔGbind	−38.180 ± 5.985	−37.257 ± 7.588	−41.425 ± 3.869	−36.998 ± 3.097

ΔEelec electrostatic energy, ΔEvdW van der Waals energy, ΔGbind total binding free energy, ΔGsolv solvation free energy, and ΔGgas gas-phase free energy.

Interaction plots of the phytocompounds with the host cell proteins

The inhibitory characteristics of the promising or lead compounds, judging by their binding affinities and those of the respective reference standards against the evaluated targets as a function of interactions with the amino acid residues at the active site of each protein are presented in Fig. 1, Fig. 2, Fig. 3 . Different degrees of bond interactions such as van der Waals (vdW) overlaps, halogen, hydrogen bonds, alkyl, π-alkyl, π-π stacked interaction, and π-π T-shaped were observed. Specifically, Fig. 1 revealed that ABA bound strongly to hACE2 with 4 conventional hydrogen bonds (with Gly187, Glu190, Asn192 and Lys544) and 13 van der Waals’ interactions (with Leu55, Glu80, Gln84, Tyr178, Tyr184, Val191, Val194, Leu373, Leu374, Asn376, Tyr492, Ser545 and Gln546) while 6 conventional hydrogen bonds (with Ala330, Glu357, Asp364, Glu380, Asn490 and Arg496) and 14 van der Waals interactions (with Hie327, Pro328, Trp331, Asp332, Ile361, Tyr367, Asn376, Gly377, Hie387, Glu384, Phe486, Hie487, Asp491 and Tyr497) were observed with KOR. MLN on the other hand had a total of 13 hydrogen and van der Waals interactions with hACE2. In addition, ABA had four alkyl interactions with Phe22, Trp51, Phe372 and Ala81 compared to two strong π-π stacked interactions with Tyr492 and Tyr184 for KOR and two alkyl interactions with Phe372 and Tyr492 as well as one π-π stacked interaction with Phe22 for MLN (Fig. 1). The interactions observed with both ABA and KOR towards hACE2 could be a probable justification for their higher binding affinities relative to MLN. Hydrogen bonds as well as van der Waals and other non-covalent interactions such as alkyl interactions are known to significantly add to the binding energy values of ligands after binding to a receptor [49,50] and based on this observation, ABA and KOR could be identified to have had good interactions with the protein in a manner that enhanced affinity suggestive of their potential inhibitory effect on hACE2.

Fig. 1

2D interaction plots of MLN, Abrusoside A (ABA) and Kaempferol-3-O-rutinoside (KOR) with the active site amino acid residues of hACE2.

Fig. 2

2D interaction plots of Nafamostat (NFM), Verbascoside (VBS), Friedelin (FDL), and Abrectorin (ABC) with the active site amino acid residues of TMPRSS2.

Fig. 3

2D plots of Bafilomycin A1 (BFA), Myricetin-3-O-galactoside (MOG), Quercetin-3-O-rutinoside (QOR), and Hemiphloin (HPN) with active site amino acid residues of cathepsin L.

Fig. 2 shows the 2D interaction plots of VBS, ABC and FDL with the active site amino acid residues of TMPRSS2 in comparison with NFM. More interactions including 24 hydrogen and van der Waals forces and 3 strong amide-stacked bonds with His158, Val330, and Pro161 were observed between VBS and the active site amino acid residues of TMPRSS2, justifying its higher negative binding energy value and better affinity for the protein relative to others (Fig. 2). Although, strong amide-stacked and alkyl interactions involving either carboxylic group or benzyl group were observed in all the compounds including NFM, however, compared to VBS, NFM, FDL and ABC had lower number of interactions with the amino acid at the binding site of the protein (Fig. 2), and this also correlates with their respective lower binding energy values.

The structure of Cath L showed the fold of the papain-like enzyme composed of two domains (left (L-) and right (R-)), and between the two domains, is a V-shaped active site cleft on which the L- and R-domain catalytic residues C25 and H163 are positioned [51], and the ligands binds at this active site. For Cath L, the ligand-interaction plot revealed that the reference standard, BFA had a total of 11 interactions consisting of 4 π-alkyl interactions with Met70, Leu69, Ala210, and Ala135, 2 conventional hydrogen bonds, as well as 5 van der Waals forces (Fig. 3). Unlike the BFA, the study compounds (MOG, QOR and HPN) had higher number of hydrogen bonds in addition to van der Waals forces. Furthermore, BFA did not show any π-alkyl interaction with the protein as observed with the three study compounds. Thus, of the three studied ligands (MOG, QOR and HPN), HPN by virtue of the increased number of electrostatic interactions exhibited, will have better binding affinities than the reference drug BFA. This finding is affirmed by the fact that the magnitude of the binding is a measure of how strong the interactions are between the ligand and the protein, and this is usually increased by interactions such as van der Waals forces, electrostatic interactions, and hydrogen bond between two molecules [49].

Dynamic stability, compactness, and flexibility of hACE2, TMPRSS2, and cathepsin L bound and unbound complexes

The data obtained with respect to structural analyses of hACE2, TMPRSS2, and Cath L complexes are presented in Fig. 4, Fig. 5, Fig. 6 , respectively alongside the calculated average values of RMSD, RoG, SASA, and RMSF (Table 4 – 6). The RMSD is a parameter that assesses how stable a complex is. The lower the average RMSD value, the more stable the complex [52,53]. The RMSD plot for the hACE2 complexes (Fig. 4a) showed that the bound and unbound (apo-enzyme) complexes converged at approximately 5 ns, and thereafter both the apo-enzyme and the respective ligand-bound complexes displayed a favorable stability throughout the simulation with overall average values of 1.771 Å, 1.756 Å, and 1.873 Å), for MLN, ABA and KOR, respectively (Table 4). Although, binding of KOR revealed an insignificant marginally higher RMSD value when compared to ABA and MLN, the binding of the three compounds do not alter the overall structural stability of hACE2. In the case of TMPRSS2, the average RMSD values for VBS (2.487 Å), FDL (2.238 Å), and ABC (2.409 Å) complexes are relatively higher than that of NFM (1.809 Å) but lower than the value for the apo-enzyme (2.589 Å) while the reference drug NFM had an average value of 1.809 Å (Fig. 5a, Table 5 ). This observation on the RMSD values for TMPRSS2 indicates that the binding of NFM, VBS and FDL induced more structural stability on TMPRSS2, which does not only agree with the results of the binding energy values and suggestive of favorable interactions between the protein and both VBS and FDL but also identifying them as promising prospect against TMPRSS2 inhibition. For the Cath L complexes, the binding of the reference standard, BFA (2.287 Å) and MOG (1.723 Å) raised the average RMSD value of the complexes, while compounds HPN (1.012 Å) and QOR (1.382 Å) lowered the RMSD values when compared with the apo-enzyme (1.572 Å) (Table 6 ). The RMSD plots revealed that both BFA and MOG complexes at approximately 55 ns induced unstable conformational changes on the protein structure as evidenced by the relatively high RMSD values (Fig. 6a). Although, the result of the binding energy suggests that the phytochemicals might be effective inhibitors of Cath L, however, the inhibition mechanism displayed by HPN and QOR may differ to that of BFA and MOG. Generally, the observations on the RMSD values and patterns for the druggable targets in this study are good indications of the test compounds as prospective drug candidates judging by their values which were <3.5 Å acceptable limit and agrees with a previous study [54], where plant secondary metabolites induced significant structural stability at the druggable sites of SARS-CoV-2.

Fig. 4

Comparative plots of C-α atoms of hACE2 with MLN, ABA and KOR displayed as (a) RMSD, (b) RoG, (c) SASA, and (d) RMSF, post-100 ns MD simulation.

Fig. 5

Comparative plots of C-α atoms of TMPRSS2 with NFM, VBS, ABC, and FDL displayed as (a) RMSD, (b) RoG, (c) SASA, and (d) RMSF, post-100 ns MD simulation.

Fig. 6

Comparative plots of C-α atoms of Cathepsin L with BFA, QOR, MOG, and HPN, displayed as (a) RMSD, (b) RoG, (c) SASA, and (d) RMSF, post-100 ns, MD simulation.

Table 4

Calculated average values of RMSD, RoG, SASA, and RMSF of hACE2 complexes.

Complex	RMSD (Å)	RoG (Å)	SASA (Å2)	RMSF (Å)
hACE2	1.822 ± 0.101	24.012 ± 1.733	25065.454 ± 800.34	1.119 ± 0.021
hACE2 + MLN	1.771 ± 0.081	24.116 ± 2.034	25475.381 ± 723.43	1.199 ± 0.143
hACE2 + ABA	1.756 ± 0.123	24.076 ± 2.120	25040.439 ± 523.33	1.180 ± 0.474
hACE2 + KOR	1.873 ± 0.072	24.161 ± 3.012	25722.981 ± 635.24	1.219 ± 0.352

Table 5

Calculated average values of RMSD, RoG, SASA, and RMSF of TMPRSS2 complexes.

Complex	RMSD (Å)	RoG (Å)	SASA (Å2)	RMSF (Å)
TMPRSS2	2.589 ± 0.212	21.703 ± 2.423	17007.673 ± 903.423	1.382 ± 0.324
TMPRSS2 + NFM	1.809 ± 0.039	21.617 ± 3.024	16556.672 ± 950.245	1.194 ± 0.154
TMPRSS2 + VBS	2.487 ± 0.322	21.610 ± 1.948	16460.933 ± 746.324	1.347 ± 0.424
TMPRSS2 + FDL	2.238 ± 0.313	21.661 ± 3.242	16807.211 ± 932.394	1.295 ± 0.322
TMPRSS2 + ABC	2.409 ± 0.301	21.738 ± 2.452	17032.312 ± 732.421	1.294 ± 0.324

Table 6

Calculated average values of RMSD, RoG, SASA, and RMSF of Cathepsin L complexes.

Complex	RMSD (Å)	RoG (Å)	SASA (Å2)	RMSF (Å)
Cath L	1.572 ± 0.422	16.791 ± 2.341	9435.845 ± 232.494	1.093 ± 0.322
Cath L + BFA	2.287 ± 0.221	17.019 ± 3.494	9857.129 ± 195.353	1.305 ± 0.234
Cath L + QOR	1.382 ± 0.123	16.816 ± 2.344	9490.176 ± 302.355	1.036 ± 0.156
Cath L + MOG	1.723 ± 0.164	16.634 ± 1.344	9143.020 ± 235.462	1.126 ± 0.253
Cath L + HPN	1.012 ± 0.423	16.678 ± 3.223	8907.722 ± 321.452	0.820 ± 0.012

The RoG is a parameter used to evaluate the structural compactness of proteins when they bind with molecules [55]. A lower and higher RoG values indicate more stable and unstable systems, respectively [56]. The hACE-2 RoG plot results correlates with its RMSD plot revealing the that binding of the molecules do not alter the hACE-2 structural stability. Fig. 4b as well as Table 5 showed a relatively close average RoG values of 24.012 Å, 24.116 Å, 24.076 Å and 24.161 Å for the apo-enzyme, MLN, ABA and KOR, respectively. For the TMPRSS-2 complexes, lower RoG values and more stable complexes were observed for the apo-enzyme as well as ligand bound complexes (Fig. 5b). The average RoG values of 21.703 Å, 21.617 Å, 21.610 Å, 21.661 Å and 21.738 Å were recorded for the apo-enzyme, NFM-, VBS-, FDL- and ABC- complexes, respectively. These results validate the structural stability observed with the RMSD plot. As observed in the RMSD plot for BFA complex with Cath L, the binding of BFA increased its average RoG values (17.019 Å) compared to the apo-enzyme (16.791 Å) and other ligands with average RoG values of 16.816 Å (QOR), 16.634 Å (MOG) and 16.678 Å (HPN) (Fig. 6b, Table 6). Results revealed that the apo-enzyme together with QOR, MOG, and HPN had low and relatedly close values, suggesting an increased structural compactness and stability, while the increased average RoG value observed with BFA-complex indicates a decrease in Cath L structural compactness.

The SASA plot measures the exposure of the protein structure to solvent environment, and the lower the SASA value, the more exposed the hydrophobic amino acid residues of the proteins are, and the systems stability also increases [57]. From Fig. 4c and Table 4, the average SASA values of the hACE2 complexes were 25040.439 Å2, 25722.981 Å2 and 25475.381 Å2 for ABA, KOR and MLN, respectively. ABA has a lower value than the standard drug (MLN) (25475.381 Å2) while KOR has a higher value than the apo-enzyme as well as the MLN-hACE2 complex. These values are relatively close to the value for the apo-enzyme (25065.454 Å2). This result showed that, the binding of ABA and KOR to hACE2 compared favorably with that of MLN and does not alter the exposure of the buried hydrophobic residues of hACE2, and ultimately did not adversely impact the systems’ stability. Unlike hACE2, the SASA values for the TMPRSS2 complexes showed that, binding of NFM (16556.672 Å2), VBS (16460.933 Å2), and FDL (16807.211 Å2) marginally lowered the SASA values compared to the apo-enzyme (17007 Å2) (Fig. 5c, Table 5), which is an indication of increased exposure of hydrophobic amino acid residues of TMPRSS2, suggestive of increased structural stability. However, the binding of ABC (17032.312 Å2) showed no effect on the exposure of the hydrophobic residues of TMPRSS2 to solvent environment and the overall observation regarding the SASA values of the TMPRSS2 complexes in this study agrees with those of RMSD and RoG following binding of VBS, FDL, and ABC and a further attestation that, the structural integrity, which is significant for inhibitory activity of TMPRSS2, was not compromised. It was observed that BFA displayed a similar trend in its SASA, RMSD and RoG plots, where the binding of BFA increased the average SASA value (9857.129 Å2) compared to the apo-enzyme (9435.845 Å2) (Fig. 6c, Table 6). This suggests that the interactions of BFA with Cath L reduced the exposure of hydrophobic residues of Cath L and its structural stability. The binding of HPN (8907.722 Å2), and MOG (9143.020 Å2) lowered the average SASA values, while the binding of QOR (9490.176 Å2) does not affect the SASA values compared to the apo-enzyme. These observations suggest that the binding of HPN and MOG increased the exposure of hydrophobic amino acid residues of Cath L leading to its increased stability. A recent report by Hassan [58] indicated MOG to be a promising inhibitor against SARS-CoV-2 spike glycoprotein (Sgp) in silico. This result adds support to the recommendation of MOG as a potential inhibitor of SARS-CoV-2 based on the findings of our study.

The RMSF is an assessment of how the amino acid residues of a receptor move or fluctuates as a result of a binding of a drug [50,59,60]. Increased RMSF value is an indication of increase in flexible movements of the alpha-carbon atoms [8]. In this study, the binding of KOR (1.219 Å) and MLN (1.199 Å) marginally increased the average RMSF value compared to the apo-enzyme (1.119 Å) for the hACE2 complexes (Fig. 4d, Table 4), indicating the two ligands caused an overall amino acid residues flexibility of hACE2. However, the binding of ABA lowered the average RMSF value (1.180 Å), and this could be indicative of restricted movement of the hACE2 active site amino acid residues. However, relatively high fluctuations were observed at residues 75–85, 305–310, 535–545. Similarly, a relatively close, low average RMSF values were observed with both the unbound and bound complexes in the TMPRSS2 systems (Fig. 5d, Table 5). Nevertheless, high residual fluctuations at residues 150–190, 250–275 recorded. The binding of NFM and the studied phytochemicals (VBS, FDL, and ABC) for TMPRSS lowered the RMSF values compared to the apo-enzyme (Table 5), indicating an overall flexible movement and stable complexes. This stable ligand-enzyme complexes correlates with the observations with RMSD, RoG, and SASA plots for TMPRSS2. For the Cath L systems, the binding of BFA (1.305 Å) and MOG (1.126 Å) increased the average RMSF values of the complexes as compared to the unbound enzyme (Cath L) (1.093 Å) (Fig. 6d, Table 6), suggesting their binding leads to a decrease or restricted movement of the amino acid residues, making the protein less flexible. In contrast, the binding of HPN and QOR resulted in the structural flexibility of the enzyme's amino acid residues as indicated by the lowered RMSF values of 0.082 Å and 1.036 Å, respectively (Fig. 6d, Table 6).

Pharmacokinetics

As shown in Table 7 , nine compounds including the standards (KOR, VBS, ABC, QOR, MOG, HPN, NFM, BFA, and MLN) are either moderately soluble or soluble in water. This will enhance the bioavailability of the nine compounds. Four compounds (ABA, FDL, ABC, and HPN) passed the drug-likeness test. Lipinski's rule of five is a tool used in assessing the drug-likeness compounds whereby a determination is made as to whether any compound that possesses certain pharmacological activities has oral drug properties that will allow it to be administered in humans [61]. On the other hand, VBS, KOR, QOR and MOG did not pass Lipinski's by violating three rules.

Table 7

Predicted ADME parameters for the most prominent compounds.

Compound	BioavailabilityScore	Lipophilicity (iLOGP)	BBBPermeability	GITAbsorption	WaterSolubility	Drug Likeness (Lipinski)a
ABA	0.55	7.44	No	Low	Poor	1 (Mwt)
KOR	0.17	−0.73	No	Low	Soluble	3 (Mwt,HBD, HBA)
VBS	0.17	−0.43	No	Low	Soluble	3 (Mwt, HBD, HBA)
FDL	0.55	7.44	No	Low	Poor	1 (LogP)
ABC	0.55	2.44	No	High	Moderately soluble	0
QOR	0.17	−1.12	No	Low	Soluble	3 (Mwt, HBD, HBA)
MOG	0.17	−0.96	No	Low	Soluble	2 (HBD, HBA)
HPN	0.55	−0.25	No	Low	Soluble	1 (HBD)
NFM	0.55	2.14	No	Low	Moderately soluble	0
BFA	0.55	4.10	No	Low	Soluble	1 (Mwt)
MLN	0.55	2.06	No	High	Soluble	0

Mwt = molecular weight; HBD = number of hydrogen bond donor groups; HBA = number of hydrogen bond acceptor group Parts; LogP = Partition Coefficient.

All the compounds and reference standards except for ABC and MLN were predicted to exhibit low gastrointestinal tract absorption (GIT). Nevertheless, low GIT absorption does not alter the therapeutic activity of drugs as many drugs (such as Nafamostat, Bafilomycin A1, and Lopinavir) with low GIT absorption are commercially available and therapeutically potent. The level and rate at which the active moiety or metabolite enters systemic circulation are measured by drug bioavailability [61]. Compared with the reference standards, four compounds, KOR, VBS, QOR, and MOG exhibited low bioavailability scores of 0.17, while the bioavailability scores predicted for the other compounds (ABA, FDL, ABC, HPN) and the three reference standards was 0.55.

Conclusion

The modulation of the entry steps and inhibition of these host cell druggable targets could help achieve the inhibition of the viral entry by preventing the interactions between the host cell and the viral proteins. In this study, two compounds (ABA and KOR) have strong affinity for the exopeptidase site of hACE2, interacted tenaciously with the essential amino acid residues (Glu384, Asn376, His383, Arg496, Thr329, Asp364) needed for the catalytic activity of hACE2. The study further revealed each of the three compounds [(VBS, FDL and ABC) and (QOR, MOG and HPN)] as potential inhibitors of TMPRSS2 and Cath L proteins, respectively, as they do not compromise the structural integrity of the proteins, but rather stabilized and established catalytic interactions with the vital amino residues needed for inhibition of the respective targets. In a nutshell, the results obtained in this study are suggestive of the structural mechanisms of inhibition of the identified leads against the proteins critical for SARS-CoV-2 to enter the human host cell and cause infection. Further refinement and development of the compounds for subsequent in vitro and other preclinical and clinical evaluations are underway.

Data availability

The data used to support the findings of this study are included within the article.

Author's contribution

SS conceptualized and supervised the study. UJO, IKA and SS generated and analyzed the data. SFO co-supervised the project, while UJO wrote the manuscript. All authors read and contributed to critical review of the manuscript for intellectual content and approved the submission for publication.

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.