Theoretical study on inhibitability of some natural alkaloids against influenza virus hemagglutinin and SARS‐CoV‐2 main protease

Abstract Berberine (V1), lycorine (V2), hemanthamine (V3), aloperin (V4), dendrobine (V5) possess structural frameworks resembling known anti‐influenza and anti‐SARS‐CoV‐2 drugs, thus subjected for a computational screening. Their quantum properties were examined using density functional theory (DFT); the ligand‐protein inhibitability was evaluated using molecular docking simulation; physicochemical properties were obtained from QSARIS‐based analysis in reference to Lipinski's rule of five; pharmacokinetic parameters were assessed by ADMET‐based analysis. DFT calculations indicate that there are no abnormal bonding constraints observed; NBO analysis suggests all possessing favorable electric configurations for intermolecular inhibition. Regarding ligand‐2VIU, the order for static inhibitability is V3‐2VIU > V2‐2VIU > V1‐2VIU > V5‐2VIU > V4‐2VIU; Regarding ligand‐6LU7, the corresponding order follows: V2‐6LU7 > V3‐6LU7 > V1‐6LU7 > V5‐6LU7 > V4‐6LU7. An exceptional hydrophilic bonding (π‐cation) with the associated Gibbs free energy of ‐10.9 kcal.mol‐1 is detected in inhibitory complex V1‐2VIU. QSARIS‐based analysis reveals that all the candidates are highly bio‐compatible. ADMET‐based analysis specifies V2 and V3 as being safe and suitable for the use as orally administrated drugs. The results encourage further investigations for more in‐depth mechanisms and experimental validations, such as molecular dynamics simulation and in vitro enzyme assays.


