An integrated virtual screening of compounds from Carica papaya leaves against multiple protein targets of SARS-Coronavirus-2.

The pandemic of SARS-Coronavirus-2 (Coronavirus-19) has been progressing by the increasing trend of the cases as well as deaths with neither vaccine nor drug is rationally used to stop the viral spread over. This study aims to perform an integrated virtual screening of compounds that had been identified from Carica papaya leaves, which are proposed to be a herbal treatment for SARS-Coronavirus-2. The screening was initiated by evaluating the 40 compounds from Carica papaya leaves for their drug-like likeness property. The selected compounds were then secondly screened using carcinogenic and toxicity filters. Further selected compounds were thirdly screened for their pharmacokinetic profile and the screening was lastly performed by docking the third selected compounds against multiple protein targets of SARS-Coronavirus-2 employing 3-chymotrypsin-like protease (3CLpro), papain-like protease (PLpro), RNA-dependent-RNA-polymerase (RdRp), endonuclease (EndoU), S1 and S2 region of spike protein. The results show that 20 of 40 compounds, which meet the requirements of drug-like likeness, carcinogenicity-toxicity filter, and pharmacokinetic profiles, can interact with the multiple protein targets of SARS-Coronavirus-2 with the order from high to low affinity as follows: S1 > 3CLpro > EndoU > RdRp > PLpro > S2. In conclusion, Carica papaya leaves are worth to be proposed for further in vitro study against SARS-Coronavirus-2 at both molecular and cellular levels.

Introduction

The outbreak of severe acute respiratory syndrome (SARS) Coronavirus 2019 (Covid-19) has been extending across the world along with the number of victims. As reported by WHO in 6th November 2020, there had been 219 countries, 47,596,852 cases, and 1,216,357 deaths affected by the coronavirus [1]. The fast transmission from human to human-made the victims increase day by day without any specific antiviral agents being applied to the infected patients [2]. The treatment using HIV antiviral agent (lopinavir) [3], antimalaria (chloroquine) [4], and anti-influenza (oseltamivir) [5] have been a little bit helping the urgent situation, however, those could be still in trial and error since no selective drug has been discovered, up to now. Currently, remdesivir is an approved re-purposed drug from ebola and marburg antiviral agent to Coronavirus-19 which is indicated for adults and adolescents (12 years old or older) with body weight at least 40 kg [6]. The viral infection also could not be prevented since the vaccine is still under assessment [7]. Although China’s Sinovac Biotech appeared to be safe in a late-stage clinical trial in Brazil, however, it still undergoes monitoring of the adverse side effects in the few months ahead [8].

The coronavirus is structurally made of an enveloped, positive sense, and single-stranded RNA that belongs to the family Coronaviridae [9]. Like many other coronaviruses (genera alpha, beta, and delta), the large replicase polyproteins pp1a and pp1ab are encoded by the partially overlapping 5′-terminal orf1a/b within the 5′ two-thirds of the genome is proteolytically cleaved into 16 putative nonstructural proteins (nsps) [10]. These putative nsps included two viral cysteine proteases, namely, nsp3 (papain-like protease) and nsp5 (chymotrypsin-like, 3C-like, or main protease), nsp12 (RNA-dependent RNA polymerase [RdRp]), nsp13 (helicase), and other nsps which are likely involved in the transcription and replication of the virus [9]. Other proteins of SARS-Coronavirus-2 that could be targeted are EndoU, S1, and S2. EndoU (nsp15) is a part of the SARS-Coronavirus-2 replicase-transcriptase system, responsible for virus replication and transcription system. S1 and S2 are regions of SARS-Coronavirus-2 spike protein that plays a key role in the receptor recognition and cell membrane fusion process. S1 domain can bind with human angiotensin-converting enzyme 2 (ACE2) receptor to initiate the fusion process by changing its conformation to pre-hairpin intermediate. This state enables the assembly of the fusion core in the S2 region of spike protein and bringing viral and cellular membranes into proximity for membrane fusion [11].

To date, there is no clinically approved inhibitor of the SARS protease yet they remain in development [12]. The protease inhibitor for a diverse virus has been existing [13] such as nelfinavir [14], amprenavir for HIV [15], and lopinavir-ritonavir [16] for HCV. However, the production of those protease inhibitors required a multistep reaction which is very expensive, while we need an emergency drug that is effective, but not costly.

A study in the past time found that biflavonoid amentoflavone from  Torreya nucifera showed significant inhibition towards SARS-Coronavirus 3CLpro with IC50 8.3 µM. At the subsequent study, flavonoids including apigenin, luteolin, and quercetin also demonstrated IC50 280.8, 20.2, and 23.8 µM, respectively, toward the protease proving that flavonoid is a potential scaffold for SARS-Coronavirus 3CLpro inhibitor [17]. Recently, flavonoid namely 5,6,7-trihydroxy-2-phenyl-4H-chromen-4-one (baicalein) has been co-crystallized with SARS-Coronavirus-2 3CLpro which confirmed the such compound’s binding site with the protease [18].

Flavonoid is known to have a pleiotropic effect meaning that is not only a single protein that can be targeted but also can affect the multiple protein targets in the one disease pathogenesis [19]. On the other hand, RdRp is an essential protease that catalyzes the RNA replication while encoded in the genomes of all RNA-containing viruses with no DNA stage. A compound mimicking biflavonoid namely theaflavin was marked to suppress SARS-Coronavirus-2 replication through inhibiting RdRp [20].

Papaya (Carica papaya) is one of the tropical fruits known containing amino acid, protein, carbohydrate, fiber, vitamin C, and other nutrients [21]. In particular leaves, there are at least 40 compounds (Table S1) identified in which flavonoid and its analogs present approximately 27% [22], [23], [24], [25]. These flavonoids and its analogs are including apigenin, catechin, deoxyquercetin, hesperitin, isorhamnetin, kaempferol, myricetin, naringenin, protocatechuic acid, quercetin, and rutin.

