key: cord-0016686-i50lps17
authors: Afjalus Siraj, Md; Rahman, Md. Sajjadur; Tan, Ghee T.; Seidel, Veronique
title: Molecular Docking and Molecular Dynamics Simulation Studies of Triterpenes from Vernonia patula with the Cannabinoid Type 1 Receptor
date: 2021-03-30
journal: Int J Mol Sci
DOI: 10.3390/ijms22073595
sha: 4465bf32edf389d094722c60ebe8ceea40f31f13
doc_id: 16686
cord_uid: i50lps17

A molecular docking approach was employed to evaluate the binding affinity of six triterpenes, namely epifriedelanol, friedelin, α-amyrin, α-amyrin acetate, β-amyrin acetate, and bauerenyl acetate, towards the cannabinoid type 1 receptor (CB1). Molecular docking studies showed that friedelin, α-amyrin, and epifriedelanol had the strongest binding affinity towards CB1. Molecular dynamics simulation studies revealed that friedelin and α-amyrin engaged in stable non-bonding interactions by binding to a pocket close to the active site on the surface of the CB1 target protein. The studied triterpenes showed a good capacity to penetrate the blood–brain barrier. These results help to provide some evidence to justify, at least in part, the previously reported antinociceptive and sedative properties of Vernonia patula.

The cannabinoid receptors (CB) belong to the superfamily of G protein-coupled receptors and are divided into two major types: CBR type-1 (CB1) and CBR type-2 (CB2). The CB1 receptors are commonly found in the central nervous system (CNS) and mostly control pain, movement, and neurotic parameters [1] [2] [3] . CB1 can be also found in peripheral tissues including retina [4] , colon [5] , testis [6] , sperm cells [7] , and adipocytes [8] . In contrast, CB2 receptors are mostly found in peripheral tissues [9] . Cannabinoids have already demonstrated great potential for the treatment of pain, inflammation, and neurodegenerative disorders [10] [11] [12] . These include the phytocannabinoids from Cannabis sativa and other plant-derived cannabinoid-like molecules called cannabimimetics that can interact with the endogenous cannabinoid system [13] .

Vernonia patula (Dryand.) Merr. (Asteraceae) (VP) is an annual herb widely distributed throughout Southeast Asia. It is used medicinally for malaria, colds, fevers, respiratory ailments, and convulsions [14, 15] . The leaves are used for their analgesic properties [16] . The aerial parts of VP have displayed antioxidant and anti-inflammatory activity [17] . Several Vernonia species, including VP, have also demonstrated antinociceptive and sedative properties [18] [19] [20] . The aerial parts of VP are known to contain flavonoids, simple phenolics, and terpenoids [21] [22] [23] . Previous reports have indicated that triterpenes could act as CB1 receptor agonists [24] . In the present work, we conducted molecular docking and molecular dynamics simulation studies on six triterpenes previously reported in the aerial parts of VP against the CB1 receptor with a view to (i) validating the medicinal properties of this plant and (ii) identifying new plant-derived cannabimimetic drug templates.

The results of the HPLC-DAD-MS analysis (Figures S1-S3) indicated the presence of epifriedelanol (1), friedelin (2) , α-amyrin (3) and α-amyrin acetate (4), β-amyrin acetate (5) , and bauerenyl acetate (6) in VP.

All triterpenes had a molecular weight of less than 500 g/mol. Their Moriguchi LogP (MLogP) values were higher than those of the tetrahydrocannabinol (THC) control, while their polar surface areas were lower than THC (Table 1) . All compounds showed high lipophilicity and insolubility in the Bioavailability Radar plots ( Figure S4 ). Friedelin and α-amyrin demonstrated a similar human intestine absorption prediction score to that of THC (0.99). All triterpenes showed a good capacity to penetrate the blood-brain barrier. Friedelin showed the highest blood-brain barrier penetration prediction score (0.99), comparable to THC (1.00). Table 1 . Drug-likeness parameters prediction for tetrahydrocannabinol (THC), epifriedelanol (1), friedelin (2), α-amyrin (3) and α-amyrin acetate (4), β-amyrin acetate (5), and bauerenyl acetate (6) . 