INTRODUCTION
Seasonal flu or influenza is a serious public health concern causing not only severe illness even deaths but also substantial economic burden also. [1] Annually, this acute respiratory infection is estimated to cause more than 3 million severe cases, and up to 650,000 deaths worldwide. In Vietnam, this infectious disease is responsible for more than 16,000 deaths in 2018 according to latest reports by the World Health Organization. Influenza is an enveloped, negative-sense single-stranded RNA virus, which belongs to the Orthomyxoviridae family. [2] They are classified and named by type: A, B, C and D of which A, B and C cause seasonal flu in humans [3] ; while, newly-characterized type D is currently reported to infect exclusively cattle, but its zoonotic potential remains questionable. [4] Type A viruses are of highest clinical interest for being responsible for 4 influenza pandemics over the last century. [5] This highly infectious type is further categorized into different subtypes based on the properties of two surface proteins, hemagglutinin (HA) and neuraminidase (NA). Up-to-date studies document 18 different HA subtypes (H1-H18) and 11 different NA subtypes (N1-N11), theoretically resulting in 198 potential subtype combinations. [6,7] While only 131 subtypes have been identified in nature, the risk of reassortment between HA and NA genes persists, potentially leading to novel antigen combinations. And a novel influenza subtype eventually causes a pandemic. [8] Although vaccines still remain the core strategy for influenza prevention, they are confronted with a number of challenges. First of all, annual reformulation and yearly administration are always needed due to the constant antigenic drift and shift of seasonal influenza viruses. [9] The rapid "mutation escapes" possibly lead to ineffectuality of the vaccine protection, posing rapid emergence of pandemic threats. [10] In addition, the predominant use of embryonated chicken eggs for influenza vaccine manufacture greatly curbs their worldwide production, regarding both timeline and scale. [11] COVID-19 pandemic also worsens the supply shortage of other vaccines. Nevertheless, recent data suggests that influenza vaccination potentially provides addition protection in SARS-CoV-2 positive patients. [12] Scientists have observed constant mutation of SARS-CoV-2 and rapid rises of its variants. The severity of COVID-19 pandemic is likely associated with the genetic variance. It was observed that among the different variants, delta variant is highly contagious and likely leads to more severe cases. US Food and drug administration (FDA) has authorised many vaccines that are proven remedies against delta and other known variants, but it is not sure whether these vaccines will remain effective against coming variants. By mid-2021, national regulatory agencies (NRA) have approved twelve vaccines for the public, e.g. Pfizer-BioNTech, [13] Moderna, [14] Sputnik V, [15] Oxford-AstraZeneca, J&J. [16] According to June 2021 report of World health organization (WHO), vaccines like AstraZeneca/Oxford, J&J, Moderna, Pfizer/BioNTech, and Sinopharm Vero Cell met the necessary criteria for safety and efficacy. [17] Most of these vaccines are currently in use in different countries.
Other than vaccination, scientists worldwide are trying to find a medicine (i.e. anti-viral drugs or inhibitors) that can complete treatment for viral infection in general, and for influenza and SARS-CoV-2 in particular. Besides the cost-effective immunoprophylaxis, antiviral agents also play a pivotal role in the anti-influenza therapy, FDA has been approved only 6 drugs for the chemoprophylaxis and treatment of influenza, including M2 ion channel protein inhibitors (amantadine and rimantadine), NA inhibitors (oseltamivir, zanamivir, peramivir) and cap-dependent endonuclease inhibitors (baloxavir). [18] However, type A viruses have developed multiple-level resistance against current medications. The first-generation anti-influenza agents, amantadine and rimantadine, are actually no longer recommended due to widespread resistance. [19] Resistant strains against newer agents have also been reported recently. [20] In consequence, continuous discovery and development of novel anti-influenza are urgently required. In terms of SARS-CoV-2, the unavailability of a proper drug against the virus and the emergence of its variants created a global crisis. Some commonly repurposed drugs and their combination used to treat COVID-19 infection are: Remdesivir, Favipiravir, Lopinavir, Ivermectin, Doxycycline, Nitazoxanide, etc. [21] However, these potential drugs and their combination have possible side effects. This can be considered as the reason for a positive impact on the herbal medicine market globally as herbal extracts are the best ways to boost up immunity in the outbreaks of COVID-19. The fact is that people are turning towards herbal medicines and herbal medicine therapy, which can also be supplementary products for the treatment. Chemists are exploring the natural possibilities of herbal medicines for the treatment of COVID-19. Many medicines comprising herbal extracts like turmeric, ginger, Triphala ® (a mixture of dried fruits), and garlic are taken in the medicinal form to assist the immune system. [22][23][24][25] Structurally diverse chemical scaffolds from natural products could serve as screening libraries for lead compound identification in drug discovery. [26] While oseltamivir, zanamivir and peramivir are synthetic compounds, their research and development process is initiated by Mother Nature with the starting materials are of natural origin. [27] Additionally, natural-based agents have long been used for the treatment of acute respiratory infection, especially in children, thanks to medication unrestrictedness, costeffectiveness and safeness. [28] Among active components of medicinal herbs, alkaloids and flavonoids are of intense interest because of their high diversity, widespread distribution, and broadness of biological activity, especially including antiinfluenza property. [29,30] Compared with laboratory experiments, molecular data-based approach provides the possibility to identify compounds with potent inhibitory activity, and to design novel active agents in a rapid and high-throughput manner. Classic experimental procedures in-turn help to validate the results predicted through such computational methods. [31] This study investigates the binding affinity of five alkaloid derivatives, including berberine (V1; PubChem CID), lycorine (V2; PubChem CID), hemanthamine (V3; PubChem CID), aloperin (V4; PubChem CID), dendrobine (V5) (figure 1), towards influenza A/X-31(H3N2) virus hemagglutinin 2VIU (DOI: 10.2210/pdb2VIU/pdb) and SARS-CoV-2 main protease 6LU7 (DOI: 10.2210/pdb6LU7/pdb) was investigated. The structural formula of the ligands and their targeted proteins structures are presented in figures 1 and 2, respectively. The compounds, in general, possess the molecular framework similar to those of well-known antimicrobial medicines. Their quantum properties were determined by DFT calculations. The evaluation for their inhibitability towards the targeted protein structures is based on molecular docking simulation. The bio-compatibility of the ligands is justified by physicochemical properties obtained from QSARISbased analysis in reference to Lipinski's rule of five. The applicability of medicinal development is predicted by ADMET-based pharmacokinetics and pharmacology properties.