A molecular docking study had been conducted by Muhammad et al., presenting seven compounds from the whole  Carica papaya  trees against SARS-Coronavirus-2 3CLpro, PLpro, and RdRp [26]. In this present study, we performed  in silico  studies presenting 40 phytoconstituents of  Carica papaya  leaves against the two proteases (3CLpro and PLpro), RdRp, EndoU, S1, and S2 targeted protein of SARS-Coronavirus-2. The study was initiated by computationally screening those 40 compounds for their drug-like likeness, carcinogenicity-toxicity, and pharmacokinetics profiles using pkCSM online tools. The 20 hits from that were then simulated separately against 3CLpro, PLpro, RdRp, EndoU, S1, and S2 proteins of SARS-Coronavirus-2 using molecular docking.

Materials and methods

Materials

The protein model was using the SARS-Coronavirus-2 3D crystal structure of 3CLpro in complex with 5,6,7-trihydroxy-2-phenyl-4H-chromen-4-one with PDBID 6M2N [27], PLpro in complex with Ac-hTyr-Dap-Gly-Gly-VME (PDBID 6WX4) [28], RdRp in complex with cofactors (nsp7 and nsp8) (PDBID 6M71) [29], apo-EndoU (PDBID 6W01) [30], apo-S1 (PDBID 6VXX) [31], and apo-S2 (PDBID 6VSB) [32]. The 3D structure of the 40 ligands was downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/). The software being used was Biovia Discovery Studio 2020 (www.accelrys.com), AutodockTools1.5.6 (www.scripps.edu), AutodockVina which is embedded in PyRx version 0.8 (https://pyrx.sourceforge.io/), and pkCSM online tool (http://biosig.unimelb.edu.au/pkcsm/prediction). The hardware is with specifications as followed: HP Notebook 14 CM-0006-AU, processor AMD Ryzen 3 2200AU, HDD 1 TB, RAM 4 GB, and OS Windows 10.

Methods

The  in silico  prediction was initiated by screening 40 compounds of  Carica papaya  leaves collected from the PubChem database using MW, LogP, followed by AMES test and lastly human gastrointestinal absorption. The final screened compounds were then simulated using molecular docking against 3CLpro, PLpro, RdRp, EndoU, S1 dan S2 (see Fig. 1 ).

Fig. 1

The screening steps of 40 compounds from Carica papaya leaves in identifying the potential of this leaves to be a natural source for combating SARS-Coronavirus-2.

Drug-like likeness study

Forty ligands were downloaded from https://pubchem.ncbi.nlm.nih.gov/ and converted into SMILES files. Its Lipinski Rule properties were individually predicted by inputting its SMILES string and the prediction was done by the server. Instead of molecular weight and LogP were the filter representing the drug-like likeness, there are a few more parameters in this study including the number of hydrogen bond donor (HBD), the number of hydrogen bond acceptor (HBA), the number of rotatable bonds and the surface area, were also predicted.

Carcinogenic and toxicity studies

Using the same protocol in section 2.1, the carcinogenic profile of the selected compounds was represented by the AMES test result. Instead, other parameters such as maximum tolerated dose (human) (hMTD), hERG I inhibitor, hERG II inhibitor, oral rat acute toxicity (LD50), oral rat chronic toxicity (LOAEL), hepatotoxicity, skin sensitization,  T. pyriformis  toxicity, and minnow toxicity.

Pharmacokinetics study

Using the same protocol in 2.2, the pharmacokinetic profiles (absorption, distribution, metabolism, and excretion) of the selected compounds were used as the filter. Subsequently, the absorption is influenced by water solubility, Caco2 permeability, skin permeability, P-glycoprotein substrate, P-glycoprotein I inhibitor, and P-glycoprotein II inhibitor, instead of human gastrointestinal absorption. The distribution is represented by VDss (human), fraction unbound (human), blood–brain barrier (BBB) permeability, and central nervous system (CNS) permeability. The metabolism is represented by the CYP2D6 substrate, CYP3A4 substrate, CYP1A2 inhibitor, CYP2C19 inhibitor, CYP2C9 inhibitor, CYP2D6 inhibitor, and CYP3A4 inhibitor. Lastly, the excretion is represented by total clearance and renal OCT2 substrate.

Molecular docking study

Protein preparation

The 3D protein crystal structures were downloaded from www.rcsb.org and then uploaded in Discovery Studio. The complexed ligand from the individual protein target was taken out from the complex and the protein was saved as pdb file. The pdb protein file was then uploaded into Autodocktools1.5.6 and added by polar hydrogen followed by giving Kollman charge. The protein was then saved as pdbqt and ready for use.

Ligand preparation

The ligand for control docking was prepared by uploading it to AutodockTools1.5.6, given polar and non-polar hydrogen, Gasteiger charge and then saved as pdbqt. The available ligands in PubChem database (https://pubchem.ncbi.nlm.nih.gov/) were directly downloaded in its 3D structure followed by uploading them into Discovery Studio and then saved as pdb file. The optimized structure in PubChem server was reported in the published paper [33], [34]. Each pdb ligands were then uploaded into Autodocktools1.5.6 and then given gasteiger charges. The ligands were then saved as pdbqt and ready for use.

Control docking

The control docking was carried by redocking the individual native ligand into the respecting protein target with the parameters as followed: gridbox size (25×25×25) with the center of mass x = –33.718, y = −65.831, z = 41.2267; x = 9.4133, y = −28.3762, z = −38.0803; and x = 113.1457, y = 115.4740, z = 125.2150, for 3CLpro, PLpro, RdRp. The coordinate ligand being docked into EndoU, S1 and S2, respectively: x = 113.146, y = 115.474, z = 125.2150; x = 228.417, y = 208.3656, z = 256.074; and x = 203.530, y = 239.600 z = 186.327. The docking was run using Autodock Vina embedded in PyRx program with exhaustiveness = 64 covering 9 conformations for each ligand. The binding energy result was collected in csv file, whereas the best docking pose was selected and saved in pdbqt file ready for analysis. The control docking parameters were accepted when the RMSD value was not greater than 2.0 Å [35].

Virtual screening

The virtual screening through molecular docking of 20 ligands against 3CLpro, PLpro and RdRp, EndoU, S1, and S2 spike proteins of SARS-Coronavirus-2 was carried out using the same procedure of the individual control docking.

Analysis

The ligands were ranked according to the binding energy from the lowest to the highest value. Ten top lowest free energies of binding were then selected as the virtual hits. These virtual hits were then analyzed its molecular interaction with the binding site of the individual protein target using Discovery Studio 2020.

Results

The drug-like likeness is one of the criteria for a compound to be a drug candidate. This is well-known with the Lipinski Rule of Five postulating that a drug should be maximum having 500 g/mol in the molecular weight, <5 in the partition coefficient (log P), maximum of 5 in the number of hydrogen bond donor (HBD), and lastly maximum of 10 in the number of hydrogen bond acceptor (HBA) [36]. This rule is not 100% ensuring, however, it is guiding the drug design process. The log P is the ratio of the concentration of the compound in  n-octanol over its concentration in water, therefore it associates with the balance of the compound’s solubility in water during oral dissolution steps with its oral bioavailability of the compound in the blood system [36]. The log P < 5 is estimated to be the ideal value when the compound dissolves in our body fluid as well as its absorption through the gastrointestinal cell membrane, and then to be transported into the blood system. In this study, the first filter of 40 compounds from  Carica papaya  leaves was shortlisting 31 compounds having the maximum logP value less than 5 (Table S2). The second filter was, 23 of them passed the carcinogenicity effect by showing non-responsive toward the AMES test (Table S3). Finally, 20 compounds having at least 30% gastrointestinal absorption as the 3rd filter were shortlisted. Table 1 presents the 20 final selected compounds having log P are less than 5 along with other Lipinski Rule of Five criteria.

Table 1

The drug-like likeness profile of 20 compounds selected from Carica papaya leaves.

Ligands	Lipinski Rule				Rotatable Bonds	Surface area
	MW	log P	HBD	HBA
2S-sambunigrin	295.291	−0.93222	4	7	4	121.142
5,7-dimethoxycoumarin	206.197	1.8102	0	4	2	86.036
anthraquinone	208.216	2.462	0	2	0	92.536
apigenin	270.24	2.5768	3	5	1	112.519
ascorbic acid	176.124	−1.4074	4	6	2	67.321
caffeic acid	180.159	1.1956	3	3	2	74.381
caffeoyl alcohol	166.176	1.1033	3	3	2	70.219
catechin	290.271	1.5461	5	6	1	119.662
deoxykaempferol	270.24	2.5768	3	5	1	112.519
deoxyquercetin	286.239	2.2824	4	6	1	117.313
dimethoxyphenol	154.165	1.4094	1	3	2	65.183
ferulic acid	194.186	1.4986	2	3	3	81.065
kaempferol	286.239	2.2824	4	6	1	117.313
niacin	123.111	0.7798	1	2	1	51.972
p-coumaric acid	164.16	1.49	2	2	2	69.587
p-coumaroyl alcohol	150.177	1.3977	2	2	2	65.425
protocatechuic acid	154.121	0.796	3	3	1	62.341
riboflavin	376.369	−1.72356	5	9	5	152.292
R-prunasin	295.291	−0.93222	4	7	4	121.142
thiamine	265.362	0.60774	2	5	4	109.957

A drug should have MW which is lower than 500 Da, log P < 5, the number of HBD ≤ 5, the number of HBA ≤ 10, rotatable bonds ≤10, and surface area ≤140 Å [36], [37]. MW will affect the potency of the drug which is commonly expressed in IC50. The higher of MW would be the lower IC50, therefore, the more potent the drug. However, the MW should not be greater than 500 considering the drug permeability during intestinal absorption [36]. The number of HBD or HBA reflects their polarity to interact with water during the dissolution process as well as their molecular interaction during the pharmacodynamic step [36]. The rotatable bonds may influence their stability during pharmacokinetics and the receptor binding, thus, the less rotatable chain in the molecule should be the more stable drug to perform their activity [36], [37]. The polar surface area (SA) associates with the permeability of drugs across the cell membrane in which the higher SA might be poorer in cell permeability (oral bioavailability) [36], [37]. Therefore, catechin, deoxyquercetin, and deoxykaempferol are among compounds which meet the Lipinski rule for drug-like likeness.

A physical or chemical agent that exposes to an individual by causing cancer is named carcinogen, which in some carcinogenic agents, they are associated with increasing the risk of developing specific types of cancer [38]. For example, construction workers that are frequently exposed to asbestos, a carcinogenic agent, have been strongly linked to the development of a specific type of lung cancer called mesothelioma [39]. Importantly, it was identified some carcinogens from drug bearing phenacetin and azathioprine structures [40]. Table 2 presents the AMES test result of the 20 final selected compounds for carcinogenicity prediction along with other toxicity profiles.

Table 2

The AMES test result of the 20 final selected compounds for carcinogenicity prediction along with other toxicity profiles.

Ligands	Toxicity
	AMES*	hMTD	hERG I **	hERG II **	LD50	LOAEL	Hepato*	SS	Tp*	Minnow*
2S-sambunigrin	No	1.117	No	No	2.714	3.316	No	No	0.285	3.396
5,7-dimethoxycoumarin	No	0.711	No	No	2.137	2.59	No	No	0.609	1.06
anthraquinone	No	0.083	No	No	2.316	2.2	No	Yes	0.922	0.514
apigenin	No	0.83	No	Yes	2.423	1.753	No	No	0.386	1.114
ascorbic acid	No	1.987	No	No	1.434	3.376	No	No	0.285	3.612
caffeic acid	No	0.326	No	No	1.992	2.028	No	No	0.034	1.587
caffeoyl alcohol	No	0.289	No	No	1.994	1.959	No	Yes	0.326	1.941
catechin	No	1.072	No	No	2.261	3.212	No	No	0.285	1.833
deoxykaempferol	No	0.479	No	No	2.428	1.977	No	No	0.474	1.536
deoxyquercetin	No	0.837	No	No	2.651	1.592	No	No	0.304	1.263
dimethoxyphenol	No	1.345	No	No	1.809	2.866	No	No	0.252	1.212
ferulic acid	No	0.475	No	No	2.076	3.046	No	No	−0.011	1.492
kaempferol	No	0.676	No	No	2.698	1.658	No	No	0.3	1.407
niacin	No	1.168	No	No	2.008	2.862	No	No	−0.446	2.329
p-coumaric acid	No	0.338	No	No	2.099	2.908	No	No	0.01	1.739
p-coumaroyl alcohol	No	0.729	No	No	2.117	2.08	No	Yes	0.033	2.118
protocatechuic acid	No	0.607	No	No	1.951	2.341	No	No	−0.136	1.955
riboflavin	No	0.509	No	No	1.91	3.81	Yes	No	0.285	3.828
R-prunasin	No	1.117	No	No	2.714	3.316	No	No	0.285	3.396
thiamine	No	0.238	No	No	2.635	1.204	Yes	No	0.247	2.599

* = toxicity; ** = inhibitor

Toxicity predictions are used to predict whether these compounds having dangerous or toxic properties towards human physiology. AMES toxicity is employed to assess compounds carcinogenicity/mutagenic properties, therefore drug candidates should not be mutagenic or carcinogenic [38], [41]. The human maximum tolerated dose (hMTD) used to predict the value of the toxic dose threshold in human with the value of log hMTD for more than 0.477, are considered to be acceptable. hERG I and II are potassium channels that mediate the repolarizing of a cardiac action potential in humans. Inhibition of these proteins is the main cause of long QT syndrome development that might lead to fatal arrhythmia, hence drugs should not inhibit these ion channels. LD50 value represents the value of the dose given to cause 50% death of a group of rats and shows toxic potency of a compound in which a higher value of LD50 is considered to be safer. LOAEL (lowest-observed-adverse-effect level) value shows the lowest concentration of a compound to cause an adverse effect in human physiology indicated by the alteration of morphology, function, growth, or development. The safety of a compound improves as the LOAEL value increases. Hepatotoxicity shows the toxicity of a compound to cause liver injuries which disrupt its functions. It is expected that drug candidates should be non-hepatotoxic as considered to be a major safety factor in drug development. Potential dermal adverse effects are tested using a skin sensitization test, thus a compound should not induce allergic skin dermatitis.  T. pyriformis  and minnow toxicity are both utilized to measure the value of toxic endpoints. Drug environment safety is now taken into account to reduce the environmental damage by drugs with the acceptable value of  T. pyriformis  and minnow toxicity to be respectively higher than 0.5 and −0.3. Based on the overall toxicity predictions, 2S-sambunigrin, catechin, deoxyquercetin, dimethoxyphenol, and  R-prunasin are considered to have the best safety profiles among 20 compounds (http://biosig.unimelb.edu.au/pkcsm/prediction).

Before entering the binding site of the receptor, the drug will undergo pharmacokinetic steps including absorption, distribution, metabolism, and excretion (ADME) [42]. The molecule should be absorbed from the gastrointestinal into the blood vessel through the lipid bilayer membrane reflecting its bioavailability. Next, the molecule will be distributed by the blood and further carried out into the targetted cell bearing the receptor, while being interrupted by the plasma-protein binding. Some drugs will have the first pass effect by the liver’s enzymes which are early metabolized into either active or inactive drugs. Upon the receptor binding, the drug will give the therapeutic effect while gradually be metabolized by cytochrome P450 big family into their inactive metabolites. These metabolism products will be excreted throughout the body. Therefore, the prediction of pharmacokinetics profiles before docking using the software will help to design a compound that in one hand meets the ADME criteria, and on the other hand, it has a good binding with the receptor. Table 3, Table 4, Table 5 and Fig. 2 present the pharmacokinetics profile of 20 hits selected from 40 compounds identified from the  Carica papaya  leaves using pkCSM online tool.

Table 3

The absorption profiles of the 20 final selected compounds as predicted by the software.

Ligands	Absorption
	Water solubility (log S)	Caco2 permeability	Intestinal absorption (human)	Skin Permeability	P-gp substrate	P-gp I inhibitor	P-gp II inhibitor
2S-sambunigrin	−3.574	0.212	40.072	−2.744	No	No	No
5,7-dimethoxycoumarin	−3.343	1.281	99.221	−2.401	No	No	No
anthraquinone	−2.387	1.041	99.925	−2.365	Yes	No	No
apigenin	−2.079	1.136	91.403	−2.736	Yes	No	No
ascorbic acid	−0.429	−0.395	39.716	−3.478	No	No	No
caffeic acid	−1.773	−0.046	55.525	−2.735	Yes	No	No
caffeoyl alcohol	−1.538	1.091	68.137	−3.047	Yes	No	No
catechin	−2.910	−0.225	62.740	−2.735	Yes	No	No
deoxykaempferol	−3.405	1.139	92.943	−2.735	Yes	No	No
deoxyquercetin	−2.987	0.109	84.972	−2.735	Yes	No	No
dimethoxyphenol	−2.101	0.841	95.368	−2.504	No	No	No
ferulic acid	−1.719	0.100	94.737	−2.709	Yes	No	No
kaempferol	−2.986	−0.053	84.952	−2.735	Yes	No	No
niacin	−0.515	1.219	86.389	−2.790	No	No	No
p-coumaric acid	−1.610	1.142	93.183	−2.566	Yes	No	No
p-coumaroyl alcohol	−1.325	1.491	92.284	−2.815	Yes	No	No
protocatechuic acid	−1.675	0.805	67.889	−2.736	Yes	No	No
riboflavin	−1.952	−0.385	54.626	−2.740	Yes	No	No
R-prunasin	−2.079	0.212	40.072	−2.744	No	No	No
thiamine	−2.903	0.872	92.302	−2.963	No	No	No

Table 4

The distribution profile of 20 compounds identified in Carica papaya leaves as predicted by the software.

Ligands	Distribution
	log VDss (human)	Fraction unbound (human)	BBB permeability (log BB)	CNS permeability (log PS)
2S-sambunigrin	−0.706	0.439	−0.929	−4.57
5,7-dimethoxycoumarin	−0.141	0.37	0.274	−2.405
anthraquinone	0.251	0.118	0.252	−1.474
apigenin	−0.361	0.143	−1.082	−2.195
ascorbic acid	−0.264	0.808	−1.233	−4.332
caffeic acid	−0.409	0.358	−0.851	−3.326
caffeoyl alcohol	0.011	0.514	−0.347	−2.526
catechin	0.675	0.187	−1.182	−3.449
deoxykaempferol	0.354	0.148	−0.864	−2.092
deoxyquercetin	0.141	0.112	−1.219	−2.378
dimethoxyphenol	−0.112	0.486	−0.133	−2.092
ferulic acid	−0.535	0.397	−0.263	−2.93
kaempferol	−0.178	0.087	−1.361	−2.357
niacin	−0.794	0.601	−0.353	−2.928
p-coumaric acid	−0.602	0.421	−0.234	−2.379
p-coumaroyl alcohol	0.124	0.464	−0.183	−1.835
protocatechuic acid	−0.332	0.373	−0.985	−3.333
riboflavin	−0.818	0.558	−1.743	−5.243
R-prunasin	−0.706	0.439	−0.929	−4.57
thiamine	0.474	0.574	−0.441	−3.022

Table 5

The interaction between 20 ligands identified from Carica papaya leaves with a diverse CYP subfaimilies.

Parameter	Metabolism
	CYP2D6 substrate	CYP3A4 substrate	CYP1A2 inhibitor	CYP2C19 inhibitor	CYP2C9 inhibitor	CYP2D6 inhibitor	CYP3A4 inhibitor
2S-sambunigrin	No	No	No	No	No	No	No
5,7-dimethoxycoumarin	No	No	Yes	No	No	No	No
anthraquinone	No	Yes	Yes	No	No	No	No
apigenin	No	No	Yes	Yes	Yes	No	Yes
ascorbic acid	No	No	No	No	No	No	No
caffeic acid	No	No	No	No	No	No	No
caffeoyl alcohol	No	No	No	No	No	No	No
catechin	No	No	No	No	No	No	No
deoxykaempferol	No	No	Yes	Yes	Yes	No	No
deoxyquercetin	No	No	Yes	No	Yes	No	Yes
dimethoxyphenol	No	No	Yes	No	No	No	No
ferulic acid	No	No	No	No	No	No	No
kaempferol	No	No	Yes	No	Yes	No	No
niacin	No	No	No	No	No	No	No
p-coumaric acid	No	No	No	No	No	No	No
p-coumaroyl alcohol	No	No	Yes	No	No	No	No
protocatechuic acid	No	No	No	No	No	No	No
riboflavin	No	No	No	No	No	No	No
R-prunasin	No	No	No	No	No	No	No
thiamine	No	No	Yes	No	No	No	No

Fig. 2

The histogram plot of 20 compounds identified from Carica papaya leaves against their total clearance.

According to Table 3, 17 compounds demonstrate more than 50% human intestinal absorption, in which 5,7-dimethoxycoumarin, anthraquinone, apigenin, deoxykaempferol, dimethoxyphenol, ferulic acid,  p-coumaric acid,  p-coumaroyl alcohol, and thiamine are among having the best absorption profile for more than 90%. The drug properties should have water solubility represented as log S value is higher than −4, therefore, all compounds are predicted among the good drug properties during the dissolution step. Caco2 is also the  in vitro  cell model to predict the absorption of an orally administered drug [43]. The drug should have Caco2-permeability for more than 0.90, thus, only eight compounds are predicted to have a high human gastrointestinal absorption. For a transdermal route of administration, a compound should have a skin permeability for at least −2.5 or lower. Among 20 compounds, 18 of them are suitable for a transdermal dosage form (excluding 5,7-dimethoxycoumarin and anthraquinone). P-glycoprotein (P-gp) is a protein transport essential during pharmacokinetics steps. This could have either advantage in therapeutic effect or even the contradictive ones [44]. A compound should have no P-gp inhibition for either P-gp I or P-gp II. Therefore, no compounds are predicted to either slow down the therapeutic effect or its contraindication. In conjunction, this should be proportional with their P-gp substrate inhibition profiles.

During the distribution, a drug should have a steady state uniform concentration while reaching up the whole tissues rather than plasma [45]. This is defined by the number of VDss ≥ −0.15. Among 20 compounds, 8 compounds (anthraquinone, caffeoyl alcohol, catechin, deoxykaempferol, deoxyquercetin, dimethoxyphenol, p-coumaroyl alcohol and thiamine) are the ones which meet this criteria. Fraction unbound describes the amount of drug which is free from plasma protein binding associated with the total concentration ready for the receptor binding. This should have the value of ≥0.15, therefore, except antraquinone, apigenin, deoxykaempferol, and deoxyquercetin, the remained compounds are available for the receptor binding which may increase the drug activity. The BBB permeability defines the possibility of compound to cross the brain membrane which may affect the CNS [46]. For a safe drug, these values should be <−1 and <−3 (poorly distributed to the brain and unable to penetrate CNS) for BBB and CNS permeability, respectively. Therefore, compounds like 5,7-dimethoxycoumarin and anthraquinone should be taken in a good control to avoid the CNS either depression or excitation. Table 4 presents the distribution profile of 20 compounds identified in Carica papaya leaves.

Drug metabolism occurs mainly in the liver which will undergo chemical alteration to improve its excretion or elimination. There are main six subfamilies (CYPs) that are being responsible for most of human drug metabolism with CYP1A2, CYP2C9, CYP2C19, CYP2D6, CYP3A4, are some of the most important within these subfamilies [47]. CYP2D6 and CYP3A4 are the most significant enzymes in overall drug metabolism, and they are also found in extrahepatic sites such as in the brain for CYP2D6, which is responsible for drug activity in the brain. Other example is CYP3A4 in the intestines which is mainly affecting oral bioavailability by first-pass metabolism. Lower bioavailability and activity are expected for compounds that have a high affinity towards both enzymes as substrate. CYP450 enzymes activity can be modified by drugs, either inhibition or activation. This is potentially causing clinical drug-drug interactions, leading to adverse reactions or therapeutic failures. Therefore, inhibition of these enzymes increases the possibility of adverse drug reactions and interactions. Most likely, all compounds do not act as the substrate for CYP2D6 and CYP3A4, except anthraquinone which may act as the substrate of CYP3A4. More compounds are likely inhibit the CYP1A2 but less inhibit other CYP as shown in Table 5.

Drug excretion is a process of either eliminating or removing drugs directly in unchanged form or its inactive form. Most drugs are excreted by the kidneys and drug excretion becomes less efficient and dosing adjustments may be needed as kidney function declines. Total clearance is the rate at which a compound is removed from the body i.e. excreted in the urine as a compound with a higher clearance value leads to shorter half-time and hence its effective duration [42]. Thiamine is the fastest compound to be eliminated from the body due to its highest total clearance. In contrast, caffeoyl alcohol is the slowest compound to be eliminated from the body due to its lowest total clearance.

Drugs may be passively excreted by the kidney through glomerular filtration or actively by tubular secretion. OCT2 transporter is one of the main renal uptake transporter to actively remove drug from blood. It plays a key role in the removal and renal clearance of mostly cationic drugs and endogenous compounds [48]. Inhibition of OCT2 (such as by cimetidine) decreases OCT2-dependent renal clearance drugs, such as metformin, hence altering pharmacokinetics and pharmacodynamics profiles that may lead to undesireable adverse effects. Interestingly, none compounds having inhibition towards renal OCT2 substrate that might be not having undesirable side effects. Based on the prediction, every ligands excretion will not depend on the renal OCT2 transporter.

The control docking results for 3CLpro and PLpro against their co-crystallized ligands are considered to be acceptable with a RMSD value of its internal ligand being 0.7146 Å and 1.2141 Å, respectively. On the other hand, the other proteins have no co-crystallized ligand, therefore, in this case, no control docking has been carried out for internal validation. However, we tried to do an external validation using remdesivir as this drug is primarily targetting RdRp, Glisoxepide for EndoU [49],  N-(9,10-dioxo-9,10-dihydroanthracene-2-yl)benzamide for S1 spike protein [50], and captopril for S2 spike protein [51] of SARS-Coronavirus-2. PyRx has a facility to be a site finder that predicts the binding site of apoprotein. The binding affinity of remdisivir into RdRp SARS-Coronavirus-2 is −6.6 kcal/mol (catalytic site), and −8.6 kcal/mol (allosteric site) representing that the parameters are capable to calculate the binding affinity of the selected 20 compounds from  Carica papaya  leaves. Table S4 presents the control docking results of 3CLpro, PLpro, RdRp, EndoU, S1 dan S2 protein. At the same time, control docking of glisoxepide (EndoU), N-(9,10-dioxo-9,10-dihydroanthracene-2-yl)benzamide (S1 spike), and captopril (S2 spike) demonstrate the binding affinity of the ligand into the corresponding target as follows: −8.2, −9.7, and −4.5 kcal/ mol, respectively. Table 6 lists down the binding affinity of the selected 20 compounds from  Carica papaya  leaves.

Table 6

The binding affinity of 20 compounds identified from Carica papaya leaves upon molecular docking study into 3CLpro, Plpro, RdRp, EndoU, S1 and S2 spike of SARS-Coronavirus-2.

Ligands	Binding Affinity (Kcal/mol)
	3CLpro	PLpro	RdRP	EndoU	S1	S2
2S-sambunigrin	−7.0	−6.0	−6.5	−6.5	−7.3	−5.3
5,7-dimethoxycoumarin	−5.8	−5.9	−5.4	−5.4	−6.1	−5.4
anthraquinone	−6.5	−5.7	−6.4	−6.8	−7.0	−5.9
apigenin	−7.2	−6.6	−6.5	−7.2	−7.7	−6.3
ascorbic acid	−5.0	−4.9	−5.4	−4.7	−5.6	−4.8
caffeic acid	−5.8	−5.2	−5.7	−6.1	−6.1	−4.9
caffeoyl alcohol	−5.4	−4.9	−5.1	−6	−5.9	−4.8
catechin	−7.5	−6.3	−6.9	−7.7	−7.8	−6.1
deoxykaempferol	−7.4	−6.8	−6.5	−7.2	−7.8	−6.0
deoxyquercetin	−7.8	−6.7	−6.7	−7.5	−8.0	−6.0
dimethoxyphenol	−4.7	−4.4	−4.7	−5.7	−4.9	−4.2
ferulic acid	−5.6	−5.2	−5.4	−5.7	−6.1	−4.8
kaempferol	−7.2	−6.7	−6.9	−7.2	−8.2	−5.9
niacin	−4.5	−4.6	−4.6	−4.6	−4.9	−4.7
p-coumaric acid	−5.3	−4.8	−5.2	−5.6	−5.9	−4.6
p-coumaroyl alcohol	−5.0	−4.7	−4.7	−5.2	−5.3	−4.4
protocatechuic acid	−5.1	−4.9	−5.6	−5.1	−5.8	−4.6
riboflavin	−7.4	−7.5	−6.6	−6.6	−8.1	−6.2
R-prunasin	−6.8	−6.1	−6.3	−6.6	−7.1	−5.6
thiamine	−5.8	−6.1	−6.0	−5.3	−7.0	−5.1
Mean	−6.3	−5.7	−5.9	−6.1	−6.6	−5.3

From the docking results, 20 ligands show affinity with 3CLpro, PLpro, RdRp, EndoU, S1 and S2 ranging in −4.5 to −7.8 kcal/mol, −4.6 to −7.5 kcal/mol, −4.6 to −6.9 kcal/mol, −4.6 to −7.7 kcal/mol, −4.9 to −8.2 kcal/mol, and −4.2 to −6.3 kcal/mol, respectively. These describe the most favorable interactions of the 20 ligands are to 3CLpro, EndoU, and S1, rather than PLpro, RdRp, and S2. This is in agreement with the average of the binding affinity, in which interaction with 3CLpro (−6.3 kcal/mol), EndoU (−6.1 kcal/mol), and S1 (−6.6 kcal/mol) has the lowest binding energy as mentioned in parentheses. Fig. 3 illustrates the overlapping 20 ligands docked conformation in the binding site of individual 3CLpro, PLpro, RdRp, EndoU, S1, and S2.

Fig. 3

The overlapping docked pose of 20 compounds identified from Carica papaya leaves in the binding site of a) 3CLpro, b) PLpro, c) RdRp, d) EndoU, e) S1, and f) S2. The proteins were presented in a surface model, whereas the ligands were presented in a yellow stick model, with C, H, O, and N are colored by yellow, white, red and blue, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In particular, on one hand, dexoyquercetin demonstrates the lowest binding energy toward 3CLpro (-7.8 kcal/mol), whereas catechin performed the most favorable binding toward EndoU. Furthermore, kaempferol shows the lowest binding energy toward S1 (-8.2 kcal/mol). On the other hand, riboflavin, kaempferol, and apigenin demonstrate the lowest binding energy towards PLpro, RdRp and S2, respectively. Fig. 4, Fig. 5 illustrates the best binding pose of ligands into individual protein binding sites.