Alpha-amyrin, friedelin, and epifriedelanol displayed the best binding affinities for CB1 with scores of −8.3, −8.1, and −8.1 kcal/mol, respectively. The THC control had a predictive binding affinity of −9.4 kcal/mol. The amino acid residues of CB1 involved in the binding with THC, epifriedelanol (1), friedelin (2) , and α-amyrin (3) with their bonding distances are depicted in Table 2 . The predictive binding affinity of the remaining triterpenes towards CB1 are shown in Table S1 . The intermolecular binding interactions of epifriedelanol-CB1, friedelin-CB1, and α-amyrin-CB1 are depicted in Figure 1 . THC showed 16 interactions, including one conventional H-bond (THC-OH·Met384-S), one Pi-Sigma, one amide-Pi stacking, and six alkyl and seven Pi-Alkyl bonds. It formed nonbonding interactions with both Phe102 and Phe379 residues within the active site of the CB1 protein. Epifriedelanol formed one conventional H-bond (Arg182-NH 2 ·Epifriedelanol-O), and 13 hydrophobic (9 alkyl and 4 Pi-alkyl) non-bonding interactions. Friedelin formed 13 hydrophobic (9 alkyl and 4 Pi-alkyl) non-bonding interactions. Both epifriedelanol and friedelin were found to bind to a pocket containing residues Ile175, Tyr172, and Phe180 that were close to the active binding site of the THC control. Alpha-amyrin formed 12 hydrophobic (6 alkyl and 6 Pi-alkyl) non-bonding interactions with residues Ala120, Phe177, His178, and Val179, also close to the active binding site. 

Molecular dynamics simulations were conducted to further analyze the stability and binding affinities of the triterpene-CB1 complexes compared to the THC-CB1 complex. As both friedelin and epifriedelanol are structurally similar, only friedelin was selected for molecular dynamics simulations. Alpha-amyrin was also considered as it showed the best docking score among all triterpenes. Alpha-amyrin acetate was selected as a representative for acetate-containing compounds. The total potential energies of CB1, CB1-THC, CB1friedelin, CB1-α-amyrin, and CB1-α-amyrin acetate were calculated for a period of 100 ns. Values obtained for the energy-minimized initial conformations at 0 ps were −1,024,595.105, −1,090,152.423, −1,045,113.241, −1,056,698.855, and −992,631.612 kJ/mol, respectively. The values obtained for each studied system after a few picoseconds were −790,249.746, −842,086.509, −804,131.370, −815,489.160, and −761,594.499 kJ/mol, respectively. The total potential energies of the studied systems remained stable throughout the 100 ns simulation period (Figure 2 ). The CB1-THC complex showed the lowest potential energy, successfully binding at the active site of CB1. The triterpenes were found to bind strongly to a pocket at the surface of the target protein, close to the active binding site of THC. The CB1-friedelin and CB1-α-amyrin complexes showed lower potential energy values than CB1 alone. The CB1-α-amyrin acetate complex showed higher energy than CB1 alone. All complexes were stable and in compliance with the docking scores obtained. 

The PCA scores plot reveals different clusters formation of the protein and the proteinligand complexes based on their structure and energy profiles (Figure 3a ). The cluster for the α-amyrin-CB1 complex overlapped the one for the friedelin-CB1 complex. On the PC2 axis, the friedelin-CB1 cluster showed a similar pattern to the CB1 cluster, but more compactness. The α-amyrin acetate-CB1 complex formed a distinct cluster showing negative correlation on the PC1 axis differing from the other triterpene-CB1 and THC-CB1 complexes. The cluster of the THC-CB1 complex was further away from the CB1 cluster compared to the other complexes. The loading plot generated from the molecular dynamics (MD) energy profiles, and structural data showed the bond, angle, and Van der Waals energies and RMSD-Cα values as the major contributing factors for the stabilizing of the protein and the complexes (Figure 3b) . A scree plot showing the eigenvalues further validated the analyses (Figure 3c ). 