Quantum chemical calculation
Molecular quantum properties of the natural alkaloids (V1-V5) and their optimised geometry were examined by density functional theory (DFT) using Gaussian 09. The computation was under no symmetry constraints [32] and at level of theory M052X/6-311++G(d,p). [33] Structural global minimum on the potential energy surface (PES) was confirmed by calculation of vibrational frequencies on the respective molecules. The basis set def2-TZVPP [34] was used for frozen-core approximation of non-valence-shell electrons, which converged to optimised geometries with single-point energies at the level M052X/6-311++G(d,p). Each optimisation run was under resolution-of-identity (RI) approximation. NBO 5.1 available was utilised for frontier orbital analysis at level of theory M052X/def2-TZVPP. [35] The highest occupied molecular orbital (HOMO) energy, EHOMO, can be interpreted as the tendency of intermolecular electron donation; on the other hand, the lowest unoccupied molecular orbital (LUMO) energy, EHOMO, represents electron-accepting ability. By exhibiting molecular electronic excitablity, energy gap E = ELUMO -EHOMO can suggest intermolecular reactivity of the host molecule. The frontier molecular orbital is simulated for the molecule using Koopman's theorem which was used for calculation of ionisation potential (I) and electron affinity (A) expressed as: I = -EHOMO and A = -ELUMO. Parameters I and A are useful in analyzing the electronegativity (χ), hardness (η), softness (S), and other frontier molecular orbital parameters, which are significant in assessing the reactivity of a molecule. The molecule electronegativity () is further inferred via equation:  = (I + A)/2 while the hardness (η) is calculated as: η = (I -A)/2 and the softness S = 1/η.
-Determination of bonding position: The important amino acids were defined as within a ligand-protein radius of 4.5 Å. Redundant residues, e.g. water and small molecules, were removed before reestablishing the enzyme action areas reestablishing the enzyme action areas.
-Lowest energy state of ligands: The 3D molecular structure of the compound was optimised following the purpose to correct the values of bond lengths, bond angles, bending angles, and unusual non-bonding interactions. This correct the coordinates of the molecular atoms in different parts of the compounds.
-Configuration of protein structures: Configuration for preparation of the protein protonation structures: Tether -Receptor with the strength of 5000; Refine of 0.0001 kcal•mol -1 •Å -1 . The achieved structures were saved in *.pdb format. Configuration for preparation of the ligand structures: Conj Grad for minima energy; termination for energy change = 0.0001 kcal‧mol -1 ; max interactions = 1000; modify charge: Gasteiger-Huckel. The configuration with minimum steric energy was determined from five automatic iterations.
Step 2: Docking investigation -Configuration for intermolecular interaction simulation: poses retaining for intermolecular interaction probing = 10; maximum solutions per iteration = 1000; maximum solutions per fragmentation = 200. The obtained data for inhibitory structure of the ligand-protein complexes was saved in format *.sdf.
Step 3: Re-docking -Re-docking of protein-compound co-crystal structure: the re-docking of protein-ligand complex co-crystal structure aimed to evaluate the reliability of docking parameters. The process was carried out with three structures of compounds as follows: (1) Separation of compounds from homogenised complexes in protein.
(2) Separation of compounds from homogeneous complexes and re-preparation.
(3) Preparation of the investigated compounds (structure drawing, structural optimisation parameters on minima energy).
Step 4: Post-docking analysis -Docking score (DS) energy indicates the Gibbs free energy of ligand-protein inhibitory system, which is considered as the primary indicator of the duo-system inhibitability. Intermolecular interactions formed between the ligand and in-pose amino acids of the protein include hydrophilic binding, electrontransferring (H-acceptor/donor), cation-arene (H-π), arene-arene (π-π), and ionic and hydrophobic interaction, van der Waals forces. The simulation provides data on in-bonding amino acids, bonding lengths, and their Gibbs free energy in regard of these interactions. Static conformation of an inhibitory complex can be assigned to root-mean-square deviation (RMSD) as it represents the average between neighbouring atoms. In-pose arrangement of the ligands was rendered on 2D and 3D planes.

