Molecular Docking and Molecular Dynamics Simulation Studies of Triterpenes from Vernonia patula with the Cannabinoid Type 1 Receptor

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.


Introduction
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.

Molecular Docking
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.

Total Potential Energy Calculations
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.

Principal Components Analysis of MD Simulations
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).

Stability Analysis
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.

Residue Mobility Analysis
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).

MM-PBSA Binding Free Energy Analysis
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.

Discussion
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].

Plant Collection and Extraction
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).

Prediction of Pharmacokinetic and Drug-Likeness Properties
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].

Protein Optimization
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.

Determination of Ligand-Protein Binding Affinities and Non-Bonding Interactions
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.

Total Potential Energy Calculations
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 of MD Simulation Data
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.

Stability and Residue Mobility Analyses
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.

MM-PBSA Binding Free Energy Calculations
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.