Abstract
The coronavirus disease 2019 (COVID-19) is a pandemic that has severely posed substantial health challenges and claimed millions of lives. Though vaccines have been produced to stem the spread of this disease, the death rate remains high since drugs used for treatment have therapeutic challenges. Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the virus that causes the disease, has a slew of potential therapeutic targets. Among them is the furin protease, which has a cleavage site on the virus’s spike protein. The cleavage site facilitates the entry of the virus into human cells via cell–cell fusion. This critical involvement of furin in the disease pathogenicity has made it a viable therapeutic strategy against the virus. This study employs the consensus docking approach using HYBRID and AutoDock Vina to virtually screen a pre-filtered library of 3942 natural product compounds of African origin against the human furin protease (PDB: 4RYD). Twenty of these compounds were selected as hits after meeting molecular docking cut-off of − 7 kcal.mol−1, pose alignment inspection, and having favorable furin-ligand interactions. An area under the curve (AUC) value of 0.72 was computed from the receiver operator characteristic (ROC) curve, and Boltzmann-enhanced discrimination of the ROC curve (BEDROC) value of 0.65 showed that AutoDock Vina was a reasonable tool for selecting actives for this target. Seven of these hits were proposed as potential leads having had bonding interactions with catalytic triad residues Ser368, His194, and Asp153, and other essential residues in the active site with plausible binding free energies between − 189 and − 95 kJ/mol from the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations as well as favorable ADME/Tox properties. The molecules were also predicted as antiviral, anti-inflammatory, membrane permeability inhibitors, RNA synthesis inhibitors, cytoprotective, and hepatoprotective with probable activity (Pa) above 0.5 and probable inactivity values below 0.1. Some of them also have anti-influenza activity. Influenza virus has many similarities with SARS-CoV-2 in their mode of entry into human cells as both are facilitated by the furin protease. Pinobanksin 3-(E)-caffeate, one of the potential leads is a propolis compound. Propolis compounds have shown inhibitory effects against ACE2, TMPRSS2, and PAK1 signaling pathways of SARS-CoV-2 in previous studies. Likewise, quercitrin is structurally similar to isoquercetin, which is currently in clinical trials as possible medication for COVID-19.
Introduction
Coronaviruses (CoVs) are a family of viruses that cause several diseases in humans. These diseases range from common cold infections to severe acute respiratory syndromes [1, 2]. In recent times, coronaviruses have caused major outbreaks in the world. Among these are the 2003 severe acute respiratory syndrome coronavirus (SARS-CoV) outbreak in China, the 2012 Middle East respiratory syndrome coronavirus (MERS-CoV) outbreak in Saudi Arabia [3], and recently, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) outbreak that emerged from the Wuhan province of China in December 2019 [4]. According to recent studies, SARS-CoV-2 is a new strain of the virus and is responsible for the coronavirus disease 2019 (COVID-19) [5]. The fast spreading rate of the virus across the world prompted the World Health Organization (WHO) to declare the outbreak as a pandemic [6]. As of January 2022, case numbers above 300 million were reported worldwide with over 5.5 million deaths. When infected with the disease, symptoms such as dry cough, sore throat, fever, and tiredness are commonly experienced. Ageusia, anosmia, and runny nose are also common symptoms of the disease [7]. In more severe cases, patients bear serious symptoms such as difficulty in breathing, chest pains, among others [8]. In some instances, oxygen is administered through ventilators to augment patient’s survival [9, 10]. The severity of this disease has instigated the demand to identify new antiviral drug candidates or repurpose existing drugs for its treatment. Currently, a number of vaccines have been rolled out in the quest to provide immunity against the virus; however, availability and acceptance have become major issues as people continually get infected every day [11]. Thus, effective therapy is required, necessitating the use of medications as an additional therapeutic option in the fight against the disease.
CoVs share a lot of similarities in their genome organization [12]. Recent studies have shown that SARS-CoV-2 uses similar host cell receptors as SARS-CoV to enter human cells. This is attributable to the similarity in sequence of the spike proteins of both viruses. The entry of CoVs into host cells is facilitated by host-cell membrane fusion that involves several processes of receptor binding and proteolytic cleavage of the spike (S) protein. The spike (S) protein consists of two subunits namely S1 and S2 that are responsible for the attachment and entry of the virus into the host cell. While S1 is involved in the receptor binding on host cell surfaces, S2 facilitates the whole fusion mechanism. During the S protein proteolytic activity, several host cell proteases like furin, cathepsin B, trypsin, elastase, plasmin, and cell surface transmembrane protease/serine (TMPRSS) cleave the S protein to facilitate the viral entry [13]. Inhibition of these proteases has been suggested as plausible therapeutic targets for viral infections and can decrease viral infectivity [13]. Fusion between the virus and host cells occurs either via cytoplasmic or endosomal membrane fusion [14]. The first important step for target cell entry is the interaction of the virus spike protein with angiotensin-converting enzyme 2 (ACE-2) [15]. Alongside ACE2, the virus also employs the enzyme TMPRSS2 for priming [16].
Though the virus enters cells using these important receptors, emphasis has been placed on furin protease cleavage sites found on the virus spike protein for its role in facilitating its entry into host cells [17]. Furin belongs to a family of serine secretory proteases known as proprotein convertases (PCs) [18]. PCs are responsible for the regulation of majority of biological processes by activating precursor forms of a wide range of receptors, hormones, and cell surface protein [19]. In viral disease processes, furin and other PCs activate cell surface glycoproteins in the pathogenicity of several family of viruses including CoVs, paramyxoviruses, herpesviruses, togaviruses, bornaviruses, flaviviruses, bunyaviruses, filoviruses, orthomyxoviruses, retroviruses, and pneumoviruses, facilitating their entry into target cells [20,21,22]. Other pathogenic examples include furin activation of precursor proteins by influenza virus—a virus that has been shown to have several similarities with SARS-CoV virus, the Ebola virus, distemper virus, and many more [19]. These and other roles played by furin proteases make them important in the viral maturation process, SARS-CoV-2 pathogenesis, and viral transmission in humans. CRISPR-Cas9 knockout of furin has been reported to significantly reduce the production of infectious SARS-CoV-2 virus [23]. In ferrets, SARS-CoV-2 virus which lacks the furin cleavage site (FCS) was observed to have low transmission to other animals as compared to the wild-type virus [24]. Also, very low frequencies were observed for SARS-CoV-2 mutants that had FCS deletions in human tissues [24]. The spike FCS has a total of ten amino acid residues, of which the 682RRARSVAS689 region is highly conserved [25] and the 681PRRA684 region has been reported to be unique to SARS-CoV-2 [26]. Inhibiting the furin protease has been shown to prevent SARS-CoV-2 binding to the human furin protease, thereby suppressing viral production [22].
Drug discovery experiments in cell culture and animal models targeting the inhibition of furin protease as therapeutic intervention for specific diseases have been promising. In 2015, some novel furin inhibitors were tested via cell culture experiments against influenza virus, anthrax, and diphtheria toxins. Their findings showed that, in the presence of these inhibitors, the spread of the avian influenza viruses, H5N1 and H7N1, was strongly inhibited [19]. Anthrax and diphtheria toxins which are not viruses but depend on furin for their propagation showed signs of protective effect in the presence of the inhibitors [19].
Several types of synthetic inhibitors have been designed for furin protease since the human and mouse protein crystal structures were solved. Majority of these inhibitors are either peptide-based or non-peptide based derivatives such as nona-d-Arg-amide (D9R) and 2,5-dideoxystreptamine [27, 28]. However, only a few of these inhibitors have entered human clinical trials. Diminazene, an antiparasitic drug was identified to be a potent furin inhibitor via structure-based virtual screening and in vitro enzyme-based assay with an IC50 of 5.42 ± 0.11 μM [29]. Another study reported an IC50 of 13.2 μM for diaminazene against the furin protease [30]. It has also been reported that furin inhibitors, decanoyl-RVKR-chloromethylketone (CMK), and naphthofluorescein showed antiviral activity against the SARS-CoV-2 virus in VeroE6 cells by blocking viral entry and suppressing viral RNA transcription [22].
Few studies have targeted the inhibition of the enzyme using natural product sources [31]. Natural products and their associated moieties have historically been great sources of therapeutic agents [32]. Due to the resistance most viruses have towards antiviral therapy, there is a growing interest in natural products as one of the best resource for finding new chemically diverse leads that can be used to develop therapeutically new antiviral agents [33]. Natural products are rich sources of diverse chemical compounds from which drugs can be isolated [34]. Compared to synthetic drugs, drugs from natural sources have lower side effects, are less expensive, and are mostly less toxic [35]. In view of this, the exploration of its chemical space for scaffold purposes has been an essential part of drug discovery.
In this study, we aimed to identify putative antiviral inhibitors of furin protease using natural product-derived compounds of African origin as potential therapeutic agents for COVID-19 disease. We carried out consensus docking using pre-filtered libraries of African natural products against the binding site of the protein. The binding mechanisms, active site residue interactions, and binding free energies of the potential leads were evaluated using molecular dynamics simulations and the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations. The biological activity and pharmacological profiles of the compounds were predicted [36,37,38].
Methods
The study employed consensus docking approach by using OEDocking HYBRID (version 3.5.0.4) [39] and AutoDock Vina [40] for the molecular docking studies of furin protease and natural product-derived compounds from African sources (Fig. 1). The compounds with binding energies of − 7 kcal.mol−1 or less from both docking studies were considered [41]. The best poses of each compound from the two docking applications were compared and those with RMSD better than 2 Å were considered for downstream analysis (Fig. 1) [42]. The molecular interactions between the furin protease and the selected ligands were investigated using PyMOL (v 2.0.6) and Ligplot + [43]. Molecular dynamics simulations and MM-PBSA computations were also performed to support the selection of the potential lead compounds (Fig. 1).
Database preparation
A library of 7391 natural product compounds were obtained from three databases comprising Northern African Natural Products Database (NANPDB) [44], East African Natural Products Database (EANPDB) [45], and AfroDb database [34], a catalog of ZINC15 [46]. NANPBD and EANPDB have recently been merged to form the African Natural Products Database (ANPDB) [45].
The compounds were then filtered using FILTER (v4.0.0.4, OpenEye Scientific Software, Inc., Santa Fe, NM, USA; www.eyesopen.com) to eliminate all compounds with undesirable properties from the library especially non-drug-likeness, non-lead-likeness, and toxicity. FILTER is a molecular filtering application whose algorithm works based on physical property calculations and functional group knowledge. Its default drug/lead-like parameters were used for the filtration. After filtering, a total of 3942 compounds passed and were used for the virtual screening.
Protein structure retrieval and preparation
Crystal structure of the human furin protease was retrieved from the protein databank (PDB ID: 4RYD). The structure was solved using X-ray diffraction at a resolution of 2.15 Å [19]. The hexameric protein structure of 4RYD is composed of six similar chains namely A, B, C, D, E, and F, each having 482 amino acids sequence length [19, 47]. As a result of the similarity, chain A was selected for this study. Within the active site of the crystal structure was the bound ligand para-guanidinomethyl-Phac-R-Tle-R-Amba (MI-1148). The bound ligand and ions such as calcium and sodium found within the structure were removed using PyMOL (v 2.0.6).
Consensus docking
Consensus docking methodology is an approach that combines multiple docking programs by comparing their top scoring poses. A previous study [48] applied the consensus docking methodology on co-crystallized complexes using AutoDock [49], AutoDock Vina [40], and DOCK6 [50]. A success rate of 82% was observed by using more than one docking tool as compared to individual accuracies of 55%, 64%, and 58%, respectively [48]. Furthermore, because docking tools differ so greatly in their search and scoring algorithms, putting more emphasis on their intersection should compensate for their weaknesses [51,52,53]. OEDocking HYBRID [39] and AutoDock Vina [40] were combined in this study to predict potential inhibitors of the human furin protease. Only ligands which had binding energies equal to or less than − 7.0 kcal.mol−1 from both OEDocking and AutoDock Vina were selected for further studies since this threshold has been shown to distinguish between putative and non-putative binders of proteins [41]. Although, a more negative binding energy does not imply a better inhibition [54, 55], previous studies have shown that ~ 97.7% of known inhibitors have binding energies of − 7.0 kcal.mol−1 or less [41, 56] and this threshold filters ~ 95% of non-inhibitors. Other studies have also used the − 7.0 kcal.mol−1 threshold to prioritize compounds [57, 58].
Docking validation
HYBRID and AutoDock Vina were both validated before employed in the virtual screening process. A dataset of 52 furin active compounds and 224 decoys were used to assess the performance of the tools in the enrichment of actives for the furin protease. Fifty of the active compounds were retrieved from BindingDB [47] while the remaining two were the co-crystallized ligands obtained from complexes (PDB IDs: 4RYD and 4OMC). All 52 actives were submitted to Database of Useful Decoys (DUD-E) [59] and 224 decoys were obtained for them. For HYBRID, conformers of the compounds were generated using OMEGA [39, 60]. An OEDocking receptor was also created. An OEDocking receptor usually consists of the protein structure, its co-crystallized ligand, and descriptions of the active site. OpenEye Spruce4Docking (OpenEye Scientific Software, Inc., Santa Fe, NM, USA, www.eyesopen.com) utility program was used for generating the OEDocking receptor with default parameters. The OEDocking receptor of the furin protease was generated using the obtained complex (PDB ID: 4RYD) devoid of the other cofactors and water molecules. In other words, the furin protease and the bound ligand para-guanidinomethyl-Phac-R-Tle-R-Amba were used for the OEDocking receptor. The presence of the bound ligand guides the tool to select similar poses for the compounds used in the virtual screening process based on the reference ligand’s shape and 3-dimensional arrangement of its chemical features. HYBRID (version 3.5.0.4) was used to dock the conformer-generated active-decoy dataset into the active site of the produced receptor.
For AutoDock Vina, the set of actives and decoys were initially energy minimized before docking was carried out. The protein’s active site was defined using grid box dimensions set at X = 26.86 Å, Y = 27.66 Å, and Z = 0.38 Å for the center while box size dimensions were set at X = 23.31 Å, Y = 23.79 Å, and Z = 26.44 Å.
The validation of each of the docking protocols was assessed via Screening Explorer by evaluating the ROC curves and Boltzmann-enhanced discrimination of ROC (BEDROC) [61]. Results from this docking validation were used as guide for the subsequent virtual screening of our pre-filtered natural product library against furin in HYBRID and AutoDock Vina.
Virtual screening via OEDocking HYBRID
The first phase of the molecular docking process was performed using the OEDocking HYBRID (version 3.5.0.4) [39]. HYBRID has an improved scoring algorithm with a higher mean area under the curve (AUC) value of 0.78 than its predecessor, FRED (0.75) [39, 60]. Using the pose mode in OMEGA [62] with default parameters, different conformers were generated for the pre-filtered library of compounds. OMEGA implements an algorithm that breaks small molecules into fragments and reassembles them into 3D conformations subjecting each conformer to energy evaluations [62]. Conformers are finally clustered after meeting a particular energy threshold [62, 63].
The conformer-generated library was then docked into the active site of the prepared receptor using HYBRID 3.5.0.4 docking program of the OEDocking tools. The same OEDocking receptor prepared during the docking validation was used for the docking. The Chemgauss4 scoring function (OpenEye Scientific Software, Inc., Santa Fe, NM, USA, www.eyesopen.com) was used to score the best binding conformations. This was the first phase of the structure-based virtual screening in the consensus docking methodology.
Virtual screening with AutoDock Vina
The pre-filtered library of 3942 compounds was also virtually screened against the human furin protease using AutoDock Vina [64]. The structure-data file (SDF) format of the 3942 compounds was loaded into PyRx (version 0.8) [65] and energy minimized using its integrated Open Babel [66]. The energy minimization was done using the United Force Field (UFF) with a limit of 200 iterations for each ligand. The energy minimized ligands were each converted to PDBQT format, an acceptable input for AutoDock Vina. PyRx was again used to prepare the receptor into an acceptable AutoDock format. The same grid box dimensions used for the docking validation were selected for the screening of the hits. These dimensions alongside a specified exhaustiveness of eight were registered in a configuration file. The configuration file was retrieved along with the PDBQT files of the receptor and the ligands. The docking was performed on a Dell EMC high-performance computing (HPC) system which consists of CentOS 7 operating system, 6 nodes, 12 GPUs, 216 CPUs, and 277 TB of storage situated at the West African Centre for Cell Biology of Infectious Pathogens (WACCBIP), University of Ghana, Accra. After docking, nine different poses were generated for each compound. For each compound, the pose with the most negative binding energy and zero RMSD value was selected. Compounds which had binding energies equal to or below − 7 kcal.mol−1 were considered for further analysis [41].
Root mean square deviation comparisons
Docking poses of hit compounds obtained from AutoDock Vina screening were compared to their respective poses obtained from screening with HYBRID docking tool using the LigAlign [67] program interfaced with PyMOL (v 2.0.6). The root mean square deviation (RMSD) calculation was used to score how well poses align with each other. An RMSD cut-off of 2 Å was used [42] since this cut-off value is widely regarded as the most effective threshold value for validating correctly posed molecules [42, 48]. Previous studies have used this threshold to validate docking protocols [42, 68, 69], rank ligands [70], align ligands to reference structures [71, 72], and for comparative studies [73,74,75]. Herein, compounds with pose alignment scoring below the RMSD cut-off were retained for further analysis.
Selection of potential lead compounds
Compounds were selected as potential leads based on three criteria. Firstly, the compounds must have reasonably good binding energies from both docking protocols and be considered top hit after each screening. Secondly, their ability to align well in the pose comparison analysis, that is, compound alignment RMSD score must meet the cut-off of 2 Å [42]. Lastly, occupation of active site regions based on pose visualization analysis and compound interactions with critical residues (His194, Ser368, and Asp153) of the active site. The selected potential leads were then subjected to stability and binding free energy analysis using molecular dynamics (MD) simulation and MM-PBSA calculations in GROningen MAchine for Chemical Simulations (GROMACS) version 2018 [76].
Molecular dynamics simulation
The potential lead compounds in complex with the receptor were subjected to 100 ns MD simulation. This was done to study the stability of the interactions between the compounds and the receptor. GROMACS (version 2018) was used to carry out the MD simulations for the unbound protein and protein–ligand complexes [77, 78]. The GROMOS96 43a1 force field was employed for each simulation. The ligand topology files were generated using PRODRG2 server [79]. The complexes were solvated within an SPC water model in a dodecahedron box, with box-solute distance of 1.0 nm. Due to the negative charges on the protein, the system was neutralized using 14 sodium ions. The steepest descent algorithm was then used to energy minimize the system over 1000 steps. The system was equilibrated using the constant-temperature constant-pressure (NPT) and constant-temperature, constant-volume (NVT) ensembles. A 100 ns production MD was finally run on each complex system. The simulations were performed on the same computing resource used for the docking procedure. Xmgrace was used to plot the simulation graphs after the production run was completed [80].
MM-PBSA calculations
The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) binding energy calculations method was used to calculate the binding free energies of each protein–ligand complex of the potential leads with g-mmpbsa [81, 82] after they had been subjected to MD simulations via GROMACS (version 2018) [76] using the GROMOS96 43a1 force field. The total binding free energies of the complexes were calculated from the van der Waals, electrostatic, polar solvation, and non-polar solvation energies. The polar solvation and non-polar solvation energies are estimated from the Poisson–Boltzmann equation and the solvent accessible surface area (SASA) methods. The total binding energy (ΔGGbindΔ������) equation is given [83] as Eq. 1.
ΔGGbind=ΔGGcomplex−[ΔGGreceptor+ΔGGligand]Δ������=Δ���������−[Δ����������+Δ��������]
(1)
Structural similarity search and activity prediction of hits
A search for compounds similar in structure to the predicted hits was done via DrugBank [84] with a threshold of 0.7 in order to find close analogues of the predicted hits. Literature was then searched for the use of these analogues in the treatment of similar disease conditions or in any experimental studies relevant to our findings. Furthermore, the biological activities of the predicted hits were determined using Prediction of Activity Spectra for Substances (PASS) [37, 85]. PASS uses Bayesian models to predict biological activities of compounds and their mechanisms of action [37, 85]. Hence, PASS can estimate the likelihood of a compound belonging to a specific class of active compounds using the probable activity (Pa) and probable inactivity (Pi) measures, which range between 0 and 1. The SMILES of the predicted hits were used as inputs for the PASS predictions.
Pharmacokinetic and toxicity profiling of potential leads
The absorption, distribution, metabolism, excretion, and toxicity (ADME/Tox) properties of the potential leads were profiled using SwissADME [86], pkCSM [87], and ProTox-II [88]. It is crucial to determine the ADME/Tox profiles of compounds since exposure to some compounds and their combinations may be harmful [88]. In silico toxicity models complement existing experimental models by predicting the effects of compounds, which in turn reduce the time, cost, and number of animal models required for testing [88]. SwissADME employs multiple linear regression, BOILED-Egg, and various binary classification models for compound pharmacokinetic predictions [86]. Likewise, the pkCSM server uses distance-based graph modeling approach to predict and optimize the pharmacokinetic properties and toxicity of compounds [87]. Graph modeling is an approach where atoms of compounds are represented as nodes and their bonds as edges [89]. However, ProTox-II is solely for toxicity prediction of chemicals focusing mainly on cytotoxicity, acute toxicity, hepatotoxicity, carcinogenicity, adverse outcomes pathways (Tox21), mutagenicity, immunotoxicity, and toxicity targets. ProTox-II relies on molecular similarity, pharmacophores, fragment propensities, and machine-learning predictive models for the various toxicity predictions [88].
Results and discussion
The availability of a crystallized 3D furin protease structure and a bound ligand makes it convenient to employ structure-based virtual screening (SBVS) methods in this study. SBVS is a computational approach used at the early stage of drug discovery that makes it easier and faster to sample out potential bioactive compounds from a large library of compounds/chemical space. With SBVS, the best interaction mode between compounds and their targets are easily predicted [90]. In view of this, top-ranked hit compounds which bind to a target must achieve a particular conformation, position, and orientation (pose) in order to attain the desired interaction with residues of the protein active site [91]. Improvement of poses has become very important now in virtual screening because the accuracy of most scoring functions of docking algorithms depends on them. Due to this, the consensus docking methodology which is similar to the consensus scoring approach was proposed [48]. This study used the approach to predict inhibitors for the human furin protease, combining OEDocking HYBRID [39] and AutoDock Vina [40]. We compared poses of 45 hit compounds obtained through consensus docking from both docking programs. Twenty compounds had an RMSD below 2 Å, a threshold below which poses are said to align properly [42]. Subsequently, seven potential leads were selected based on their interactions with active site residues.
Docking validation
This study evaluated the docking performance of HYBRID and AutoDock Vina via Screening Explorer [61] in order to ascertain how effectively these docking tools can distinguish between active and inactive molecules for the furin target. Their performances were thus evaluated using metrics including the area under the curve (AUC) of the receiver operating characteristic curve (ROC) and Boltzmann-enhanced discrimination of ROC (BEDROC). The ROC AUC metric, with values ranging between 0 and 1, measures the relationship between sensitivity (fraction of selected true active compounds out of the total actives) and specificity (fraction of false positives out of the total inactive compounds). An AUC value close to 0 and less than 0.5 means poor discrimination [67, 68]. The greater the AUC, the higher the likelihood of the virtual screening protocol to discriminate between active and inactive compounds. AUC values of approximately 0.50 and 0.72 were achieved for HYBRID and AutoDock Vina, respectively.
Early recognition of active molecules by both docking protocols was assessed using BEDROC. BEDROC is an enrichment metric that uses an exponential function to assign more weights to early ranked molecules than late ones. BEDROC values range between 0 and 1, with 0.5 being the value expected for random selection. The BEDROC parameter α was set to 20.0 which means that 80% of the maximum contribution to the metric comes from 8% of the data list [67]. BEDROC values of 0.024 and 0.649 were observed for OEDocking and AutoDock Vina, respectively. BEDROC values greater than 0.5 are considered best performance for early enrichment. HYBRID had relatively low AUC and BEDROC values because the tool was unable to generate poses for 32 of the active compounds and 81 decoys using the bound MI-1148 as reference. The AUC value for AutoDock Vina which was above 0.7 (Fig. 2), implies that AutoDock Vina has a highly acceptable discrimination rate for our data compared to HYBRID.
HYBRID has the advantage of using the reference ligand to enhance the docking performance by selecting compounds with similar poses while AutoDock Vina has been predicted to have a high active-decoy discrimination rate. Thus using these two docking protocols complement the other’s strengths. Values obtained for the assessment analysis are suggestive that AutoDock Vina is a good docking tool for selecting actives against the protein target.
Molecular docking studies
Virtual screening with HYBRID
A library of 3942 pre-filtered natural product compounds was docked into the active site of the crystal structure of human furin protease. The spatial disposition of the known ligand (MI-1148) was used to define the active site by using Spruce4Docking (OpenEye Scientific Software, Inc., Santa Fe, NM, USA, www.eyesopen.com). A bound ligand is important when using HYBRID because its docking and scoring algorithm work based on how well docked molecules mimic the shape and 3D positioning of chemical features of the co-crystallized ligand in the protein active site. Docked poses of compounds are usually biased towards the pose of the crystallographic ligand. This is one reason why this docking application was employed for our initial docking protocol to obtain compounds which might dock with similar binding poses as MI-1148 (Fig. 3).
HYBRID 3.5.0.4 was used to run the initial molecular docking taking as inputs the receptor-ligand complex and the pre-filtered library. This formed the first phase of the consensus docking process. After this phase of virtual screening with HYBRID, compounds with binding energies below − 7 kcal.mol−1 were considered [41]. A total of 89 compounds were found in all, with the most negative binding energy being − 12.04 kcal.mol−1, implying the highest binding affinity. Upon visual inspections with VIDA (OpenEye Scientific Software, Inc., Santa Fe, NM, USA, www.eyesopen.com), all 89 compounds were found to dock in the active site. Compound ZINC000001530775 had the most negative binding energy with the furin protease with a value of − 12.035076 kcal.mol−1, followed by smirnovine and tubastrine with binding energies of − 10.141556 and − 9.564701 kcal.mol−1, respectively. Naringenin 7-p-coumaroylglucoside and ( +)-lariciresinol-9-O-beta-D-glucopyranoside had binding energies of − 9.401727 and − 9.345255 kcal.mol−1, respectively. Caulindole A, 6-bromo-N-methylaplysinopsin and pinobanksin_3-(E)-caffeate also had strong binding affinity to the furin protease with binding energies of − 9.247932, − 9.177244, and − 8.666552 kcal.mol−1, respectively.
Virtual screening with AutoDock Vina
The pre-filtered library was also docked against the furin protease using AutoDock Vina. A total of 3911 out of the 3942 compounds were successfully screened against the furin protease. In all, 1754 compounds had binding energies of − 7.0 kcal.mol−1 or lower. ZINC000095486208 demonstrated the most negative binding energy of − 10.0 kcal.mol−1 against the human furin protease. The known inhibitor, MI-1148, 10′-hydroxyusambarensine and 3-O-(beta-D-glucopyranosyl)_etioline had a binding energy of − 9.7 kcal.mol−1 when docked against the protease. Comparatively, a previous study docked 163 ligands against the furin protease and the TMPRSS2 using AutoDock Vina [92]. The top 23 compounds found from the study were reported to bind with the catalytic residues of furin at the active site with binding energies ranging from − 8.7 to − 7.0 kcal.mol−1 [92]. Also, in the quest to find SARS-CoV-2 inhibitors, 29 phytocompounds from Acacia pennata (L.) Willd were docked against the SARS-CoV-2 Mpro and the furin protease [93]. The 29 compounds had binding energies ranging from − 9.0 to − 2.0 kcal.mol−1 [93]. Another study docked mozenavir with the furin protease using AutoDock and reported a binding energy of − 12.05 kcal.mol−1 as against the reference furin inhibitor decanoyl-RVKR-chloromethyl ketone which had a binding energy of − 6.89 kcal.mol−1 [94].
Selecting hits from the docking protocol
Compounds which had met the − 7.0 threshold when docked using both OEDocking HYBRID and AutoDock Vina were shortlisted for further studies. From the top 89 and 1754 compounds from HYBRID and AutoDock Vina, respectively, 45 compounds demonstrated binding energies of − 7.0 kcal.mol−1 or lower with both docking protocols and were thus selected as potential hits. The binding energies of the top 45 compounds ranged from − 9.7 to − 7.0 kcal.mol−1. The compound, 10′-hydroxyusambarensine had the least binding with the furin protease with a binding energy of − 9.7 kcal.mol−1, followed by ZINC000095486212 with − 9.4 kcal.mol−1. ZINC000085967772 and apigenin-7-O-6″-E-p-coumaroyl-beta-D-glucopyranoside_beta_2 both had a binding energy of -9.3 kcal.mol−1. The compounds identified herein have relatively high binding affinity to the furin protease than those predicted using AutoDock Vina in previous studies [92, 93] and are worthy of further exploration. Furthermore, visualization of the selected hit compounds in PyMOL (v 2.0.6) showed that all the selected compounds docked well in the active site.
RMSD comparison
The poses of the top 45 compounds from AutoDock Vina were compared to their respective HYBRID poses. Out of the 45 selected compounds from the docking with AutoDock Vina, 20 of them had RMSD values ≤ 2 Å when compared to their corresponding HYBRID pose. The compound with the best aligned pose (malvinidin-3-arabinoside) had a RMSD score of 0.616, whereas pinobanksin_3-(E)-caffeate had the highest RMSD value of 1.993 Å among the top 20 compounds (Table 1). Compound 10′-hydroxyusambarensine which had the highest affinity to the furin protease (as predicted via AutoDock Vina) also had a RMSD of 0.772 Å (Table 1).Table 1 Top 20 potential hit compounds seeded for analysis after consensus docking of human furin protease. Shown in the table are their binding energy scores obtained from HYBRID and AutoDock Vina docking. RMSD values obtained from pose comparison are also presented
Furin-ligand interactions
Inhibitors which bind to the human furin protease were observed to bind at certain interaction sites on the catalytic domain of the protein (Fig. 4) [95]. These interaction sites within the active site contain the catalytic triad Ser368, His194, and Asp153 which are important for ligand interactions. Aside these catalytic residues, other ligand binding residues include Trp254, Pro256, Val231, Asn295, Glu257, Gly255, Ala267, Asp233, Asp236, Tyr308, and Asp264 [93, 94]. Mozenavir was previously shown in silico to interact with Tyr308, Gly265, Gly255, Asp154, and Val231 via hydrogen bonds and formed hydrophobic contacts with Glu236, Pro256, Trp254, Leu227, His194, Gly229, Asp264, Asp191, Glu230, and Gly229 [94]. Imatinib forms two hydrogen bonds with Glu236 and Gly255, and interacts with the furin protease via weak hydrophobic contacts with Val231, Pro256, Trp254, and Gly294 [29].
After visual inspection, most of the compounds occupied sites S1 to S4 (Fig. 4). None of the compounds was found to occupy position S5. Also, none of the compounds was found to occupy all the five sites. This may be due to the smaller sizes of the compounds when compared to the co-crystallized ligands.
A total of seven compounds were shortlisted based on the following criteria: (i) visual inspection of compound occupation of the active site mainly in the areas S1, S2, S3, S4, and S5; and (ii) compound interaction with important active site residues, with emphasis being placed on the catalytic triad. The seven shortlisted compounds include quercitrin, teucrol, malvinidin-3-arabinoside, N-E-caffeoyl tyramine, ZINC000085967772, pinobanksin 3-(E)-caffeate, and abyssinone IV. It is worthy to note that all seven selected potential leads (Fig. 5) occupied mostly S1 and S2 as these are the areas where majority of the interacting residues are found including the catalytic triad. Two compounds namely quercitrin and abyssinone IV were also found to form hydrogen bond interactions with critical residue Ser368 and hydrophobic bond interactions with His194 (Table 2 and Fig. 6A and B). The other five compounds did not share hydrogen bond interactions with any of the catalytic triad; however, they shared hydrophobic bond interactions with Ser368 and His194 (Figure S1). They also shared strong hydrogen bond contacts with other important surrounding residues (Table 2).
Table 2 Potential lead compounds and interacting residues via hydrogen bond (H-bond) interactions. Also shown is their binding energies obtained via AutoDock Vina comparable to that of MI-1148 which was used as control
Molecular dynamics simulation
Understanding the dynamic behavior of the proposed potential leads within the target’s binding pocket involves studying factors such as stability, structure compactness, and residue fluctuations. These factors were assessed via the root mean square deviation (RMSD) and radius of gyration (Rg) plots for 100-ns simulation time.
Root mean square deviation (RMSD)
The RMSD measures the flexibility/stability of the protein. This was monitored both for the unbound protein and the protein–ligand complexes. At the initial stages of the simulation, both the unbound protein and the complexes had their RMSD values peak to about 0.25 nm. The RMSDs of all the complexes rose until about 80 ns where they were observed to remain stable till the end of the 100-ns simulation time (Fig. 7). These results were consistent with the RMSD plots for furin-mozenavir, furin-folic acid, and furin-folinic acid complexes [94, 96]. However, same observation cannot be said about the unbound protein as the RMSD continued to rise even after 80 ns to levels around 0.3 nm. Similar observations have been reported in a 10-ns MD study, where the RMSD of the furin protease rose towards the end of the simulation time to about 0.4 nm [94]. This is indicative that the presence of the ligands may have brought stability to the protein. Nonetheless, complexes with ligands pinobanksin_3-(E)-caffeate (yellow colored) and teucrol (cream colored) achieved stability at RMSDs of 0.3 nm and higher.
Radius of gyration (Rg)
Rg is used to evaluate the compactness of the protein structures throughout the simulation period. A steady Rg for a particular structure indicates that the protein is stably folded [97, 98]. Most of the complexes including the unbound protein had relatively stable Rg values revolving between 2.15 and 2.2 nm (Fig. 8), consistent with previous Rg studies on furin-folic acid and furin-folinic acid complexes [96]. Only quercitrin and teucrol complexes had their Rg values below this range; however, they maintained a steady value of about 2.14 and 2.125 nm, respectively, beyond 60 ns (Fig. 8).
Binding free energy
MM-PBSA calculations have been deemed to be more accurate in estimating binding free energies of compounds than scoring functions of most docking programs [99]. Estimating binding free energies for protein–ligand interactions has been made robust using g-mmpbsa [81]. The g-mmpbsa integrates both GROMACS and APBS for its energy calculations. Herein, binding energies for each of the eight complexes were estimated separately after each undergoing molecular dynamics simulation. After the calculations, pinobanksin 3-(E)-caffeate had the strongest binding with a binding free energy score of − 189.219 ± 25.602 kJ/mol (Table 3 and Fig. 9). Comparatively, the same compound had a better binding affinity score from the molecular docking results with a value of − 9.1 kcal.mol−1 (Table 2). ZINC000085967772 which had the most negative binding energy (− 9.3 kcal.mol−1) in the docking results among the potential leads also showed strong binding from the MM-PBSA calculations with a binding free energy score of − 172.892 ± 24.913 kJ/mol (Table 3 and Fig. 9). All the other potential leads had reasonable binding free energies (Table 3) which correlated well with their binding affinities obtained after molecular docking. A previous study investigated the molecular interactions between the furin protein and two ligands comprising folic and folinic acids using the MM-PBSA method. Folic and folinic acids had binding free energies of − 27.90 and − 12.84 kcal.mol−1, respectively [96]. Although, the interacting residues of the furin protease with folic and folinic acids are different from those observed herein, the predicted compounds in this study have higher affinities to the furin protease than the two folate analogues [96].Table 3 Contributing energies of the potential leads estimated from molecular Mechanics/Poisson-Boltzmann Surface Area (MM-PBSA) calculations. The energy values are stated in average with their standard deviations
Herein, MM-PBSA calculations on the control (MI-1148) showed very strong binding with the target detailing a binding energy of − 2178.767 kJ/mol ± 67.643 kJ/mol. This is particularly not surprising as the compound was synthetically designed to specifically target the protein [19]. Nonetheless, the binding free energies exhibited by the predicted compounds are reasonable. Therefore, structural modifications or utilization as scaffolds for the design of furin-targeted compounds could elicit stronger binding with the furin protease.
Biological activity prediction and structural similarity search
Quercitrin, one of the potential leads is a glycoside formed from the flavonoid quercetin and deoxy sugar rhamnose. Quercetin has also been found to show therapeutic effects against Enterovirus 71 protease, therefore blocking the virus replication in hosts [100]. Enterovirus 71 is also a single-stranded RNA virus known to cause hand, foot, and mouth disease in under-aged children (< 5 years) [100]. Quercitrin is found in some medicinal plants used by folk medicine in the treatment of inflammation [101]. Additionally, it was predicted to be anti-inflammatory with Pa of 0.754 and Pi of 0.010. Inflammation is a key feature in COVID-19 disease, accounting for some unforeseen circumstances or even death in patients [102]. A study identified anti-inflammatory compounds which targeted p38 MAPK receptor in the quest to salvage high concentrations of pro-inflammatory cytokines in COVID-19 mechanisms [103]. Quercitrin was also predicted to be a membrane permeability inhibitor (Pa = 0.986 and Pi = 0.001), RNA synthesis inhibitor (Pa = 0.558 and Pi = 0.004), as anti-influenza virus (Pa = 0.683 and Pi = 0.007), and anti-herpes (Pa = 0.518 and Pi = 0.008).
Teucrol, abyssinone IV, and ZINC000085967772 were all predicted to be membrane permeability inhibitors with Pa values of 0.759, 0.644, and 0.714 and Pi values of 0.018, 0.065, and 0.033, respectively. Teucrol was also predicted to be a membrane integrity agonist (Pa = 0.948 and Pi = 0.004) with anti-inflammatory properties (Pa = 0.569 and Pi = 0.039). Abyssinone IV and ZINC000085967772 were predicted to be cytoprotectant (Pa = 0.694 and Pi = 0.005) and hepatoprotectant (Pa = 0.589 and Pi = 0.013). Abyssinone IV has antiviral properties against the rhinovirus (Pa = 0.600 and Pi = 0.006). Pinobanksin 3-(E)-caffeate which is a propolis component was predicted to be anti-inflammatory (Pa = 0.821, Pi = 0.005) and hepatoprotectant (Pa = 0.829, Pi = 0.004). Interestingly, propolis components have recently been studied as possible therapeutics for COVID-19, and it showed inhibitory effects on ACE2, TMPRSS2, and PAK1 signaling pathways [104, 105]. Propolis is also used in traditional medicine worldwide due to its reported biological activities which include antibacterial, antiviral, anti-inflammatory, and anticancer [106,107,108].
In addition, rosmarinic acid, which was part of the 20 selected hits, has also been shown to have anti-influenza effects when its activity was tested on mice infected with A/FM/1/47 H1N1 virus. Significant improvements were observed in the infected mice treatment [109]. Rosmarinic acid is a compound isolated from a plant known as Sarcandra glabra (Thunb.) Nakai. Pursuing these compounds for antiviral, anti-inflammatory, and cell protective activities hold promising future especially against single-stranded viruses like SARS-CoV-2. Therefore, the potential leads and some of the hit compounds could be explored via in vitro and in vivo studies for antiviral purposes.
Structural similarity search was performed for each of the selected compounds via DrugBank using a similarity threshold of 0.7. This procedure was performed to decipher if any of the potential leads have similar structures or analogues studied in vitro, in vivo, or in clinical tests. Quercitrin (quercetin 3-O-α-l-rhamnopyranoside) had a DrugBank similarity index (SI) of 1.0 with isoquercetin, a quercetin derivative which is monoglycosylated. Malvinidin-3-arabinoside which was also identified as a potential lead compound had a similarity index of 0.822 with isoquercetin. Isoquercetin is a drug that is currently being tested clinically as medication for moderate to severe COVID-19 patients [110]. This supports the possible evaluation of quercitrin with a broad antiviral activity as a potential COVID-19 therapeutic agent. Likewise, teucrol was structurally similar to amiloxate with similarity index of 0.722. Amiloxate is a derivative of cinnamic acid which is known to have anti-inflammatory properties [111]. Abyssinone IV and pinobanksin 3-(E)-caffeate were also found to be structurally similar to taxifolin, a compound known to be an anti-inflammatory agent [112, 113] with similarity indexes of 0.749 and 0.767, respectively. Abyssinone IV also has similarities to hesperetin (SI: 0.805). Hesperetin is known to have anti-cancer and anti-inflammatory effects [114, 115]. Pinobanksin 3-(E)-caffeate is also structurally similar to silibinin with SI of 0.767. Silibinin is a flavonolignan which is said to have hepatoprotective effects and is used in managing hepatitis [116, 117]. Hepatitis is caused by the hepatitis virus which is known to exploit furin family of proteins (proprotein convertases) in its disease mechanism [21].
Pharmacokinetic and toxicity profiling of potential leads
SwissADME and pkCSM predictions showed that all the potential lead compounds have reasonable solubility and intestinal absorption with an absorption rate greater than 45% (Table 4 and Fig. 10). Also, none of the compounds were predicted to cross the blood–brain barrier (BBB). This is specifically significant because earlier studies have established that drugs that target sections of the human body other than the nervous system should not cross the BBB in order to avoid psychotropic side effects [118]. The compounds were also favorable towards cytochrome P450 enzymes, having shown none or few possible inhibitions to CYP2D6, CYP2C19, CYP2C9, and CYP1A2. However, four of the ligands namely Teucrol, N-E-caffeoyl tyramine, Pinobanksin 3-(E)-caffeate, and Abyssinone IV were predicted as possible substrates of CYP3A4 enzyme. Hepatic and renal clearance were low for all the compounds. For the toxicity predictions, the compounds passed the Ames test for toxicity and none but N-E-caffeoyl tyramine was found to be hepatotoxic (Table 4). Protox-II toxicity predictions revealed that the compounds were immunotoxic with probability greater than 0.6. Since furin is induced in T cells and is crucial for maintaining peripheral immune tolerance [119], this may be expected. All the compounds, except for N-E-caffeoyl tyramine had plausible lethal doses (LD50) making them favorably less toxic (Table 5). The higher the lethal dose per kilogram of a compound, the lesser its toxicity. Three of the compounds namely Quercitrin, ZINC000085967772, and Pinobanksin 3-(E)-caffeate were predicted to have about 50% chance of being carcinogenic. Teucrol was predicted to have three other active targets with each having a probability above 0.7 (Table 5). However, the fact that most of these compounds were predicted to be within the toxicity classes of 4 and 5 gives strong indications that these compounds may nonetheless be non-toxic.Table 4 ADMET property predictions of potential leads using pkCSM server
Table 5 Toxicity prediction of potential lead compounds using Protox-II
Conclusion
The study used the consensus docking approach via the structure-based virtual screening of human furin protease to effectively identify seven potential novel anti-COVID-19 compounds consisting of quercitrin, teucrol, malvinidin-3-arabinoside, N-E-caffeoyl tyramine, ZINC000085967772, pinobanksin_3-(E)-caffeate, and abyssinone IV. This is in a bid to support ongoing efforts to identify effective therapeutics against the SARS-CoV-2 by targeting the furin protease [29, 30]. Furin protease is a plausible COVID-19 target due to its cleavage site on the spike protein and its role in facilitating the entry of SARS-CoV-2 into host cells. The molecules showed strong active site interactions with the catalytic residues of the furin with plausibly high binding affinity than previously identified compounds from earlier studies [92, 93]. They also had favorable results when subjected to MD simulations including MM-PBSA calculations, biological activity prediction, and structural similarity search [78, 94,95,96]. The compounds were predicted as antiviral, anti-inflammatory, anti-cancer, hepatoprotective, cytoprotective, RNA synthesis inhibitors, and membrane permeability inhibitors with reasonable ADMET properties. Isoquercetin, a structurally similar compound to quercitrin is currently undergoing clinical trials as COVID-19 drug. Once the biological activities of these compounds are reinforced experimentally, they could nonetheless serve as scaffolds for the design of new and improved furin protease inhibitors with potent antiviral properties.
Data availability
Not applicable.