key: cord-0055279-dunvc8se
authors: de Sá, Ézio R. A.; Costa, Allan N.; Costa, Rayla K. M.; Souza, Janilson L.; Ramos, Ricardo M.; Lima, Francisco das C. A.
title: In silico study of the interactions of Pilocarpus microphyllus imidazolic alkaloids with the main protease (M(pro)) of SARS-CoV-2
date: 2021-01-12
journal: nan
DOI: 10.1080/08927022.2021.1873321
sha: 021032984c089a8fb888a2f6ecf43e070563ef95
doc_id: 55279
cord_uid: dunvc8se

The disease outbreak caused by SARS-CoV-2 continues to rise worldwide, even in countries which have considered it controlled. As new cases appear daily, infecting millions of people and causing thousands of deaths, the current in silico study aims to investigate the imidazolic alkaloids of the species Pilocarpus microphyllus (Jaborandi) as a potential inhibitory activity against the M(pro) protease from SARS-CoV-2, since it plays a fundamental role in the processing of polyproteins that are translated from viral RNA. Jaborandi is distributed in some Brazilian biomes, being easily identified, yet little researched, with proven anti-inflammatory, contraceptive, anti-diabetic and gastroprotective activities. In this work, DFT calculation of thermodynamic properties, electrostatic potential surface, frontier molecular orbitals and descriptors of chemical reactivity of imidazolic alkaloids were associated with the use of molecular docking techniques, molecular dynamics and ADMET predictions. One can verify a good reactivity chemistry and energetic stability of epiisopiloturine, epiisopilosine, isopilosine and e pilosine with some residues of amino acids present in the active site of the main protease of COVID-19. In this sense, the results point out to the imidazolic alkaloids of Jaborandi as promising targets for in vitro and in vivo tests, as possible candidates for inhibitors of the enzyme M(pro).

The disease known as COVID-19 is of zoonotic origin and has caused the severe acute respiratory syndrome by coronavirus-2 (SARS-CoV-2) in human beings [1] . In December 2019, the disease was diagnosed for the first time in patients from Wuhan, a city with about 11 million inhabitants, in the capital of the province of Hubei, located in the People's Republic of China [2] . Due to the high rates of the new coronavirus transmission among human beings, on March 11, 2020, the World Health Organization (WHO), classified its contagion as a worldwide pandemic [3] .

Hence, SARS-CoV-2 is highly contagious, and easily transmitted through respiratory droplets released into the air by infected people through their coughing or sneezing [1] . According to the Pan American Health Organization, PAHO/WHO Brazil (https://www.paho.org/bra/), by November 30 of 2020, there are a total of 62.363.527 cases (496.892 new cases compared to the previous day) and 1.456 .687 deaths (7. 697 new deaths in relation to the previous day), resulting from COVID-19 [4] .

After investigating the coronavirus organism, the researchers found that its main SARS-CoV-2 protease (M pro ) is an essential enzyme for processing the polyproteins translated from viral RNA. Thus, by blocking the process in one of the main steps (the cleavage of precursor polyproteins), the vital cycle of viral reproduction would be impaired [5] ; indicating it to be a promising target for the development of drugs to combat SARS-CoV-2 [6, 7] .

The crystallographic structure of M pro complexed to the N3 ligand was recently deposited in the Protein Data Bank (PDB) with a resolution of 2.16 Å [8] . Since N3 is a Michael acceptor inhibitor, developed by a computer aided drug design, it is able to specifically inhibit M pro from multiple coronaviruses, including SARS-CoV and MERS-CoV [9, 10] . It is covalently linked to the receptor and is being investigated in silico and in vitro to discover the mechanism of action against this protease [8] .

It is well known that natural herbal medicines are important sources of active ingredients for the treatment of diseases. Affected by the traditional medicine, which effectiveness comes from medicinal plants use, the modern medicine also uses many of these effective healing agents [11] . In this scenario, the researches involving extracts of medicinal plants are very relevant, as there is an enormous biodiversity in Brazil, which enables the large-scale investigation of substances with pharmacological activities not unknown by science [12] .

For instance, it is possible to mention the use of the plant species Pilocarpus microphyllus Stapf ex Wardleworth, belonging to the Rutaceae family, originally from the North and Northeast regions of Brazil, popularly known as Jaborandi, with important medicinal properties; but with applications still little known [13] .

The extract of P. microphyllus is indicated for the treatment of the flu, fever, inflammation, pneumonia, asthma, glaucoma, diabetes, rheumatism, and others [14] . It also has antidiabetic activities [15] , and among its constituents which have already been investigated are: pilocarpine, in the use of glaucoma and xerostomia treatment [16, 17] ; epiisopiloturine, in the fight against schistosomiasis [18] [19] [20] ; as well as anti-inflammatory, contraceptive [21] and gastroprotective activities [22] ; and the molecules of epiisopilosine and isopilosine with pharmacological action against S. mansoni [23] . Under those circumstances, the phytomolecules of Jaborandi become promising targets for the study of the properties not known or little described in the literature.

Given the importance and the urge of the discovery, in conjunction with the planning of new drugs, molecular modelling emerges as a powerful technique, which provides a set of tools to be used in the processes of identification, selection, manipulation, optimisation and characterisation of new drug candidates, due to its high capacity for forecasting and obtaining the thermodynamic, energetic and structural properties of the complexes (receptor-ligand) analysed [24] .

Accordingly, the molecular docking method has stood out as an important fitting tool for predicting interactions between macromolecular ligands and receptors [25] . In addition to it, the molecular dynamics has also shown expressive importance through the investigation of variations in the movement of ligands over time, the mechanism of diffusion, the folding of molecular chains and structural changes caused by interactions at receptor binding sites [24] . Regarding the pharmacokinetic properties in silico, it has been an investigative approach widely used in the initial study of the properties of Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET), reducing costs with biological tests of unattractive substances, and providing greater agility in drug selection [26] .