Fig. 4

The docking pose of a) deoxyquercetin (−7.8 kcal/mol) into 3CLpro, b) riboflavin into PLpro (−7.5 kcal/mol), and c) kaempferol (−6.9 kcal/mol) into RdRp. The proteins were visualized in ribbon model, whereas the ligands were presented in a yellow stick model, with C, H, O, and N are colored by yellow, white, red and blue, respectively. The H-bond interaction and the hydrophobic interaction are presented in black and orange dashed lines, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5

The docking pose of a) catechin (−7.7 kcal/mol) into EndoU, b) kaempferol into S1 (−8.2 kcal/mol), and c) apigenin (−6.3 kcal/mol) into S2. The proteins were visualized in ribbon model, whereas the ligands were presented in a yellow stick model, with C, H, O, and N are colored by yellow, white, red and blue, respectively. The H-bond interaction and the hydrophobic interaction are presented in black and orange dashed lines, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The binding pose of deoxyquercetin in the 3CLpro binding site occurs through the H-bond interaction with GLU166 and PHE140. Instead, the ligand pose might be stabilized by the hydrophobic interaction with HIS163, MET165, and GLU166. In the binding pose of riboflavin in PLpro binding site, the H-bond interaction occurred at ASP302, VAL165, ARG166, TYR273, and SER245, whereas the hydrophobic interactions were performed by interacting with PRO247, TYR268, and ARG166. Kaempferol interacts with LYS545, THR556, ARG553, and ARG624. The hydrophobic interactions of this complex occur at the residues identical to H-bond interactions.