QSARIS-based analysis
Drug-likeness properties of phytochemical are screened in order to verify the pharmacokinetic ability of the drug like molecule. Thus, physicochemical properties were obtained by QSARIS system via Gasteiger-Marsili method. [39] The parameters include molecular mass (Da), polarizability (Å 3 ) and volume or size (Å), and dispersion coefficients (logP and logS), subjected for assessment of orally pharmacological compatibility based on Lipinski's rule of five. [40] According to the criteria, a well membrane-permeable candidate is predicted following the requirements: (1) Molecular mass < 500 Da; (2) no more than 5 groups for hydrogen bonds; (3) no more than 10 groups receiving hydrogen bonds; (4) the value of logP is less than +5 (logP < 5). [41,42]

ADMET-based analysis
ADMET properties are analysed as a vital step for drug discovery. A web environment developed and maintained by the Molecular Modeling Group, Swiss Institute of Bioinformatics, i.e. SwissADME (http://www.swissadme.ch/), was utilised to yield pharmacokinetic parameters including chemical absorption, distribution, metabolism, elimination and toxicity (ADMET) of the studied compounds. Along with pharmacokinetic rules, the inhibitor must also follow the ADMET properties. Pires et al. proposed the theoretical interpretations for the parameters, [43] which are powered by the University of Melbourne and University of Cambridge for public reference (http://biosig.unimelb.edu.au/pkcsm/theory). Figure 3 shows the structures of V1-V5 optimised C1 symmetry in the ground state using level of theory M052X/6-311++g(d,p); The ground state electronic energy and dipole moment values of the structures are presented in table 1. Overall, the optimisation reaches a convergence for all the compounds, validating their natural stability (i.e. existing in natural sources). There are no abnormal bonding constraints, either angles or lengths. All the calculated figures for the latter are in the characteristic values, e.g. 1.54 Å for C-C, 1.09 Å for C-H, 1.43 Å for C-O (oxatriquinane), 1.23 Å for C=O (carbonyl compounds), ~1.5-2.5 Å for N-H, and 1.47 Å for N-H (imine). The ground state electronic energies of the optimised structures vary from -708177.49 (V1) to -436517.10 (V4) kcal.mol -1 . The corresponding figures for the dipole moment are from 1.010 (V4) to 4.332 (V5) Debye. In essence, higher dipole moment values are conducive to the formation of a bond or complex between the ligand and the target protein. [44,45] The highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of the studied compounds V1-V5 are shown in Figure 4, and in-detail quantum chemical parameters are given in Table 2. Regarding spatial localisation, the electron densities of V1 are distributed rather evenly over its molecular plane, cf. those of other compounds. This preliminarily implies that the structure might be more flexible in terms of charge-transferring, which are highly conducive to intermolecular interaction. The argument is due to theoretical charge-transferring properties of HOMO (representing intermolecular electron donation tendency) and LUMO (representing electron-accepting ability. However, this does not guarantee a predominant stability as the property is also related to chemical/physical affinity of the ligands and their targeted structures. Regarding parameterisation, all structures (V1-V5) are considered exhibiting electronic stability given their significantly low EHOMO (commonly agreed under -5 eV), especially that of V1 (-9.028 eV). In a preceding work, theoretical complexes of a tetrylone family with the corresponding figures ranging from -3 to -7 eV were already assessed as highly stable. [46] Besides, the band-gap energy (EGAP) varies between 5.066 and 6.833 eV, lying on the transition of an insulator (> 9 eV) and a semiconductor (< 3.2 eV). [47] The property of conductivity is thought to signify the intermolecular binding capability towards protein structures since the conductivity of the polypeptide molecules was well-observed and explained based on super-exchange theory (or electron tunneling) [48,49] and the electron hopping model. [50] On the other hand, the corresponding ionisation potentials (aka. the opposite of EHOMO) of significance signify the difficulty of chemical ionisation, thus suggesting the tendency of intermolecular binding (via hydrophilic and hydrophobic bonds) over compositional changing (by chemical reactions). Therefore, given justification based on quantum chemical properties, all the studied are proposed as promising for protein inhibitability in general, and towards structures of SARS-CoV-2 main protease (6LU7) and influenza virus hemagglutinin (2VIU) in particular.    The complexes regarding these inhibitions are selected for more in-depth discussion on their inhibitability.   In another works, the inhibitory effectiveness of semi-synthesised derivatives from Dolichandrone spathacea iridoids and Dipterocarpus alatus dipterocarpol towards diabetes-related proteins were experimentally justified of profound significance by DS values at ca. -13 to -15 kcal.mol -1 (corresponding to IC50 values < 10 μM; control drug IC50 ca. 300 μM). [51,52] This indicates the promising potentiality of V3 and V2 as their corresponding figures are not in substantial inferiority; also, they are expected to form rigid conformations with the protein (except for V1-2VIU), which are widely accepted if < 2. [53] Otherwise, the docking-based investigation on V1-2VIU inhibitability finds an exceptional hydrophilic bonding (π-cation) with the associated Gibbs free energy of -10.9 kcal.mol -1 , between an aromatic cyclohexene of the ligand and a nitrogen atom of Lys A326. The value of significance is firstseen among all our resembling researches. This means although V1-2VIU does not register a highly stable and conformation-rigid structure in overall, it is still possible that V1 could channel sufficient distortion forces to the secondary-tertiary structure of the protein via the bonding, thus inducing denaturation and loss of shape-based activity. However, whether the amino acid (aka. Lys A326) contributes significantly to the deformity of the protein structure is still an open contemplation, which can be answered by more powerful implementations such as molecular dynamics simulation, or experimental trials such as in vitro enzymatic assays. If the hypothesis is verified, a very promising inhibitor towards influenza virus hemagglutinin (2VIU) can be introduced for further development.
Regarding 3D in-pose morphology and 2D interaction map of all ligand-protein inhibitory complexes, the sites are observable rather closed and small. The former is likely conducive the confinement on the ligands yet limit their entry; while, the latter indicates that the unavailableness for inhibition of multiple or larger ligands. Besides, all the ligands are still expected to be size-fitting (from 3D observation) and shape-complementary (by 2D continuous proximity contours); otherwise, V1-2VIU is the exception. DS: Docking score energy (kcal.mol -1 ); RMSD: Root-mean-square deviation (Å) L: Ligand; P: Protein; T: Type; D: Distance (Å); E: Energy (kcal.mol -1 ). Table 6 summarises QSARIS-based physicochemical properties of the investigated compounds (V1-V5), which are molecular mass (amu), polarisability (Å 3 ) and size (Å) as well as the logP and logS dispersion coefficients. The parameters are in consideration as an attempt to evaluate their pre-docking biocompatibility, i.e. the interactions between the ligands and potential plasmatic components in the polarised media of biological bodies. First, none of the candidates possess molecular mass over 500 amu or dispersion coefficient over +5. According to

Physicochemical properties
Lipinski's rule of five, they are assessed with good biocompatibility (in case of development for orally administrated drugs). Besides, their hydrophilic behaviour referenced from molecular docking simulation appears to satisfy Lipinski's requirements in general. Other than that, their polarisability constants are of significance, i.e. > 28.3 Å 3 , likely inducing the formation of molecular dielectric moments. [56] Furthermore, it is noteworthy that a framework based on docking-QSARIS combination was proposed as a good model to predict the inhibitory behaviour of ligand-protein biological inhibitory systems, in general. [51] Table 7 presents chemical absorption, distribution, metabolism, excretion, and toxicity (ADMET) of the investigated compounds (V1-V5). The theoretical interpretations of the parameters were described by Pires et al. [43] and powered by The University of Melbourne and University of Cambridge for public reference (http://biosig.unimelb.edu.au/pkcsm/theory). First, all the promising inhibitors (V2-V3 and V1) are predicted to be absorbed effectively through intestine by the percentage of over 60 % (considerably high at > 30 %), they perform moderate Caco-2 permeability with the corresponding values of logPapp from 0.107×10 -6 to 1.844×10 -6 cm.s -1 (considerably high at > 8×10 -6 cm.s -1 ). The latter means their overall intestinal processes, i.e. transcellular transport, paracellular transport, and some aspects of efflux-active transport, would meet certain resistance. However, this can be overcome by transport-enhanced synergic agents for more inpractice applications. The modelling also predicts no ligand-P-glycoprotein interactions of V3, thus no inhibitory effects on transmembrane activities. Second, all the test compounds seem not to accumulate into specific parts of the body. In fact, their VDss (in log L.kg -1 ) register within 0.687 and 0.882, implying balanced tissue-plasma distributions. In particular, V2 and V3 are considered to be unable to cross the blood-brain barrier given their low BBB permeability (logBB < 0.3), thus poorly distributed to the brain; also, their CNS permeability indicates their low penetrability to the central nervous system. Third, the studied compounds register different effects towards the inhibitors and substrates relating to the cytochromes P450 family. While V3 expects to interact with CYP3A4 substrate and CYP2D6 inhibitor, V2 is predicted to interact with only the former. Otherwise, no data on the quantitative structure-activity relationship suggest serious intervention of these compounds with other structures of the cytochromes. Therefore, in general the compound can be effectively oxidised in the liver. Overall, if the compounds are in consideration for drug development, the monooxygenase-based potential interaction leads to a certain caution for the use on patients with liver-related prehistoric conditions. Regarding excretion, except for V2, all other candidates can be carried out of the body by Organic Cation Transporter 2; also, their total clearance is ca. 1 log ml.min -1 .kg -1 .

ADMET-based pharmacokinetics and pharmacology
In terms of toxicity, all the tested compounds are safe. None of them expresses apparent AMES toxicity, indicating no mutagenic potential; most do not interact with hERG inhibitors, indicating no fatal ventricular arrhythmia. Besides, hepatotoxicity and dermal sensitisation are exclusive. However, the regression model predicts all the studied compounds are considered toxic towards T. Pyriformis (pIGC50: from 0.343 to 0.679 log μg.L -1 ) yet generally safe to Flathead Minnows (log LC50: from -0.586 to 1.959 log mM), yet, the values are not significantly over the standard threshold provided. The negative logarithm of the concentration required to inhibit 50 % growth of the former is predicted as toxic if > -0.5 log μg.L -1 ; meanwhile, the lethal concentration values for the latter are assessed as toxic if < -0.3. Since the protozoa bacteria is often used as a toxic endpoint, this abnormal regressive result requires special concern from further investigations or development. Apparently, the in silico implementation supports the application of the docking-based most promising inhibitors (aka. V2 and V3) into further investigations for drug development with the maximum tolerated dose (in log mg.kg -1 .day -1 ) estimated to be -0.346 and -0.608, respectively; yet, appropriate synergic agents are still needed to reconcile the resistance given by the overall intestinal transport and sufficient dosage should be determined carefully from experimental and clinical trials.  . The findings include an exceptional hydrophilic bonding (π-cation) with the associated Gibbs free energy of -10.9 kcal.mol -1 in inhibitory complex V1-2VIU. QSARIS-based analysis reveals that all the candidates are highly biocompatible. ADMET-based analysis especially suggests V2 and V3 as being safe and suitable for the use as orally administrated drugs. The results encourage the utilisation of more powerful implementations, such as molecular dynamics simulation, to clarify the possible inhibitory effects of V1 (berberine) against the enzymatic function of 2VIU (against influenza virus hemagglutinin), and further investigation for drug development.