The atomic root mean square deviations (RMSDs) of Cα atoms of CB1 and of the selected CB1-triterpene complexes were analyzed to compare their structural stability ( Figure 4 ). All complexes reached equilibrium after 50 ns. The α-amyrin acetate-CB1 complex exhibited the lowest RMSD value among all complexes. A decrease in RMSD was observed for this complex from 80 to 100 ns, suggesting that this complex may not be stable overall. The RMSD value of the α-amyrin-CB1 complex showed a more stable pattern than the corresponding acetate complex. The RMSD values obtained for the α-amyrin-CB1 and the α-amyrin acetate-CB1 complexes were lower than those of CB1 alone, which suggest that α-amyrin and α-amyrin acetate contributed to lower the energy of the protein and made the protein more stable. The RMSDs of the CB1-THC and CB1-friedelin complexes fluctuated greatly from 10 to 60 ns and then stabilized. The value of the CB1-THC complex reached 2.0-2.2 Å after 60 ns. The RMSD of the CB1-friedelin complex fluctuated mostly between 10 and 20 ns, then showed a steady increase over the simulation period. These relatively higher deviations compared to the RMSD of the CB1 structure suggested that THC and friedelin may change the protein conformation at the binding site.

The binding of THC to CB1 was found to primarily induce local flexibility of the active site residues and revealed that the protein became more flexible in all regions. In contrast, the α-amyrin-CB1 and the friedelin-CB1 complexes showed the lowest RMSF, indicating that the binding of α-amyrin and friedelin to CB1 made the protein less flexible. The α-amyrin acetate-CB1 complex showed a similar trend to the THC-CB1 complex in terms of residue motility analysis ( Figure 5 ). 

Additional calculations of the binding free energies of the CB1-ligand complexes investigated in the molecular dynamics simulations using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method are presented in Figure 6 . Alpha-amyrin acetate showed the highest binding free energy value (−36.6 ± 5.09 Kcal/mol) among the selected triterpenes. 

It has previously been demonstrated that the aerial parts of VP are antinociceptive and can induce sedation by suppressing locomotor activity and exploratory behavior in mice. This has been linked with the presence of phenolic compounds predicted to interact with CB1 [20, 23, 25] . However, it has also been reported that triterpenes, including α-amyrin and β amyrin, could induce in vivo antinociceptive and anti-inflammatory effects via activation of the cannabinoid receptors [24, 26] . Other terpenoids with effects on the CNS include α-pinene, which suppresses locomotor activity, increases sleeping time, and produces an anxiolytic effect in vivo [27] , and phytol and terpinolene with sedative and locomotor suppressive effects [28, 29] . Similar sleep-inducing and locomotor relaxant effects have been reported for citral, limonene, and myrcene [30] . A decrease in locomotor activity can be correlated to a potential CNS depression [31] . Myrcene and α-pinene have been shown to increase the GABA A receptor activity in vitro and potentiate the release of inhibitory neurotransmitters [32] . Moreover, it was also reported that GABA A -stimulating drugs such as flurazepam can synergistically potentiate the cataleptic effects of THC in mice [33] .

The purpose of our computational studies was to evaluate the predictive binding affinity of VP triterpenes towards the CB1 receptor and carry out molecular dynamics simulations to describe the nature of the interactions of these triterpenes with CB1. Friedelin, α-amyrin, and epifriedelanol showed a strong binding affinity for the CB1 receptor in our molecular docking study. Molecular dynamics simulations, through stability and residue mobility analyses, was used to understand the structural variations and conformational flexibility of the CB1 protein and CB1 complexed with the selected triterpenes. All the studied triterpenes showed stable binding throughout the molecular dynamics simulations. The most stable interactions with CB1 were predicted for α-amyrin and friedelin. The latter had a blood-brain barrier penetration prediction score comparable to THC. Our molecular docking and molecular dynamics simulation results suggest that by interacting with CB1 receptors, VP triterpenes could contribute to the significant antinociceptive and CNS-depressant activity we have previously demonstrated for this plant [20] .