The H-bond interactions were further shown by catechin while interacting with HIS250, ASN278, LYS290, SER294, and LEU346 of EndoU. Besides, the hydrophobic interactions were performed by interacting with TYR343 and LYS345. Kaempferol was again showing the best affinity with SER371, ARG408, and GLN409 of S1 protein. The only hydrophobic interaction was performed by interacting with PRO384. S2 protein binding site was occupied by apigenin in its best binding affinity among 20 ligands performing H-bond interactions with LYS933, THR719, SER929, and ALA930. There is no hydrophobic interaction identified in this binding pose. The detail information about distance and angle of individual H-bond interactions is presented in Table 7 . H-bonds can be classified into weak, moderate and strong H-bonds depending on its angle and distance. Strong H-bonds have the distance range from 2.2 to 2.5 Å and angle between 170 and 180 degree, whereas weak H-bonds have the distance more than 3.2 Å and the angle more than 90 degree. Moderate H-bonds have the distance between strong and weak H-bonds and the angle more than 130 degree [52].

Table 7

The detail information about distance and angle of individual H-bond interactions. The atom type and numbering of each ligands are defined in Table S5.

Targets	Ligands	Interactions	HB distance (Å)	HB angles (°)	Remark
Mpro	deoxyquercetin	GLU166:NH – DOQ:O1	2.235	152.940	Moderate
		DOQ:H29 – PHE140:O	2.929	91.570	Weak