Therefore, as there are still no drugs available and effective treatments for the cure of COVID-19, the present in silico study aims to investigate the reactivity and chemical interaction of imidazolic alkaloids of the species Pilocarpus microphyllus (Jaborandi) with the M pro protease of SARS-CoV-2.

The crystal structure of the SARS-CoV-2 M pro receptor (PDB code: 6LU7) was obtained from the Protein Data Bank [8, 27] . The three-dimensional structures of the ligands Pilosine (PS), Isopilosine (IPS), Epiisopilosine (EPS), Epiisopiloturine (EPR), Pilocarpine (PC), Isopilocarpine (IPC), Pilocarpidine (PD), Isopilocarpidine (IPD), Pilosinine (PN) e 13-nor-7 (11) -dehydro-pilocarpine (13N), (11) -dehydro-pilocarpine (13N), found in the extract of Jaborandi leaves [23, [28] [29] [30] , were built with the GaussView software version 6 [31] .

In order to obtain the global minimum energies, thermodynamic properties, Molecular Electrostatic Potential Surface (MEPS), frontier molecular orbitals and descriptors of the global chemical reactivity of the ligands, calculations geometric optimisation and frequency in the gas phase (vacuum system) were performed, at the level of Density Functional Theory (DFT) and based on studies by Rocha et al. [20] and Costa et al. [32, 33] . For this purpose, it was used the Gaussian 09 software [34] and the B3LYP hybrid method [35, 36] associated with the basis set 6-311++G(d,p) [37, 38] . The descriptors of global chemical reactivity were determined based on the theorems of Janak [39] , Perdew et al. [40] and Parr et al. [41, 42] .

The preparation of the ligands and the receptor, just as the visualisation of the complexes, were carried out with the UCSF Chimera software [43] . During the protease preparation phase, water molecules and the N3 crystallographic ligand were removed. The addition of hydrogen, the calculation of the Gasteiger loads of the receptor, and ligands were performed using the AutoDock Tools (ADT) software, version 1.5.6 [44, 45] . Molecular docking calculations were performed with AutoDock Vina [46] software. During the process, the receptor was considered as a rigid structure and the ligands as flexible ones [47] . The Met165 residue was chosen because it is part of the active M pro protease site and performs interactions with the side chains of the N3 ligand [8] , being your cartesian coordinates (x = −16.406, y = 13.719 and z = 67.212) used to define the centre of the grid box [27] .

In the Autodock, the Lamarckian genetic algorithm was used as a research parameter [48] . The size of the cubic box generated by AutoDocK Vina, in the region of the receptor interaction, was defined as 22.5 Å for each cartesian axis, the number of modes was set to 50, and the exhaustiveness to 24 [49] . The analysis of molecular docking was concentrated in the lower energy clusters, with the most stable conformation being chosen. The binding energy of the N3 critalographic ligand was calculated by molecular docking as an observation model for the imidazolic alkaloids of Jaborandi. Intermolecular interactions (hydrogen bonds, electrostatic, van der Waals and hydrophobic contacts) were identified and visualised by the Discovery Studio Visualizer 2020 software [50] .

The initial coordinates for the molecular dynamics simulations were obtained from the complexes with the lowest binding energy, determined by molecular docking. The CHELPG method (CHarges from ELectrostatic Potentials using a Grid-based method) was applied to obtain the atomic charges of the ligands [51] , based on the data presented in the literature [52, 53] , and from the geometric optimisation calculations with the B3LYP/6-311++G(d,p) model and use the Gaussian 09 software [34] .

For the protonation of the SARS-CoV-2 M pro enzyme, the H++ online server was used [54] . The molecular dynamics simulations were performed with the GROMOS96 force field 53a6 included in the GROMACS package, version 2018.1 [55, 56] . The systems were simulated using the NPT set (in which the number of particles, pressure, and temperature are constant), and periodic boundary conditions (cubic) [57] . The water molecules during the process were made explicit by the Single Point Charge (SPC) model [58] . The simulation time was 50 nanoseconds (ns), with an integration step of 2 fs [59] . The systems were gradually warmed up, starting with 100 K (10 ps), 150 K (5 ps), 200 K (5 ps), 250 K (5 ps), and later being adjusted to 310 K [57] . The initial 6 ns of each simulation were considered as part of the heating (0.025 ns) and equilibrium (5.975 ns) steps, not being used in the data analysis [57] .

The temperatures of solvent and solutes (protein, ligands, water, and sodium ions) were independently coupled to a thermal bath with a relaxation time of 0.1 ps using the V-rescale thermostat [57] . The pressure in the system was weakly coupled to a pressure bath with 2 ps of relaxation time, using the Parrinello-Rahman barostat [60, 61] . Bond lengths were constrained using the LINCS algorithm [62] , and electrostatic interactions among non-ligand atoms were evaluated by the Particle Mesh Ewald (PME) method [63] . GROMACS provides the ability to calculate the energies of short-range (Lennard-Jones potential) and long-range (Coulomb potential) unbound interactions and their error estimates for the receptor-ligand complexes [64, 65] , being the sum of the energies realised posteriorly. The visualisation and graphic production of the results happened throughout the use of the UCSF Chimera software [43] .

The obtaining of the physicochemical and pharmacokinetic properties of the ligands occurred through the in silico search, using the structures (2D, 3D or SMILES) on the online platforms: PreADMET [66] , FAF-Drugs4 [67] , SwissADME [68] and PASS Oline [69] .

Some of the parameters analysed were: Human Intestinal Absorption (HIA), Penetration in the Blood-brain Barrier (BBB), Lipinski's rule, permeability in Caco2 and MDCK cells (Madin-Darby canine kidney), aqueous solubility, the mutagenicity of species Ames Salmonella (TA100 and TA1535), the carcinogenicity, mutagenicity and bioactivity assays [70] .

The 3D structures of the molecules ( Figure 1 ): pilosine, isopilosine, epiisopilosine, epiisopiloturine, pilocarpine, isopilocarpine, pilocarpidine, isopilocarpidine, pilosinine and 13-nor-7 (11)-dehydro-pilocarpine, are the most stable for the model B3LYP/6-311++ G(d,p).

In this study, it is possible to detect the existence of four isomeric forms (C 16 [20] , observed through the rotation of two conforming atoms (C4 and C9), and by the presence of three chiral carbons (C5, C7, and C8), which corroborate with the data presented in the experimental studies of Rocha et al. [23] .

It also points out Rocha et al. [20] that the chemical structure of these alkaloid isomers has a difference in stability between their dihedral angles. Because, when acquiring a cislike geometric conformation, the PS alkaloid ( Figure 1d ) tends to be less stable than the EPR (Figure 1a ), demonstrating to be the most stable one due to its trans-type geometric configuration. Then, the frequency calculations were performed to determine the thermodynamic properties of imidazolic alkaloids, as shown in Table 1 .

It can be seen in Table 1 , the proposed model describes the energetic properties of the systems well, demonstrating that these data correspond to the thermodynamically more stable structures for the respective ligands ( Figure 1 ). However, those with the lowest energy values of the parameters: zero-point, thermal energy, enthalpy, and Gibbs free energy, are considered the most stable ones. Thus, the best results are found in the compounds EPR, EPS, IPS, and PS respectively, indicating a great similarity between their energy parameters, possibly caused by the phenomenon of isomerism. Meanwhile this, the least stable is the PN molecule, due to its associated value, caused by the different geometric properties. However, it is necessary to understand that properties, such as enthalpy and entropy, can undergo energetic changes, according to the contributions caused by the vibrational, translational, rotational, and electronic movements of the molecules [71] .

MEPS is a map that represents the distribution of electronic density on the surface of a molecule. This one is widely applied in predicting sites and relative reactivity for electrophilic attack, in studies of biological recognition and interactions by hydrogen bonds. Figure 2 shows the imidazolic alkaloids, in which the negative electrostatic potential (in red) represents the attraction of the proton by the region of high concentration of electrons in the molecule, whereas the positive electrostatic potential (in blue) corresponds to the repulsion of the proton by atomic nuclei [32] .

It is possible to observe, in Figure 2 , the electronic density for each molecule is represented by red regions (polar and negatively charged) and blue regions (non-polar and positively charged). A greater electron density is present around the nitrogen atom (N3) of the imidazolic ring and the oxygen atoms (O1 and O3) of the γ-butyrolactone ring (oxolan-2one), thus representing the most likely locations for the occurrence of chemical interactions.

The energy of the frontier molecular orbitals, HOMO (Highest Occupied Molecular Orbital) and LUMO (Lowest Unoccupied Molecular Orbital), provides information about the donor and acceptor electron character, respectively, of a compound [72] . The HOMO and LUMO are mainly responsible for the biological interactions between ligands and proteins in a complex [73] [74] [75] .

The energy values of the HOMO and LUMO orbitals are fundamental to determine the electrical properties, kinetic stability, and the descriptors of global chemical reactivity in a molecule [33] . Thus, Table 2 contains some electronic properties determined from the calculation method B3LYP/6-311+ + G(d,p) used for the imidazolic alkaloids of the species P. microphyllus.

It is observed that the compounds, in Table 2 , demonstrate a similarity in the electronic properties, however, the substance EPR stands out in relation to its greater capacity for charge transfer based on its chemical potential (−3.57 eV). The EPS, IPS and PS molecules have a better donor electron character (nucleophilic), and acceptor electron character (electrophilic) in relation to the others molecules. With based in the HOMO and LUMO energies, this molecules show a good chemical reactivity. The 13N molecule has the lowest LUMO energy (−1.50 eV) and thus the least resistance to electron acceptance due to its high electronegativity (4.05 eV), low hardness (2.56 eV), high smoothness (0.39 eV) and the best index of electrophilicity (3.21 eV). Concerning the chemical hardness, a large difference in energy gap between HOMO and LUMO characterises a molecule as being hard and is related to more stable molecules, while a small difference in energy gap corresponds to a molecule as being soft and it is related to more reactive molecules [33] .

It can be seen, in Figure 3 , the energy gap values for the imidazolic alkaloids are very close; however, the molecules of 13N (5.11 eV), PS (5.37 eV), EPS (5.44 eV) and IPS (5.44 eV) have the smallest chemical gaps and hardness in relation to the others molecules and it can be considered the most reactive.

After obtaining the electronic and energetic properties of the imidazolic alkaloids present in the extract of the leaves of Jaborandi, the molecular docking was performed with the SARS-CoV-2 M pro protease, responsible for the severe acute respiratory syndrome in humans. Table 3 shows the results of the binding energies found for the imidazolic alkaloids and the N3 critalographic ligand. As a matter of demonstrating the activity with the active site of the protease M pro , Table 3 shows the ligands EPR (ΔGbind = −7.0 kcal mol −1 ) and EPS (ΔGbind = −6.8 kcal mol −1 ) have the best interaction energies with the peptidase. Therefore, EPR was the alkaloid that came closest to the binding energy carried out by the crystallographic ligand N3 (ΔGbind = −8.3 kcal mol −1 ), accompanied by the isomers EPS, IPS (ΔGbind = −6.8 kcal mol −1 ) and PS (ΔGbind = −6.7 kcal mol −1 ). Meanwhile, the PC has the highest interaction energy value (ΔGbind = −5.0 kcal mol −1 ). Because they are isomeric compounds, the results of the binding energies for these molecules are very close, with similar interactions at the active site of the receptor (Table 4 ).

In the in silico study carried out by Huynh et al. [1] , the interactions of the molecules (chloroquine, bromhexine, favipiravir, dipyridamole, ambroxol, hydroxychloroquine, montelukast, cinaserin, GS-441524, kaempferol, lopinavir, entecavir, umifenovir, quercetin, remdesivir, nelfinavir, curcumin, and N3), were analysed along with the SARS-CoV-2 protease M pro through molecular docking with the AutoDock Vina software. In which, N3 was considered as the control molecule, which showed an important interaction at the binding site with one of the best fitting scores, estimated at −7.1 kcal mol −1 . In this sense, the relevance of the binding energies presented by the imidazolic constituents of the Jaborandi leaf is perceived, mainly the molecules that are isomers to EPR.

It is possible to establish a relationship between a molecular docking study carried out by Barros et al. [49] , which was relevant when investigating the interaction activity of 24 ligands (based on commercialised drugs), and four SARS-CoV-2 receptors, presenting a M pro with a stable ΔGbind when coupled with the drugs Pemirolast (−5.8 kcal mol −1 ), Benserazide (−4.9 kcal mol −1 ) and Luteolin (−6.1 kcal mol −1 ). The results presented are as important as those of the aforementioned study, since the energy of interaction of the imidazolic alkaloids EPR, EPS, IPS and PS of Jaborandi with M pro (Table 3) were very favourable to the receptor-ligand interaction. Table 4 shows the best molecular affinity results, for imidazolic alkaloids and their respective intermolecular interactions (hydrogen, electrostatic, hydrophobic, and van der Waals bonds), with some amino acid residues in the active site of the main peptidase of COVID-19.

By drawing a parallel based on the crystallographic results of the SARS-CoV-2 protease M pro , each protomer is formed by the junction of three domains: Domain I (residues 8-101), Domain II (residues 102-184) and domain III (residues 201-303). The peptidase presents a catalytic dyad composed by the amino acid residues: Cysteine and Histidine; and the substrate binding site represented by a cavity between domains I and II [8] . In brief, it is noticed that the ligands in Table 4 are located in the same region and present similar interactions with the main residues of the active protease M pro site.

The EPS and IPS ligands (Table 4 and Figure 4 ) when interacting with the protease, present hydrogen bonds with the His41 and Gln189 residues, and electrostatic interactions with His41 and His163; the hydrophobic contacts with His41, Cys145, and His163; however, they differ regarding the van der Waals bonds with the amino acids Leu141 and Ser144, respectively. While in Figure 5 , the PS shows a similarity in relation to the hydrogen bonding to Gln189, the electrostatic interactions with His41 and His163, the hydrophobic contacts with His41 and Cys145, and the van der Waals bonds (His164, Gly143, Asn142, Phe140, His172, Glu166, Tyr54, Arg188, and Asp187). EPR makes contacts with other residues, differentiating in hydrogen bonds when interacting with Thr190 and Gln192, having no electrostatic interactions, presenting hydrophobic contacts with Met49, Met165, and Gln189, and van der Waals connections with Ala191, Leu167, Glu166, and Asp187. In Figure 5 (a and b) , all interactions performed after molecular docking are illustrated. It is noticed in the EPR (Figure  5a) , the existence of a hydrogen bond between the oxygen (O3) of the organic ester function of the γ-butyrolactone ring with the Gln192 residue, and a second bond occurs between the hydrogen (H29) of the hydroxyl group of the alcohol organic function with the amino acid Thr190. While in PS (Figure 5b) , a hydrogen bond is created between the oxygen (O1) of the organic ester function of the γ-butyrolactone ring with the Cys145 residue, and a second bond between the hydrogen (H29) of the hydroxyl group of the alcohol organic function with the amino acid Gln189. This fact is observed in the crystallographic structure of the protease M pro, in which, through electronic density, the crystallographic ligand N3 performs a 1.8 Å covalent bond with the Cys145 of the binding site [8] .

It is noticed that the EPR and EPS molecules have the best interaction affinities, based on the binding energies with the enzyme M pro (Table 3 ). This can be justified by the location, the type of connection, and the position of the ligands in the active site. According to Huynh et al. [1] , after studies of docking and molecular dynamics with 19 drug molecules, it was observed that the 'anchor' site at the protease binding site M pro , plays an important role in the stabilisation of ligands, caused by hydrophobic interactions performed. In the crystallographic structure of the M pro peptidase from SARS-CoV-2 complexed to N3, one end of the ligand occupies the 'anchor' site of the active site, which was revealed as an important region for the interaction with the receptor [1] .

Agents capable of inhibiting the SARS-CoV-2 protease M pro are essential to blocking viral replication [5] . In this binding site, there are the amino acid residues His41, Met49, Gly143, Cys145, His163, His164, Glu166, Pro168, and Gln189, as shown in the recent study on α-ketoamide inhibitors [5] . In the present study, it was possible to identify the presence of the same amino acid residues by performing different interactions with the ligands EPR, EPS, IPS, and PS, thus being considered as strong candidates for the enzymatic inhibition of M pro in COVID-19.

Molecular dynamics calculations were performed based on the most stable conformations of the ligands determined by molecular docking. The M pro -EPR and M pro -EPS complexes were selected for the 50 ns simulations, using the amino acid residue Met165 as a reference located at the active site of the protease. Thus, it can be seen, in Figure  6 , that the EPR ligand has a good interaction with the site, very close to the selected residue, in the intervals of 6-20 ns, 24-27 ns, 29-35, 42 ns and 44-50 ns, and remaining out of that location in two short time intervals 21-23, 28 ns, 36-41 and 43 ns. Throughout this simulation, a constant change in the conformation of the M pro -EPR complex is noticeable. All things considered, the results reproduce the system well, presenting a good connection affinity in most of the simulation.

It is observed, in Figure 7 , the behaviour of the binding affinity in the M pro -EPS complex during the simulation time, in which the ligand in the period from 6 to 50 ns shows a good interaction, remaining in the binding site. The occurrence of a high thermodynamic stability is shown, due to the energetic balance between the ligand and the active site of the protease, then providing the stabilisation of the complex attributable to the interactions carried out with the residues of that region.

The M pro -EPS complex, in Figure 7 , presents different conformations as the simulation time increases, with that, the ligand adjusts geometrically to the active site of the protease as it undergoes structural changes. Figure 8 shows the interaction times in nanoseconds of the M pro -EPS and M pro -EPR complexes during the simulation.

Based on the analysed data, the EPR and EPS molecules when interacting with the SARS-CoV-2 protease M pro have a dynamic and reversible character, caused by the flexibility and structural reorganisation capacity of the complexes during the simulation. Then, when the ligands interact with the peptidase M pro , they assume a new structural conformation. However, when this equilibrium condition is affected, it causes a destabilisation of the system and subsequent expulsion of the ligand, as occurred with the EPR. When establishing a new balance, the energy of the system returns to its favourable condition of interaction with the ligand; allowing to return it to the region of the active site of the protease [76] . Table 5 shows the values of the energies of electrostatic interactions (Coulomb), by van der Waals (Lennard-Jones) and the sum of these (Coul + LJ), calculated for the M pro complexes of SARS-CoV-2 with the ligands EPS and EPR, during the time of analysis of the simulation by molecular dynamics.

These results in Table 5 show that the EPS ligand features a better Coul energetic interaction (−50.4488 ± 4.9 kJ mol −1 ) in relation to EPR (−15.5741 ± 4.4 kJ mol −1 ) when interacting with the enzyme M pro , and this fact is caused by the electrostatic character of the hydrogen bonds formed in the complex. The LJ interactions are present in the two complexes, in which the EPS ligand has the lowest energy (−143.973 ± 2.9 kJ mol −1 ) compared to the EPR (−122.483 ± 6.0 kJ mol −1 ). Therefore, we can observe the same behaviour with the approximate sum of the binding energies (Coul + LJ), as the EPS (−194.422 ± 5.7 kJ mol-1) presents greater thermodynamic stability in relation to EPR (−138.057 ± 7.4 kJ mol −1 ), when the M pro protease of SARS-CoV-2 is complexed. This behaviour can be explained by the longer interaction time and better conformational balance with the binding site. Although EPS is thermodynamically more stable, EPR showed considerable interactions with peptidase M pro . However, the results for the energies in the interaction process between the ligand and the active site of the receptor are not the only determining factor for binding affinity.

The ADMET results of the physical-chemical and pharmacokinetic properties of the EPR, EPS, IPS, and PS molecules are shown in Tables 6-8. It is important to highlight that the EPR, EPS, IPS, and PS isomers have similar or even equal parameters as seen in this section.

The ideal molar solubility in water (logS) should be greater than −4 and less than 6, based on a qualitative solubility estimate provided by the logS scale [68] . Whereas the distribution constant (logD) is a descriptor of lipophilicity, adjusted to a pH using a tampon, its ideal value tends to be greater than 3 [68, 77] . In effect, it can be seen in Table 6 that the compounds have low lipophilicity and reasonable hydrophilicity. According to Storpirtis et al. [78] , a soluble molecule facilitates many activities for the development of drugs, especially in relation to handling and formulation. On the other hand, when it is insoluble it can hinder the dissolution process in biological fluids and the development of some pharmaceutical forms.

TPSA (Topological Polar Surface Area) is a descriptor commonly related to hydrogen bonds, which can cause changes in permeability and oral bioavailability. It can be used in combination with counting rotatable bonds and assessing molecular flexibility, which can also alter the oral bioavailability of many molecules [79] . The polarity or polar topological surface area must be between 20 and 130 Å2; thereby, the values in Table 6 are in accordance with the ideal value [80] .

The Hydrogen Bonding Donor (HBD) and Hydrogen Bonding Acceptor (HBA) are two important factors. Studies show that compounds with a presence of HBA, greater than HBD, are favourable for ADMET conditions, considering that a large number of HBD may result in low permeability, absorption, and bioavailability [81] . Thus, it is observed that the analysed isomers obey the rule (HBD ≤5 and HBA ≤10).

Many rules have been developed to guide the selection of drug candidate compounds [82] . The similarity of drugs was established from structural or physical-chemical inspections of the compounds, the advanced and sufficient development to be considered candidates for oral drugs, according to the parameters established in the filters.

Lipinski's rule, also known as the rule of five, outlined the relationship between pharmacokinetics and physical-chemical The parameter Fsp3 defines the molecular complexity, the number of sp3 carbons hybridised in the total count [68] . All compounds have a higher value than recommended, which can imply high aqueous solubility.

The compounds also have a bioavailability score in which the closer to 1, the higher the quality of the molecule as a potential drug candidate [83] . For this reason, the data are favourable to the analysed compounds. Table 7 reproduces the parameters ADME and Toxicity of imidazolic alkaloids.

In Table 7 , the Rule of Five and the WDI like rule (PW≤550; HBA≤ 9 and HBD≤5; logP≤5) indicate these molecules as suitable. The last is based on compounds with molecular properties with a limit equal to or greater than 90% of those found in the World Drug Index (WDI), then, the molecules are in accordance with the rule, as they have a 90% cut [84] .

The BBB indicates whether the compound is able to overcome the barrier between blood and cerebrospinal fluid [85] . Then, the present molecules have a very low penetration capacity. Caco-2 cells exhibit many of the morphological and functional properties of human intestinal cells, whereas MDCK cells can be used as a good tool for rapid permeability screening [86] . In ADME terms, the results show Caco-2 with average values and MDCK with small values.

The HIA data are calculated from bioavailability and absorption, being a crucial factor in predicting the viability of medication absorption through the small intestine [87] . Low-absorbed compounds have HIA of 0-20%, moderately absorbed 20-70%, and highly absorbed 70-100%. The values, in Table 7 , refer to compounds that have a high human intestinal absorption.

While the prediction of toxicity is an indispensable factor in the screening of new drugs, the results are favourable as the tests ran on rats and mice were negative. The risk of heart disease (hER risk) is an important factor; and the molecules here analysed present a medium risk. The assessment of in silico toxicity in drug planning is an important factor, as it assists in determining the toxic dose in studies with animal models and decreases the number of guinea pigs infected by the disease [88] .

In Table 8 , PASS Online predicts both the probability of being active (Pa) and the probability of being inactive (Pi). Pa estimates the chance of the studied compound to belong to the subclass of active compounds and Pi estimates the chance of the compound to belong to the subclass of inactive compounds [69] .

The compounds have potential activity against Picornavirus, which the best known of them is Rhinovirus, the one that causes the common cold in humans and is associated with severe respiratory infections [89] . They also have antioxidant effects, playing a crucial role in protecting healthy body cells against the oxidising action of free radicals [90] . And they are considered anticancer agents, through the negative tests ran on rats and mice.

Therefore, the ADMET in silico results for the EPR, EPS, IPS and PS molecules are promising for the continuity with the in vitro and in vivo tests with the SARS-CoV-2 strains. This is because they present activity against serious respiratory infections caused by Rhinovirus and absence of toxicity in the organisms of rats and mice.

Based on the results presented in this study, the imidazole alkaloids epiisopiloturine, epiisopilosine, isopilosine and pilosine have shown indications as possible inhibitors of the main SARS-CoV-2 protease. This is a fact observed through the high thermodynamic stability, chemical reactivity based on quantum descriptors (especially epiisopilosine, isopilosine and pilosine), the strong interactions with the main amino acid residues in the active site of the enzyme M pro , the constant interaction activity during the 50 ns of simulation, the favourable pharmacokinetic profile and the absence of toxicity for rats and mice.

However, only the epiisopilosine and epiisopiloturine molecules were investigated in a dynamic manner, what is promising regarding the energy stability and the affinity of interaction with the target. In this way, the in silico data is relevant and allows the continuation of the steps in vitro and in vivo in the search for a treatment or cure for COVID-19. 

In silico exploration of the molecular mechanism of clinically oriented drugs for possibly inhibiting SARS-CoV-2's main protease

A novel coronavirus from patients with pneumonia in China

Global coalition to accelerate COVID-19 clinical research in resource-limited settings

Fact sheet -COVID-19 (disease caused by the new coronavirus)

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved a-ketoamide inhibitors

Structure of coronavirus main proteinase reveals combination of a chymotrypsin fold with an extra α-helical domain

The crystal structures of severe acute respiratory syndrome vírus main protease and its complex with an inhibitor

Structure of M pro from SARS-CoV-2 and discovery of its inhibitors

Design of wide-spectrum inhibitors targeting coronavirus main proteases

Structure of main protease from human coronavirus NL63: insights for amplo espectro anti-coronavirus drug design

Moroccan medicinal plants as inhibitors against SARS-CoV-2 main protease: computational investigations

The use of medicinal and phytotherapy plants: An integrational review on the nurses performance

Seasonal change in main alkaloids of jaborandi (Pilocarpus microphyllus Stapf ex Wardleworth), an economically important species from the Brazilian flora

Handbook of medicinal herbs

Antidiabetic activity of pilocarpus microphyllus extract on streptozotocin-induced diabetic mice

Nanoparticle cross-linked collagen shields for sustained delivery of pilocarpine hydrochloride

Treatment of xerostomia and hyposalivation in the elderly: a systematic review. Oral Medicine Oral Pathology and Oral Surgery

Nanopharmaceutical approach of epiisopiloturine alkaloid carried in liposome system: preparation and in vitro schistosomicidal activity

Anthelmintic activity in vivo of epiisopiloturine against juvenile and adult worms of Schistosoma mansoni

Computational quantum chemistry, molecular docking, and ADMET predictions of imidazole alkaloids of Pilocarpus microphyllus with schistosomicidal properties

Anti-inflammatory and antinociceptive activity of epiisopiloturine, an imidazole alkaloid isolated from Pilocarpus microphyllus

Epiisopiloturine hydrochloride, an imidazole alkaloid isolated from Pilocarpus microphyllus leaves, protects against naproxen-induced gastrointestinal damage in rats

Anthelmintic, antibacterial and cytotoxicity activity of imidazole alkaloids from Pilocarpus microphyllus leaves

Estudo em complexos fármaco-receptor utilizando bioquímica quântica [Quantum biochemistry study in drug-receptor complexes

An investigation into the identification of potential inhibitors of SARS-CoV-2 main protease using molecular docking study

Recent Trends in drug likeness prediction: a comprehensive review of In silico methods

The protein data Bank

Pilocarpus spp.: a survey of its chemical constituents and biological activities

Characterization of the variation in the imidazole alkaloid profile of Pilocarpus microphyllus in different seasons and parts of the plant by electrospray ionization mass spectrometry fingerprinting and identification of novel alkaloids by tandem mass spectrometry

Industrial scale Isolation, structural and spectroscopic characterization of epiisopiloturine from pilocarpus microphyllus Stapf leaves: a promising alkaloid against schistosomiasis

GaussView, version 6

Structural, vibrational, UV-vis, quantum-chemical properties, molecular docking and anti-cancer activity study of annomontine and N-hydroxyannomontine β-carboline alkaloids: a combined experimental and DFT approach

Quantum chemical properties investigation and molecular docking analysis with DNA topoisomerase II of β-carboline indole alkaloids from Simaba guianensis: a combined experimental and theoretical DFT study

Development of the colle-salvetti correlation-energy formula into a functional of the electron-density

Density-functional thermochemistry. III. The role of exact exchange

Contracted Gaussian-basis sets for molecular calculations. 1. second row atoms, Z = 11-18

Self-consistent molecular orbital methods. XX. Basis set for correlated wave-functions

Proof that in density-functional Theory

Density-functional theory for fractional particle number: derivative discontinuities of the energy

Principle of maximum hardness

UCSF Chimera-a visualization system for exploratory research and analysis

Automated docking of flexible ligands: applications of AutoDock

Using AutoDock4 and AutoDock Vina with AutoDockTools: A tutorial. The Scripps Research Institute

Autodock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading

AutoDockFR: advances in protein-ligand docking with explicitly specified binding site flexibility

Automated docking using a Lamarckian genetic algorithm and an empirical binding free energy function

Interaction of drugs candidates with various SARS-CoV-2 receptors: an in silico study to combat COVID-19

Discovery studio visualizer

Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis

Atomic charges derived from electrostatic potentials: a detailed study

Numerical study on the partitioning of the molecular polarizability into fluctuating charge and induced atomic dipole contributions

Hþþ: a server for estimating pKas and adding missing hydrogens to macromolecules

The GROMACS development team

A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6

Structurefunction studies of BPP-BrachyNH 2 and synthetic analogues thereof with Angiotensin I-converting enzyme

A systematic study of water models for molecular simulation: derivation of water models optimized for use with a reaction field

Interaction of wild type, G68R and L125M isoforms of the arylamineN-acetyltransferase from mycobacterium tuberculosis with isoniazid: a computational study on a new possible mechanism of resistance

Constant pressure molecular dynamics for molecular systems

Polymorphic transitions in single crystals: a new molecular dynamics method

LINCS: a linear constraint solver for molecular simulations

Particle mesh Ewald: an N, log(N) method for Ewald sums in large systems

From proteins to perturbed hamiltonians: a suite of tutorials for the GROMACS-2018 molecular simulation package

Bioinformática da biologia à flexibilidade molecular [Bioinformatics from biology to molecular flexibility

FAF-Drugs4: free ADME-tox filtering computations for chemical biology and early stages drug discovery

The PreADME: pc-based program for batch prediction of adme properties

SwissADME: a free web tool to evaluate pharmacokinetics, drug-likeness and medicinal chemistry friendliness of small molecules

Prediction of the biological activity spectra of organic compounds using the PASS online web resource

Learning medicinal chemistry absorption, distribution, metabolism, excretion, and toxicity (ADMET) rules from cross-company matched molecular Pairs analysis (MMPA)

Termodinâmica: Uma Proposta de Ensino a partir da Química Computacional

An AM1 study on the electron-donating and electron-accepting character of biomolecules

On the limits of highestoccupied molecular orbital driven reactions: the frontier effectivefor-reaction molecular orbital concept

Advanced QSAR studies on PPARδ ligands related to metabolic diseases

Two rules on the protein-ligand interaction

Constituents of buriti oil (Mauritia flexuosa L.) like inhibitors of the SARS-coronavirus main peptidase: an investigation by docking and molecular dynamics

The influence of lipophilicity in drug discovery and design

Farmacocinética básica e aplicada

Exploring the physicochemical profile and the binding patterns of selected novel anticancer himalayan plant derived active compounds with macromolecular targets

Fast calculation of molecular polar surface area as a sum of fragment-based contributions and its application to the prediction of drug Transport properties

Computational analysis of calculated physicochemical and ADMET properties of proteinprotein interaction inhibitors

Quantifying the chemical beauty of drugs

Combinatorial library design for diversity, cost efficiency, and drug-like character

Development of a new system for prediction of drug absorption that takes into account drug dissolution and pH change in the gastro-intestinal tract

Madin-Darby Canine Kidney) cells: a tool for membrane permeability screening

In vitro models for the determination of drug absorption and a prediction of dissolution/absorption relationships

In silico evaluation of multispecies toxicity of natural compounds

Rinovírus humano em infecções respiratórias agudas em crianças menores de cinco anos de idade: fatores envolvidos no agravamento da doença [Human rhinovirus in acute respiratory infections in children under five years of age: factors

Application of several chemical methods to determine antioxidant activity in fruit pulp

 We are grateful for the institutional support of CENAPAD-UFC, IFPI,  IFPA, IFMA, UFPI and 

No potential conflict of interest was reported by the author(s).