The aerial parts of VP were collected from Chittagong (Bangladesh) in 2015. The plant material was identified by M.A. Ali at the Bangladesh National Herbarium in Dhaka where a voucher specimen (DACB: 35107) is kept for future reference. The dried powdered aerial parts of VP (100 g) were macerated in ethyl acetate (EtOAc; 1.5 L) for 5 days at 25 ± 2 • C, with occasional shaking. The resulting filtrate was concentrated under reduced pressure to afford the final extract (6.2 g, 6.2% yield).

The extract (10 mg) was subjected to solid phase extraction (SPE) using a Hypersep C18 cartridge (Thermo Scientific, Waltham, MA, USA) eluting sequentially with 70%, 80%, 90%, and 100% methanol (MeOH). Fraction 1 (1.5 mg), fraction 2 (1.8 mg), fraction 3 (2.1 mg), and fraction 4 (4.2 mg) were collected, and stock solutions (0.1 mg/mL) were prepared for further analysis. The sample solution and blank control were all stored at 4 • C and filtered through a 0.22 µm PTFE syringe filter (Thermo Scientific, Waltham, MA, USA) prior to analysis by HPLC-DAD-MS. Analysis was carried out using an Agilent HPLC 1260 binary pump and a Sunfire C18 column (2.1 × 150 mm, 3.5 µm) (Waters, Milford, MA, USA). The diode array detector (DAD) was set at λ = 210, 254, 269, and 310 nm. The known, or suspected, major adducts for all samples were registered in the positive electrospray ionization (ESI) mode either as [M+H] + or [M+Na] + . The mobile phase used was a gradient of 1% formic acid in water and acetonitrile (solvents A and B, respectively) at a flow rate of 0.2 mL/min. Elution was performed as follows; Fraction 1 (0 to 30 min-60 to 80% solvent B), fraction 2 (0 to 30 min-70 to 90% solvent B), fraction 3 (0 to 30 min-80 to 100% solvent B), and fraction 4 (0 to 30 min-100% solvent B).

Absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties of the selected triterpenes were determined using the online ADMET structure-activity relationship database (admetSAR). The specific pharmacokinetic profile for each compound was obtained using compound specific SMILES strings [34] . MLogP was used as an alternative to the regular LogP model [35] . In addition, drug-likeness and molecular properties of the triterpenes were calculated using SwissADME, and Bioavailability Radar plots were generated taking account of lipophilicity, size, polarity, solubility, flexibility, and saturation [36] .