PLpro	riboflavin	VAL165:NH – RF:O6	3.050	117.750	Weak
		ARG166:NH – RF:O6	2.317	159.880	Moderate
		TYR273:H – RF:N9	2.667	92.620	Weak
		RF:H46 – ASP302:OD1	2.323	144.660	Moderate
		SER245:HB2 – RF:O5	2.567	143.085	Moderate
RdRp	kaempferol	LYS545:HZ3 – KF:06	1.970	167.720	Moderate
		THR556:HG1 – KF:O4	2.910	115.190	Weak
		ARG624:HH12 - KF:O4	2.367	107.920	Weak
		ARG624:HH22 - KF:O4	2.283	108.620	Weak
EndoU	catechin	LYS290:HZ3 – CTC:O5	1.960	172.756	Strong
		SER294:NH – CTC:O2	2.042	160.573	Moderate
		LEU346:NH – CTC:O3	2.108	154.170	Moderate
		CTC:H34 – HIS250:NE2	2.312	129.033	Weak
		ASN278:HD22 – CTC:O3	2.208	137.886	Moderate
S1	kaempferol	GLN409:HE21 – KF:O4	2.820	94.092	Weak
		ARG408:HH12 – KF:O2	2.336	135.509	Moderate
		ARG408:HH22 – KF:O2	2.247	137.850	Moderate
S2	apigenin	LYS933:HZ3 – APG:O1	2.225	147.531	Moderate
		THR719:HB – APG:O3	2.615	136.176	Moderate
		SER929:HB2 – APG:O4	2.548	119.748	Weak
		ALA930:HA – APG:O1	2.892	132.931	Moderate

Discussion

Carica papaya leaves have been traditionally used to relieve dengue fever in some Asian countries. The leaves are sliced into smaller pieces and followed by boiling them into the water for at least 15 min and then filtered out to collect the liquid phase. This 30 mL of the aqueous extract is three times daily used for dengue patients until the fever fully recovered into a normal body temperature [53].

Scientifically, the methanolic extract of  Carica papaya  leaves showed cytotoxic effects (CC50 = 0.6156 mg/mL) to LLC-MK2 cells and it showed inhibitory activity (EC50 ≥ 1 mg/mL) against DENV-2 with a selectivity index value of ±>1 [54]. Treatments with 500 mg/kg and 1000 mg/kg of freeze-dried  Carica papaya  leaf juice increased the platelet and leukocyte counts in DENV2-infected AG129 mice. A significant decrease (p < 0.05) in viral RNA level was detected in the liver and kidney of infected AG129 mice treated with 1000 mg/kg of freeze-dried  Carica papaya  leaf juice.  Carica papaya  treatment also significantly decreased (p < 0.05) the levels of certain cytokines and chemokines in plasma, liver, and kidney tissues of infected AG129 mice [55]. A study by Sharma et al. suggested that papaya leaves extract significantly decreases the expression of the envelope and NS1 proteins in DENV-infected THP-1 cells. This marked a decrease in intracellular viral load upon the extract treatment confirmed its antiviral activity [56].