As triterpenes had been reported as CB1 receptor agonists [23] , we decided to consider six triterpenes, namely epifriedelanol (1), friedelin (2), α-amyrin (3), α-amyrin acetate (4), β-amyrin acetate (5), and bauerenyl acetate (6) ( Figure S5 ), previously reported in the aerial parts of VP [4, 5] as the putative ligands for our in silico analyses. Optimizations for the triterpenes and vibrational frequency calculations were determined at gaseous phase using the Gaussian 09 software version A.02 [37] with a semi-empirical PM6 method [38] . The three dimensional structures of each compound were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/, accessed on 26 March 2021), geometry-optimized, and saved in PDB format using GaussView v.5.0 (https://gaussian.com/gaussview6/, accessed on 26 March 2021).

The crystal structure of the human CB1 cannabinoid receptor and its active site [39, 40] was retrieved from the RSCB Protein Data Bank (PDB ID: 5U09). The heteroatoms and water molecules were removed from the crystal structure using PyMOL Molecular Graphics System v. 1.3 (https://pymol.org, accessed on 26 March 2021) [41] . The structure of the protein was further optimized using the Swiss-PDB viewer software (v.4.1.0) based on the energy minima. The protein and ligand structures were saved in the PDBQT format.

The active binding pocket of CB1 was predicted by CASTp (v. 3.0) [42] . The protein showed the highest pocket area and volume at 652.65 Å 2 and 331.44 Å 3 , respectively. The grid box was generated so as to include the CB1 binding site residues, with a center set at 12.5230, 7.2495, and 17.7743 Å and a size set at 81.5, 81.5, and 81.5 Å in x, y, and z directions, respectively. Autodock Vina (v.1.1.2) was used to perform the molecular docking study [43] . The docked pose of the lowest binding free energy conformer (highest probable binding affinity) with CB1 was analyzed using PyMOL Molecular Graphics System (v. 

The molecular dynamics (MD) simulations were performed on the CB1, CB1-THC, and three selected triterpene-CB1 complexes, namely CB1-friedelin, CB1-α-amyrin, and CB1-α-amyrin acetate. The first two complexes were chosen based on the highest docking scores obtained in the molecular docking study. The third complex was selected to observe the influence of an acetate group on the protein-ligand interactions. MD simulations were conducted using YASARA Dynamics v. 20.8.1 [44] . The AMBER14 force field was employed to simulate the macromolecular system [45] . Each protein was subject to hydrogen bond optimization prior to simulation [44] , and the transferable intermolecular potential 3-point (TIP3P) water model was used by incorporating Cl− and/or Na+ ions. Periodic boundary conditions were incorporated to perform the simulations, where the cell size was 10 Å larger than the protein size in all cases. The initial energy minimization for each system was performed using the steepest gradient approach (5000 cycles), MD simulations were carried out using the particle-mesh Ewald (PME) method to designate long-range electrostatic interactions at a cut off distance of 8 Å and defining physiological conditions at 298 • K, pH 7.4, 0.9% NaCl [46] . The simulation temperature was controlled using a Berendsen thermostat with the pressure kept constant. A multiple time step algorithm was employed with a time step of 2.00 fs [47] . Finally, MD simulations were performed for 100 ns at constant pressure and Berendsen thermostat, and snapshots were saved every 100 ps. Further analysis was conducted using the default YASARA MACRO script [48] .

Principal component analysis (PCA) was used to analyze any subtle variability among the structural and energy profile data obtained from MD simulations for the selected triterpene-CB1 complexes and CB1 alone. Bond energies, bond angle energies, dihedral angle energies, planarity energies, Coulomb energies, Van der Waals energies, and RMSD-Cα values were included as the variables [49] [50] [51] . Multivariate responses were arranged in an X matrix according to the following equation:

where T k describes how the samples relate to each other, P k demonstrates how the variables relate to each other, k is the number of factors included in the model, and E is the matrix of residuals. PCA was conducted using the OriginPro 2021 (Principal Component Analysis app v.1.50) software package.

The root mean square deviations (RMSDs) of Cα atoms on the backbone of the CB1 protein and of the selected CB1-triterpene complexes were calculated during the MD simulation period using YASARA macro file followed by OriginPro 2021. Root mean square fluctuation (RMSF) analysis using YASARA macro script and OriginPro 2021 was used to observe the regions that fluctuated during the MD simulation period.

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method [52] was used to calculate the binding free energies of the CB1-ligand complexes investigated in the molecular dynamics simulations. Default macro scripts of YASARA dynamics were employed for the calculations. Selected snapshots from the last 50 ns MD simulation were used for all CB1-ligand complexes. Protein ligand binding free energy values were calculated using the following equation:

∆G binding = ∆G complex − [∆G ligand + ∆G protein ], and ∆G binding = ∆G MM + ∆G PB +∆G SA-T∆S = (∆G elec + ∆G VdW ) + ∆G PB + ∆G SA-T∆S (2) where ∆G complex = total free energy of the protein-ligand complex in solvent, ∆G ligand = total energy of the ligand in solvent, and ∆G protein = total energy of the protein in solvent. ∆G MM = molecular mechanics interaction energy, where the ∆G elec and ∆G VdW are the electrostatic and Van der Waals interactions, respectively. ∆G PB and ∆G SA represent polar solvation and nonpolar solvation energy, respectively. T∆S (temperature = T and entropy = S) is the contribution of entropy to the free energy.

Dopamine activation of endogenous cannabinoid signaling in dorsal striatum

Cannabinoids, hippocampal function and memory

The endogenous cannabinoid system and its role in nociceptive behavior

Localization of cannabinoid CB1 receptors in the human anterior eye and retina

Differential expression of cannabinoid receptors in the human colon: Cannabinoids promote epithelial wound healing

Molecular cloning of a human cannabinoid receptor which is also expressed in testis

Cannabinoid receptors in sperm

Presence of the cannabinoid receptors, CB1 and CB2, in human omental and subcutaneous adipocytes

Species differences in cannabinoid receptor 2 (CNR2 gene): Identification of novel human and rodent CB2 isoforms, differential tissue expression and regulation by cannabinoid receptor ligands

Targeting cannabinoid agonists for inflammatory and neuropathic pain

The neurobiology of addiction: A neuroadaptational view relevant for diagnosis

Regulation of µ-opioid receptors: Desensitization, phosphorylation, internalization, and tolerance

Cannabimimetic plants: Are they new cannabinoidergic modulators? Planta

Indian Medicinal Plants: An Illustrated Dictionary

A comparative analysis of medicinal plants used by folk medicinal healers in three districts of Bangladesh and inquiry as to mode of selection of medicinal plants

Medicinal plants used by a Tonchongya tribal community at Taknatala village in Rangamati District

Anti-inflammatory and antioxidant activities of ethanolic extract of aerial parts of Vernonia patula (Dryand.) Merr. Asian Pac

Analgesic, antipyretic, anti-inflammatory effects of methanol, chloroform and ether extracts of Vernonia cinerea less leaf

Antinociceptive and anti-inflammatory effects of ethanol extract from Vernonia polyanthes leaves in rodents

Antinociceptive and sedative activity of Vernonia patula and predictive interactions of its phenolic compounds with the cannabinoid type 1 receptor

Chemical and pharmacological studies on Vernonia patula (Dry.) Merr. and Vernonia cinerea (Linn.) Less

Studies on the constituents from the herb of Vernonia patula

Chemical constituents of Vernonia patula

Activation of cannabinoid receptors by the pentacyclic triterpene α, β-amyrin inhibits inflammatory and neuropathic persistent pain in mice

Cannabinoid system involves in the analgesic effect of protocatechuic acid

Euphol, a tetracyclic triterpene produces antinociceptive effects in inflammatory and neuropathic pain: The involvement of cannabinoid system

Intracerebral distribution of α-pinene and the anxiolytic-like effect in mice following inhaled administration of essential oil from Chamaecyparis obtuse

The sedative effect of inhaled terpinolene in mice and its structure-activity relationships

Anxiolytic-like effects of phytol: Possible involvement of GABAergic transmission

Central effects of citral, myrcene and limonene, constituents of essential oil chemotypes from Lippia alba (Mill

Ficus hispida bark extract prevents nociception, inflammation, and CNS stimulation in experimental animal model

GABA A receptor modulation by terpenoids from Sideritis extracts

Drugs which stimulate or facilitate central GABAergic transmission interact synergistically with delta-9-tetrahydrocannabinol to produce marked catalepsy in mice

0: Web-service for prediction and optimization of chemical ADMET properties

SLOWER: A performance model for Exascale computing

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

Ab initio calculation of vibrational absorption and circular dichroism spectra using density functional force fields

Exploring the ligand efficacy of cannabinoid receptor 1 (CB1) using molecular dynamics simulations

Cannabidiol binding and negative allosteric modulation at the cannabinoid type 1 receptor in the presence of delta-9-tetrahydrocannabinol: An In Silico study

The PyMOL Molecular Graphics System

Computed atlas of surface topography of proteins with structural and topographical mapping of functionally annotated residues

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

Assignment of protonation states in proteins and ligands: Combining pK a prediction with hydrogen bonding network optimization

The amber lipid force field

Fast empirical pKa prediction by Ewald summation

New ways to boost molecular dynamics simulations

Increasing the precision of comparative models with YASARA NOVA-A selfparameterizing force field

A molecular modeling approach to identify effective antiviral phytochemicals against the main protease of SARS-CoV-2

Multivariate Calibration

Principal component analysis

Investigating the binding affinity, interaction, and structure-activity-relationship of 76 prescription antiviral drugs targeting RdRp and Mpro of SARS-CoV-2