Studies were carried out on natural papain inhibitor from papaya latex. The isolated fractions, identified as inhibitors I and II, showed a negative reaction with ninhydrin; however, the fraction identified as P-III showed a positive reaction with ninhydrin. Kinetics data showed non-competitive inhibition (inhibitor I) and uncompetitive (inhibitors II and P-III) [57].  In silico  anti-dengue activities of the extracts from  Carica papaya  by using bioinformatics tools were investigated. Interestingly, the flavonoid quercetin performed the highest binding energy against NS2B-NS3 protease which is evidenced by the formation of six hydrogen bonds with the amino acid residues at the binding site of the receptor [58]. This is then later proven by  in vitro  study, that flavonoid quercetin able to inhibit the DENV2 NS2B-NS3 protease through a non-competitive inhibition [59]

In this present study, an integrated  in silico  screening has been performed against 40 compounds identified from  Carica papaya  leaves. The 40 compounds have a diverse scaffold bearing flavonoid, flavanol, alkaloid, phytosterol, glycoside, phenylpropanoid, carbohydrate, coumarin, anthraquinone, and flavin class compounds. These compounds have their drug-like likeness, carcinogenicity, toxicity, pharmacokinetics, and pharmacodynamic profiles that contribute to the therapeutic effect of  Carica papaya  leaves as discussed above. There had been 20 of those 40 compounds finally selected, which might contribute to the therapeutic effects better than its toxic properties.

These 20 compounds were predicted to meet the rule of Lipinski, therefore, they should be stable during oral administration. These selected compounds should not be carcinogenic, cardiotoxic, hepatotoxic, and having drug tolerance that lowering the drug safety index. These compounds also should be well absorbed by the biological membrane and then well distributed into the site of action. Upon therapeutic effects, these compounds should be easily metabolized into an inactive metabolite and followed by the excretion of chemicals from our body system to minimize chemical retention that could be having adverse drug reactions. Finally, these 20 compounds demonstrated molecular interactions upon  in silico  docking them into the diverse protein targets of SARS-Coronavirus-2. According to the binding energy average, the 20 compounds interact with the protein in a sequence of higher to lower binding affinity as follows: S1 > 3CLpro > EndoU > RdRp > PLpro > S2. This means that most likely,  Carica papaya  might disrupt the SARS-Coronavirus-2 life cycle by interrupting the binding of spike protein into ACE2 receptor, the proteolysis of nsp4-13, and the reverse RNA transcription rather than inhibiting the RNA replication, proteolysis of nsp1-3, and the fusion core in S2 region of spike protein and bringing viral and cellular membranes into proximity for membrane fusion.

The notable compounds' performance in their protein target is most likely having flavonoid/ flavanol/ flavin scaffold. This is for sure in correspondence with the previous studies that flavonoid class compound has a pleiotropic effect upon biological targets. By blocking the S1 protein binding to its ACE2 receptor, flavonoids in  Carica papaya  may contribute to its prevention toward the SARS-Coronavirus cell invasion. The immunomodulatory effects of  Carica papaya  leaves may also work on the target, therefore, this is in agreement with the study by Norahmad et al. The second hypothesis on how the  Carica papaya  leaves would eradicate the SARS-Coronavirus-2 cell proliferation is by blocking the 3CLpro enzymatic activity. This is proven by the therapeutic effect of flavonoid quercetin to inhibit DENV2 NS2B-NS3 protease which is also a serine protease enzyme as studied by de Sousa et al [60], [61].

The  in silico  structure-based drug design has been applied in this study to speed up the discovery of SARS-Coronavirus-2 antiviral agents using an inexpensive and rational approach [62], [63].  Carica papaya  is a natural product which is a highly available and abundant resource that could lead for an effective and efficient herbal drug [64] to combat SARS-Coronavirus-2 pandemic.

Conclusion

In conclusion, Carica papaya leaves have the potential to be the SARS-Coronavirus-2 antiviral agent from herbal. This is due to the 20 compounds presenting in its leaves which have drug-like likeness structure, non-carcinogenic, non-toxic, pharmacokinetically and pharmacodynamically stable as predicted by in silico experiments. The major flavonoid compounds that are interacting with a diverse protein target of SARS-Coronavirus-2, could be the rational reason on how this herbal would be promising as this antiviral agent.

CRediT authorship contribution statement

Pandu Hariyono: Investigation, Conceptualization. Christine Patramurti: Writing - review & editing. Damiana S. Candrasari: Writing - review & editing. Maywan Hariono: Conceptualization, Writing - review & editing.

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.