A molecular informatics and in-vitro approach to evaluate the HMG-CoA reductase inhibitory efficacy of monoterpenes, carvacrol and geraniol

Abstract

Currently, statins, the β-hydroxy-β-methyl-glutaryl-CoA reductase (HMG-R) inhibitors, are widely used to lower cholesterol, nevertheless, they have several side effects. Consequently, the present study is designed to unravel the cardioprotective role of selected natural monoterpenoids (carvacrol and geraniol) via in-vitro targeting and molecular informatics study of HMG-R. Computational molecular informatics study revealed that carvacrol and geraniol efficiently occupies the catalytic site of HMG-R with the binding affinity (ΔG) of −4.60, and −1.99 Kcal/mol, respectively, and molecular mechanical-generalized Born surface area (MM-GBSA) free binding energy was depicted as −17.05 and −29.48 Kcal/mol, respectively. Further, molecular dynamics simulation was carried out for 100 ns. Carvacrol and geraniol potentially and competitively inhibit the in-vitro HMG-R activity with an IC₅₀ value of 78.23 ± 2.21 µM, and 72.91 ± 2.92 µM, respectively. Thus, both carvacrol and geraniol exhibited significant anti-hypercholesterolemic activity while the molecular simulation studies depicted that the GR complex showed better stability than the carvacrol complex.

KEYWORDS:

• Atherosclerosis
• hypercholesterolaemia
• monoterpenes
• HMG-R inhibition
• enzyme kinetics
• molecular docking

1. Introduction

Amid other progressive and non-communicable diseases; atherosclerotic cardiovascular diseases (ASCVDs) are the major universal health concern as well as the foremost reason for mortalities. However, the average number of ASCVD cases are rising globally every day [1,2]. Several cardiac health-associated risk factors i.e. hypertension, elevated levels of total cholesterol (TC), triglycerides (TG), lipoproteins and inflammation influenced the ASCVD severity [3]. Despite the existence of other risk factors, raised low-density lipoprotein (LDL-C) is still a major cause of cardiac health conditions [2]. The altered LDL-C (oxidized-LDL and glycated LDL) is more efficient in the progression of atherosclerotic plaque than unmodified LDL-C and about 70% of circulatory LDL is cleared by hepatic LDL receptor (LDL-R), and their transcriptional activators, sterol regulatory element binding proteins (SREBPs) [4,5].

Moreover, hepatic 3-hydroxy-3-methyl-glutaryl-CoA reductase (HMG-R) is the rate-limiting enzyme of cholesterol homeostasis, the foremost target for ASCVD management [3,6,7]. Presently, the statins are most widely prescribed inhibitors of HMG-R. However, statins have several side effects like statin intolerance, the most common problem in traditional preclinical practice [1,8]. Even though the mortality rate has dropped, more than 70% of ASCVD events are occurring regardless of statin therapeutic strategies. The adverse event caused by the statins indicates an urgent need for novel and safer drugs for cholesterol-lowering via targeting the HMG-R. Natural compounds, specially, plant-derived secondary metabolites received a lot of attention due to their medicinal properties including their less toxicity and cost efficiency. In this milieu, various plant extracts and phytochemicals have been screened for their hypolipidemic activity via targeting HMG-R [7,9–12].

Moreover, essential oils extracted from aromatic plants are a rich source of bioactive metabolites such as carvacrol (CR) and geraniol (GR), commonly referred to as monoterpenes (Figure 1). CR (C₁₀H₁₄O), a monoterpenoid phenol mainly found in the Carum copticum L., Origanum vulgare L., Zataria multiflora Boiss., Thymus vulgaris L., Nigella sativa L., Lepidium flavum Torr. [13–16]. The Federal Drug Administration (FDA) and the Council of Europe have approved its use as an antimicrobial, and flavouring agent in foods [17–20]. Several pharmacological properties of CR are well known like antiviral, antibacterial, anticancer, anti-inflammatory, antidiabetic and insecticides [21–37].

Figure 1. 2D-Schematic chemical structure of Geraniol and Carvacrol.

On the other hand, GR (C₁₀H₁₈O), a monoterpene alcohol, is majorly found in Monarda fistulosa, Ninde, Rose, Citronella and Palmarosa oil [38–42]. GR has distinct pharmacological properties such as antioxidative, anti-inflammatory, cardioprotective and neuroprotective [43–45]. Based on the above discussion, the present study is designed to unravel the cardioprotective role of selected natural monoterpenoids (carvacrol and geraniol) via molecular informatics study and in-vitro targeting of HMG-R.

2. Material and methods

2.1. Chemicals

Purified monoterpenoids i.e. CR, GR and HMG-R activity assay kit was procured from Sigma Aldrich, USA. Other chemicals were procured from HiMedia Lab. Pvt. Ltd., India and are of analytical grade.

2.2. Absorption, distribution, metabolism and excretion (ADME) profiling of monoterpenes, carvacrol and geraniol

The physiochemical and pharmacokinetic property of selected ligands was evaluated by using a web server preADME (https://preadmet.qsarhub.com/) [46]. On the basis of whether or not they complied with Lipinski's rule of five, Egan, Veber and Pians rules, the ligands were further examined for their drug-like qualities [47].

2.3. Predicted toxicity of the monoterpenes, carvacrol and geraniol

Toxicity analysis was completed via using the ProTox-II (https://tox-new.charite.de/protox_II/index.php?site = compound_input), a web server that forecasts the toxicity of minor compounds. It offers plentiful details of the compounds including LD₅₀, Carcinogenicity, Immunotoxicity, Mutagenicity Cytotoxicity as well as most important Hepatotoxicity [48].

2.4. Molecular docking simulation study platform

The in-silico investigation and visualizations were accomplished on a personal Lenovo machine with Intel (R) Core-i5, 8th gen processor equipped dedicated 2GB NVIDIA GeForce MX250 graphics card. The following software was utilized for this study GraphPad prism 8, Schrodinger’s (v 12.8) Maestro, and Desmond.

2.5. Ligand’s structure retrieval and preparation

The 3D structure of the CR and GR was retrieved from the PubChem (http://pubchem.ncbi.nlm.nih.gov) database. The retrieved 3D structure.sdf files were incorporated in the Maestro workspace and all the ligands were minimized with the help of the LigPrep module by applying the OPLS_2005 force field for the molecular docking studies. The Epik module was utilized for possible ionization and tautomer of selected structures originated at pH 7.0 +/− 2.0 [49,50].

2.6. Retrieval and preparation of human HMG-R for docking

The 3D structure of the Human HMG-R (PDB-ID: 1HW8, resolution of 2.10Å) [51], protein was downloaded from the Protein Data Bank http://www.rcsb.org/ [52]. We prepared our docking target protein was prepared using the protein preparation wizard's of Schrödinger's Maestro module. The assignment of hydrogen bonds (H-bonds), bond orders, addition of hydrogen atoms, optimization and minimization, removal of water beyond five angstroms from the het group, and conversion of selenomethionines to methionines were important steps in the preparation. Further, pH 7.0 +/− 2.0 was used to fill the loop and then pre-processed with the protein. To refine the targeted protein, the OPLS_2005 force field was applied and converging complete heavy atoms to a root mean square deviation (RMSD) of 0.30 with the help of the protein preparation module Maestro v12.8. Using the SiteMap tool, a highly potential binding pocket of ligands on proteins was identified.

2.7. Receptor grid generation

The sitemap tool in the Maestro v12.8 module was used by means of Goodford's GRID algorithm. This tool provides the available regions in the target binding site for ligands [53,54]. The five-ligand binding site was identified by Sitemap, for generation of the receptor grid, we elect the highest site score and size of the volume to make a grid on the binding site and chosen a van der Waal's radius scaling of 1.0 with 0.25 partial charges cut-off by default. The retrieved receptor (1HW8) has a crystallographic native ligand in the identical binding site. Since both SiteMap's generated active pocket and co-crystalized ligand with IHW8 showed similar active sites. The receptor grid generation was made by using the “pick to identify the ligand molecule” method in the Glide module [55,56].

2.8. Molecular docking studies of carvacrol and geraniol

We used the Maestro’s Glide module to access the binding pattern of CR and GR compounds against the target protein. The previously prepared receptor grid file and prepared ligand file were input in the Glide module of protein–ligand docking and by leaving all the fields in default mode the job was submitted. Post molecular docking, the ligand binding position examination and construction of data were done with the help of Maestro 12.8 [57,58].

2.9. Analysis of MM-GBSA score

The molecular mechanical-generalized Born surface area (MM-GBSA) score for the estimation of binding free energies of selected ligands with the receptor proteins by using the Prime MM-GBSA module of Schrodinger. The OPLS_2005 force field and variable dielectric surface generalized born (VSGB) solvation model were utilized in this process [59,60]. Finally, the binding free energy of a protein–ligand complex was computed using the following equation: $\begin{array}{rl}r\mathrm{\Delta }G{B}_{Bind}& =\mathrm{\Delta }{G}_{Coulomb}+\mathrm{\Delta }{G}_{vdW}+\mathrm{\Delta }{G}_{Covalent}\\ & +\mathrm{\Delta }{G}_{H\text{-}bond}+\mathrm{\Delta }{G}_{Sol\mathrm{\_}Lipo}+\mathrm{\Delta }{G}_{Solv\mathrm{\_}GB}\\ & +\mathrm{\Delta }{G}_{Packing}+\mathrm{\Delta }{G}_{Self\text{-}contact}\end{array}$Where, ΔG_Bind, ΔG_Coulomb, ΔG_vdW, ΔG_Covalent, ΔG_H-bond, ΔG_{Sol_Lipo}, ΔG_{Solv_GB}, ΔG_Packing and ΔG_Self-contact represent total free energy of protein–ligand complex, Coulombic free energy, van der Waals’ interactions free energy, covalent bonds free energy, hydrogen bonds free energy, the free energy of the surface area, Solvation free energy, packing free energy and free energy of self-contact respectively. MM-GBSA has been shown to perform better in scoring congeneric ligands than Glide [61,62].

2.10. Molecular dynamic simulation analysis

MD simulation of CR and GR were performed in complex with its respective targeted enzyme, HMG-R, using “Desmond (Schrodinger 2020, LLC, NY, USA)” as described earlier [62,63]. The receptor–ligand complex was preprocessed using Protein Preparation Wizard of Maestro, which included complex optimization and minimization. All the systems were prepared using the System Builder tool. TIP3P, a solvent model with an orthorhombic box (10 Å), was chosen. In the simulation, the OPLS_2005 force field was used [64]. To make the models neutral, counter ions were introduced. To mimic physiological conditions, 0.15 M sodium chloride (NaCl) was added. The NPT ensemble with 300 K temperature and 1 atm pressure was chosen for the entire simulation. The models were relaxed before the simulation. The trajectories were saved for examination after every 100 ps. The parameters such as RMSD, root mean square fluctuation (RMSF), Radius of gyration (Rg), Solvent Accessible Surface Area (SASA) and protein–ligand interactions were analyzed to establish the stability of protein–ligand complexes. The Principal Component Analysis (PCA) was analyzed by using Bio3D package of R [65]. PCA is calculated using a script written in R [66].

2.11. In-vitro HMG-R inhibitory activity

The HMG-R assay kit for in-vitro study includes human HMG-R catalytic domain (recombinant GST fusion protein expressed in E. coli; 0.52–0.85 mg/mL). The cardioprotective potential of CR and GR on HMG-R catalytic activity was measured as described earlier [7].

2.12. Enzyme kinetics studies

The enzyme kinetic experiments were conducted to evaluate the mechanistic manner of inhibition of HMG-R via CR and GR. Different concentrations of the substrate, HMG-CoA i.e. 66.66, 133.33 and 266.66 µM were added to monitor the HMG-R activity using distinct amounts of CR and GR (100, 200 and 300 µM) as well as pravastatin (0, 50 and 100 nM). Further, the double-reciprocal Lineweaver–Burk plot was designed to unknot the mode of inhibition as described previously [10].

2.13. Statistical analysis

All the in-vitro biochemical assays were carried out in triplicate and data was represented as mean ± SEM. The GraphPad Prism version 4.02 (GraphPad Software, San Di-ego, USA) was used for the statistical analysis.

3. Results and discussion

3.1. Pharmacokinetic profiling and drug-likeness properties of carvacrol and geraniol

Targeting the pathway responsible for the endogenous production of cholesterol is the most preferred course of action against hypercholesterolaemia. Since, HMG-R is the rate-limiting enzyme of cholesterol homeostasis [5] and statins are the most widely used inhibitors of HMG-R but their long-term use causes adverse outcomes as reported by several reports [8,67–70]. In this context, various plant extracts and phytochemicals have been screened for their hypolipidemic activity via targeting HMG-R [1,10]. Natural compounds specially plant-derived secondary metabolites received a lot of attention due to their medicinal properties including their fewer toxicity and cost-effectiveness [71].

We initially executed the pharmacokinetic profile i.e. absorption, distribution, metabolism, excretion and toxicology (ADMET), drug-likeness and predicts the human intestinal absorption, blood–brain barrier (BBB) potential, skin permeability, plasma protein binding (PPB) potential of a chemical entity [46]. while the physiochemical study gives the pharmacological feature of drug candidate such as provide ranges for the lipophilicity, water solubility, topological surface area (TPSA) as well as general feature of drug candidate i.e. molecular weight, number of atoms, number of rotatable bonds, number of hydrogen bonds donor and acceptor.

Since, drug-likeness (i.e. Lipinski’s rule of five), describes the physiochemical ranges of the bioactive compounds to being an oral drug, our results demonstrated that both CR and GR have acceptable scores of Lipinski’s rule of five. On the other hand, CR and GR does not show any violation of Veber and Egan rules and no alerts was predicted against PAINS (pan assay interference structure) rules, that identifies of potentially problematic fragments, as well as one violation (Molecular weight ≤160) was observed in Goshe rule. These properties of CR and GR indicated that they fulfil the criteria of drug-likeness and pharmacological properties (Table 1).

Table 1. Chemical and Drug likeness properties of test compounds and standard drug.

Compound Name	Carvacrol	Geraniol	*Pravastatin	*Atorvastatin
PubChem ID	10364	637566	54687	60823
M. wt.	150.22	154.25	424.53	558.65
Log P	2.24	2.52	3.01	3.81
HBA	1	1	7	6
HBD	1	1	4	4
Rotatable Bonds	1	4	11	13
Violation	0	0	0	1
Aq. Solubility (mg/L)	180.95	665.54	405.17	0.11
BBB	6.39	6.74	0.06	0.67
CaCo2	38.01	8.76	18.94	21.71
HIA	100.00	100.00	77.69	94.65
MDCK	87.33	271.03	0.08	0.04
Skin permeability (logKp)	−1.05	−1.06	−3.08	−2.42
Plasma protein Binding	100.00	100.00	85.81	90.34
CYT. p450 2D6 inhibition	Non	Non	Non	Non
YT. p450 2D6 substrate	Non	Non	Non	Non

*Standard drug.

Further, our in-silico investigation depicted that the CR and GR showed considerable aq. solubility 180.95 and 665.54 mg/L, respectively (Table 1). Moreover, the ligands with BBB scores > 2.0, between 2.0–0.1 and <0.1 show high, modest and low absorption across BBB, respectively [72,73]. While pravastatin and atorvastatin score low BBB absorption. Additionally, the valuation of colorectal adenocarcinoma (Caco-2) and Madin-Darby canine kidney (MDCK) cell penetrability has been a well-known measure for the drug development procedure [73,74]. Results exhibited that CR and GR score Caco-2 permeability with 38.01 and 8.76 nm/s, while 87.33 and 271.03 nm/s were observed in MDCK permeability, respectively. This indicates that CR and GR have the potential capability to pervade across Caco-2 and MDCK cells and are well supported by known permeability ranges across these cellular models [75].

Our selected compounds also showed 100% PPB while, pravastatin, and atorvastatin score 86 and 90% respectively. Further, the predicted human intestinal absorption (HIA) for CR and GR was observed as 100%. Whereas, pravastatin and atorvastatin exhibit 77.95 and 94.65%, respectively. The enhanced HIA exhibits the improved bioavailability of these compounds which can enhance their potent pharmacological activities. In the same way, our results depicted that CR, GR, pravastatin and atorvastatin exhibit non-reactivity toward CYP2D6 inhibitor and substrate which indicates comparatively faster metabolism, instantaneous beneficial actions as well as clearance from the body. These results are well acceptable by the available reports [76,77].

3.2. Toxicity assessment of carvacrol and geraniol compounds

The selected monoterpenes; CR and GR fall under the class 4 and 5 with LD₅₀ of 810 and 2100 mg/kg, respectively. That means these compounds were safe and inactive against all the tested toxicity parameters whereas, pravastatin was active against the immunotoxicity, carcinogenicity and predicted toxicity class was 6 as well as LD₅₀ was 8939 mg/kg. While, atorvastatin showed hepatotoxicity and falls under the toxicity class of 5 and 5000 mg/kg LD₅₀ was predicted (Table 2). To the sum of all, the selected test compounds are predicted as non-toxic and safe for further carrying out in the drug development programme.

Table 2. Toxicity assessment of selected compounds.

Compound Name	LD50 (mg/kg)	Class	Hepatotoxicity	Carcinogenicity	Immunotoxicity	Mutagenicity	Cytotoxicity
Carvacrol	810	4	Inactive	Inactive	Inactive	Inactive	Inactive
Geraniol	2100	5	Inactive	Inactive	Inactive	Inactive	Inactive
*Pravastatin	8939	6	Inactive	Active	Active	Inactive	Inactive
*Atorvastatin	5000	5	Active	Inactive	Inactive	Inactive	Inactive

*Standard drug.

3.3. In-silico molecular docking revealed that carvacrol and geraniol efficiently binds to the HMG-R active pocket and MM-GBSA analysis

The computer-aided molecular docking studies were executed to determine the binding affinity of CR and GR within the catalytic site of HMG-R. Our results illustrated that CR and GR occupied the active pocket of HMG-R with the binding energies (ΔG) of −4.60 and −1.99 Kcal/mol, respectively (Table 3). CR was comparable and showed better ΔG against pravastatin and atorvastatin (−4.21 and −3.26 Kcal/mol, respectively). The result for the standard drug was well corroborated with the previously published studies [7,78]. The involvement of 8 amino acid residues Arg⁵⁹⁰, Met⁶⁵⁵, Met⁶⁵⁷, Asn⁶⁵⁸, Asp⁶⁹⁰, Lys⁶⁹¹, Asp⁷⁶⁷ and Gln⁷⁷⁰ stabilized the CR-HMG-R−complex. The stabilized complex has one conventional H-bond with Asp⁶⁹⁰ residue along with one pi-cation bond and 2 hydrophobic interactions along with two negative and two positive charge interactions with the target protein (Figure 2) (Table 3). The pi-cation bond interaction is a crucial participant in the architecture and stability of protein [79–81]. Compared to an H-bond, pi-cation bonds are stronger energetically [82]. In studies, it has been demonstrated that the pi-cation interaction promotes recognition between proteins and ligands and that it can also be regarded as a valuable predictor of drug–receptor interactions [83]. Whereas, GR occupied the active pocket of HMG-R with the 8 identical amino acid residues i.e. Arg⁵⁹⁰, Met⁶⁵⁵, Met⁶⁵⁷, Asn⁶⁵⁸, Asp⁶⁹⁰, Lys⁶⁹¹, Asp⁷⁶⁷, Gly⁸⁰⁷. These interactions incorporated by two conventional H-bonds with Arg⁵⁹⁰ and Asp⁶⁹⁰ along with two similar hydrophobic interaction (Met⁶⁵⁵ and Met⁶⁵⁷), and two negative and two positive charged interaction (Figure 2) (Table 3).

Figure 2. 2D-Schematic representation of Carvacrol and Geraniol with the target protein HMG-R. (A) Carvacrol interaction with HMG-R in the binding pocket visualized through the 2D structure. (B) The 2D structure interaction representation of Geraniol with the active pocket amino residues of HMG-R.

Table 3. Number of residues stabilizing the carvacrol, geraniol and the standard drug atorvastatin and pravastatin-HMG-R-complex.

Compound name (PubChem ID)	Docking Score (Kcal/mol)	Glide score	H-bond interaction	Residues involved	Common amino acid involved in each ligand interaction
Carvacrol (10364)	–4.60	−4.60	Asp^690	Arg^590, Met^655, Met^657, Asn^658, Asp^690, Lys^691, Asp^767, Gln^770.	Arg^590, Met^657, Asn^658, Asp^690, Lys^691.
Geraniol (637566)	–1.99	−1.99	Arg^590, Asp^690	Arg^590, Met^655, Met^657, Asn^658, Asp^690, Lys^691, Asp^767, Gly^807.
HMG-CoA (445127)	−4.62	−5. 95	Glu^528, Arg^590, Asn^658, Ser^684, Asp^690.	Tyr^479, Lys^480, Gly^524, Ala^525, Cys^526, Glu^528, Arg^590, Met^657, Asn^658, Ser^684, Asn^686, Cys^688, Asp^690, Lys^691, Lys^692, Ala^769.
*Pravastatin (54687)	–4.21	−4.21	Ala^525, Asp^690, Lys^691.	Ala^525, Cys^526, Cys^527, Glu^528, Arg^590, Met^655, Met^657, Asn^658, Asp^690, Lys^691, Lys^692, Asp^767, Gly^807, Gly^808.
*Atorvastatin (60823)	–3.26	−3.26	Lys^691, Lys^692, Arg^590.	Arg^590, Met^657, Asn^658, Ser^661, Glu^665, Val^682, Val^683, Ser^684, Asn^686, Asp^690, Lys^691, Lys^692.

*Standard drug.

Furthermore, the molecular binding pattern of the substrate (HMG-CoA) and the standard drug pravastatin as well as atorvastatin were also analyzed to compare the binding efficiency of CR and GR with the active pocket of HMG-R. The HMG-CoA interaction with HMG-R showed the highest binding affinity (−5.95 Kcal/mol) in comparison to our tested compounds. It stabilizes the HMG-CoA−HMG-R complex with the formation of 16 amino acid residues i.e. Tyr⁴⁷⁹, Lys⁴⁸⁰, Gly⁵²⁴, Ala⁵²⁵, Cys⁵²⁶, Glu⁵²⁸, Arg⁵⁹⁰, Met⁶⁵⁷, Asn⁶⁵⁸, Ser⁶⁸⁴, Asn⁶⁸⁶, Cys⁶⁸⁸, Asp⁶⁹⁰, Lys⁶⁹¹, Lys⁶⁹² and Ala⁷⁶⁹ along with five H-bonds, one water mediated H-bond with Asp⁶⁹⁰ as well as six hydrophobic and three polar interactions were formed.

Whereas, pravastatin occupied the active site of HMG-R with a ΔG of −4.21 Kcal/mol, this complex was stabilized by the participation of 14 amino acid residues (Ala⁵²⁵, Cys⁵²⁶, Cys⁵²⁷, Glu⁵²⁸, Arg⁵⁹⁰, Met⁶⁵⁵, Met⁶⁵⁷, Asn⁶⁵⁸, Asp⁶⁹⁰, Lys⁶⁹¹, Lys⁶⁹², Asp⁷⁶⁷, Gly⁸⁰⁷ and Gly⁸⁰⁸) that lied within the active pocket of the enzyme, 4 conventional H-bonds, five hydrophobic, three positive, three negative and one polar interaction. On the other hand, atorvastatin−HMG-R−complex was formed by 12 amino acids (Arg⁵⁹⁰, Met⁶⁵⁷, Asn⁶⁵⁸, Ser⁶⁶¹, Glu⁶⁶⁵, Val⁶⁸², Val⁶⁸³, Ser⁶⁸⁴, Asn⁶⁸⁶, Asp⁶⁹⁰, Lys⁶⁹¹ and Lys⁶⁹²), this interaction stabilized by three H-bonds, three hydrophobic interaction, four polar, two positive and three negative interactions (Figure 3). Pravastatin and atorvastatin are quite bigger molecules than the CR and GR compounds that contain a greater number of carbons, hydrogen and oxygen, which provides the feasibility of the interaction with the target. Whereas, the binding affinity of standard molecules was influenced by hydrogen bonding and other interactions. The bonding pattern with the amino acid and binding energy also depends upon the bond type, number of rotatable, hydrogen donor, acceptor bonds and topological surface area. Each molecule has its own physiochemical properties that are responsible for their interaction with the targeted protein [84,85].

Figure 3. 2D-Schematic representation of pravastatin, atorvastatin and HMG-CoA with the target protein HMG-R. (A) 2D interaction representation of pravastatin with the active pocket of HMG-R. (B) Representation of atorvastatin with HMG-R (C) HMG-CoA interaction with HMG-R in the binding site is visualized through the 2D structure.

It has been found that Prime MM-GBSA investigation improves the binding energy calculations as compared with molecular docking energies. As a result, MM-GBSA analysis discloses that the ligands bind more strongly to the receptors. Our study unveiled that the MM-GBSA scores of CR and GR complexes are −17.05 and −29.48 Kcal/mol respectively (Table 4). Instead, the MM-GBSA score for the HMG-CoA-HMG-R complex is −18.70 Kcal/mol. Moreover, the binding free energy of pravastatin and atorvastatin of −39.47 and −35.89 Kcal/mol, respectively. Thus, GR showed stronger binding free energy which is closer to atorvastatin and pravastatin, while CR exhibited a lower MM-GBSA score among others.

Table 4. Components of free energy of target protein (HMG-R) with the identified inhibitor using Prime/MMGBSA approach. All the energies are reported in Kcal/mol⁻¹.

Compound name	Carvacrol	Geraniol	HMG-CoA	*Pravastatin	*Atorvastatin
PubChem ID (CID)	10364	637566	445127	54687	60823
Glide Score (Kal/mol)	–4.60	–1.99	–4.62	–4.21	–3.26
MMGBSA_ΔG_Bind	–17.0542	–29.4603	–18.6961	–39.4696	–35.8869
MMGBSA_ΔG_Bind_Coulomb	–6.21268	–15.3648	–108.965	–38.0835	0.616678
MMGBSA_ΔG_Bind_Covalent	0.791871	0.134129	4.157413	2.08821	–1.03162
MMGBSA_ΔG_Bind_H-bond	–0.90141	–1.97611	–13.5219	–2.9454	–2.17027
MMGBSA_ΔG_Bind_Lipo	–7.49686	–16.3907	7.215147	–14.2742	–15.9854
MMGBSA_ΔG_Bind_Solv_GB	7.175254	18.78108	111.3575	40.261	13.93434
MMGBSA_ΔG_Bind_vdW	–10.4104	–14.6439	–18.9387	–26.0755	–31.2506
Lig_Strain_Energy	0.764257	1.238485	38.63314	2.771479	9.476791

*Standard drug.

3.4. Molecular dynamics simulation study of carvacrol and geraniol

It is a critical approach to discover the stability of protein–ligand complexes and analyze the intermolecular interactions. The molecular dynamic simulation (MDS) provides insight into the intermolecular interactions present in the protein–ligand complexes within the given reference time interval. In our case, the complexes of the selected compounds with the target, (HMG-R), were analyzed by using the MDS approach to envisage the complex stability, conformational changes and intermolecular interactions at the atomic level of the selected CR and GR complexes, against the time interval of 100 ns. MDS trajectories of each complex were then analyzed using the simulation interaction diagram (SID) module of the Desmond interface and simulation results were analyzed including RMSD, root mean square fluctuation (RMSF) and protein–ligand interaction mapping of the complexes. The RMSD parameter is used to measure the average distance generated due to the removal of the selected atoms in the complex during the simulation interval for a specific time interval concerning the reference time frame. The RMSD of a docked complex of selected compounds with the target, HMG-R, were studied for the C-α backbone, side chain and heavy atoms. The RMSD of the protein fit ligand from all time frames concerning the reference time is calculated.

3.4.1. RMSD of protein

The RMSD values between 1 and 3 $\text{Å}$ within the reference protein structure are acceptable, but the larger value indicates significant conformational changes or fluctuations in the protein–ligand complex. The more considerable conformational changes in protein refer to the unstable and nonequilibrium state of the system during the simulation interval. The RMSD analysis of CR and GR shows that the C-α atoms of the target protein, exhibit some oscillation in RMSD values in both complexes. The RMSD values of both complexes were mostly around 0−6 $\text{Å}$, which is higher than the acceptable range. In CR−HMG-R complex protein RMSD were stable from starting to 18 ns, after that (19−30 ns) showed some fluctuation in RMSD values, and later (30−100 ns) the RMSD values of protein stabilized in between 4 and 6 $\text{Å}$. Similarly, the GR−HMG-R complex protein RMSD value mostly lies between 3 and 5 $\text{Å}$ (Figure 4).

Figure 4. The Root Mean Square Deviation (RMSD) plot of (A); carvacrol and (B); geraniol with HMG-R complex. The RMSD is used to measure the average change in displacement of a selection of atoms for a particular frame with respect to a reference frame. It is calculated for all frames in the trajectory. (C); Principal Component Analysis eigenvalue plotted versus the percentage of variance CR−HMG. The varying areas are displayed in three separate sections. Variations in PC1, PC2 and PC3 add up to 57.89%, 12.51% and 5.77%, respectively. (D); Principal Component Analysis eigenvalue plotted versus the percentage of variance (CR−HMG-R). The varying areas are displayed in three separate sections. Variations in PC1, PC2 and PC3 add up to 45.4%, 16.51% and 8.08%, respectively.

3.4.2. RMSD of ligand

After analyzing the RMSD data obtained from the protein fit ligands, we found that all the docked proteins with the selected compounds were favoured by RMSD values <3 Å. The CR and GR compounds exhibited higher RMSD values throughout the simulation of 100 ns. The CR showed higher RMSD values which extend up to 5−105 Å at some intervals when compared to favoured RMSD value, which showed structural modification and less stabilization in ligands throughout the simulation (Figure 4). A similar pattern was also observed in GR−HMG-R-complex, where the ligand RMSD values stretched over the favoured region at the initial simulation time, and it deviated between 16−20 in the RMSD values throughout the simulation period.

3.4.3. Principle component analysis

Principal Component Analysis (PCA) is a statistical technique that can be used to identify the essential dynamics of proteins, which describe the change in molecular structure or conformation over time [86]. The graph of eigenvalues (protein) against eigenvector index (eigenmode) for the first 20 modes of motion (CR−HMG-R and GR−HMG-R) (Figure 4C,D, respectively) represents stability. The eigenvalues depicted the hyperspace eigenvector fluctuations. In simulations, the eigenvectors with the highest eigenvalues control the total mobility of the target protein. The CR−HMG-R and GR−HMG-R, top 5 eigenvectors in utilized systems displayed dominant movements and had higher eigenvalues (57.9−87.1%), and (45.4−84.4%) respectively, than the other eigenvectors. All changes were recorded and plotted on three PCs (PC1, PC2 and PC3). PC1 clusters had the largest variability (57.89% and 45.4% respectively), PC2 showed variability (12.51% and 16.51% respectively) and PC3 had the lowest variability (5.77%, and 8.08% respectively). Because of its low variability, PC3 has a more compact structure than PC1 and PC2, as well as was thought to be a more stabilized protein–ligand binding complex. The simple clustering in the PC subspace revealed conformational variations across all the groups. Blue signifies the most significant mobility; white shows intermediate movement; and red suggests less flexibility [87].

3.4.4. Root mean square fluctuation (RMSF)

Typically, we observe that the N and C terminal residues fluctuate more from the rest of the protein parts. Such fluctuations can be ignored during the simulation and RMSF analysis. In the complex of CR−HMG-R, all the residues of HMG-R i.e. from Ile⁵³¹ to Val⁸⁶³ exhibited about <3 Å fluctuation in RMSF values. While, a higher fluctuation (3−6.3 Å) was observed between Leu⁴⁶² to Val⁵³⁰ throughout the 100 ns simulation time. While the amino acids Gln⁴⁶⁹, Asn⁴⁷², His⁴⁷⁵, Ile⁴⁷⁶, Glu⁴⁸², Asp⁵¹⁶, Tyr⁵¹⁷, Ser⁵²⁰, Val⁵²², His⁵²⁴, Ala⁵²⁵, Cys⁵²⁶ and Cys⁵²⁷ interacted with ligand but showed high fluctuation (5.54, 5.91, 6.3, 6.09, 5.76, 3.60, 3.05, 4.53, 4.85, 5.31, 4.51, 4.53 and 3.81Å respectively). On the other hand, most of the protein–ligand contacts were observed in the favoured RMSF region. As shown in Figure 5, in the GR−HMG-R-complex, the fluctuation in most of the residues occurs between 1 and 3.5 Å. The highest fluctuation was observed in the loop region of HMG-R 7.78 Å, while the peak showed that the following amino acid residues i.e. Ser⁵⁰⁰, Lys⁵⁰¹, Lys⁵⁰², Leu⁵⁰³, Ser⁵⁰⁴, Glu⁵⁰⁵, Pro⁵⁰⁶, Ser⁵⁰⁷, Ser⁵⁰⁸ with fluctuation of 4.31, 5.71, 6.01, 6.17, 7.78, 6.32, 5.07, 4.42 and 4.34Å, respectively was observed. Comparing the RMSF of all the complexes, it has been depicted that GR exhibited more stability than the CR complex through the simulation. The GR can be considered a stable compound as compared to CR because carvacrol has a higher RMSF value.

Figure 5. RMSF plot of representation of (A); carvacrol, and (B); geraniol with HMG-R. RMSF plot defines characterizing local fluctuations along the protein chain. In the above plot, peaks represent regions of the protein that fluctuate most during the simulation. Typically, N- and C-terminal fluctuate more than any other part of the protein. Secondary structure elements like alpha helices and beta strands are usually more rigid than the unstructured part of the protein, and thus fluctuate less than the loop regions. Dependency of (C and E); radius of gyration (Rg) of CR−HMG-R and GR−HMG-R, respectively, and (D and F); solvent accessible surface area (SASA) of target protein of HMG-R in the presence of CR and GR, respectively on simulation time.

3.4.5. The radius of gyration (Rg) and solvent-accessible surface area (SASA)

The radius of gyration measures the protein–ligand complex’s compactness, and the solvent-accessible surface area determines its exposure to solvent molecules. Each of these characteristics sheds light on the protein–ligand complex’s stability throughout the simulation in this study, the Rg of HMG-R was determined in the presence of CR and GR during 0–100 ns, as shown in Figure 5C,E. The Rg of CR-HMG-R and GR-HMG-R complexes varied in the range of 23.4−26.0 Å and 24.8−26.4 Å respectively, with an average value of 24.5 ± 0.67 Å and 25.89 ± 0.37 Å, respectively. The Rg values indicate how compact or extended the complexes are. Smaller Rg values typically suggest a more compact or folded structure, while larger Rg values suggest a more extended or unfolded structure. Similarly, the SASA of HMG-R was determined in the presence of CR and GR during 0–100 ns (Figure 5D,F). The SASA of CR-HMG-R and GR-HMG-R complexes fluctuated in the range of 18,500−20,500 Ų and 18,300−20,200 Ų, respectively, with an average value of 19,362.92 ± 469.00 Ų and 19,744.82 ± 253.38 Ų, respectively. Higher SASA values often indicate that more of the protein's surface is exposed to solvent molecules, which might be indicative of increased conformational flexibility or greater solvent accessibility. Conversely, lower SASA values can suggest a more compact or stable protein–ligand complex. The above findings display the increased conformational flexibility or greater solvent accessibility of protein–ligand complexes by showing that fluctuations in Rg and SASA of the protein of interest in the presence of CR and GR diverge significantly.

3.4.6. Protein–ligand contacts

Numerous types of contacts between protein and ligands, including H-bonding have an essential role in drug metabolism, absorption and drug specificity. Besides H-bond interactions, the π–π interactions, ionic or polar interactions, and water bridge contacts are also responsible for the stability of the protein–ligand complex. Furthermore, CR and GR compounds form many intermolecular contacts like H-bond interactions, hydrophobic interactions and water bridge interactions, thus exhibiting moderate stability throughout the simulation interval.

In the CR−HMG-R complex, there is the presence of ionic interaction established by Gly⁷⁷⁶ having an interaction fraction of 0.001. Besides H-bond, hydrophobic and water bridge interaction also formed during simulation, which might contribute towards the stability of the complex. The active site residues establish major H-bond interactions such as Lys⁷¹⁵, Glu⁷¹⁹ and Gly⁷⁴⁷ with interaction fractions of 0.011, 0.021 and 0.036, respectively. On the other hand, major hydrophobic interactions were exhibited by Val⁵³⁸, Lys⁷¹⁵, Pro⁵³⁵, Tyr⁵³³, Trp⁶⁹⁸, Val⁷¹⁶, Ile⁷⁷⁸, Pro⁸¹³ and Ala⁸¹⁶ amino acids; as well as water bridge interactions also formed during simulation (Figure 6).

Figure 6. The properties of both ligands (carvacrol and geraniol) during the 100 ns simulation. (A, B) The carvacrol and geraniol interaction with HMG-R. (B, D) The histogram of both ligands (carvacrol and geraniol) displayed bonds with amino acids of HMG-R.

In the GR−HMG-R complex, several H-bond, hydrophobic and water bridge interaction was formed during the simulation. The major H-bonds were formed in between the Asp⁵¹⁶, Tyr⁵¹⁷, Tyr⁵³³, Met⁵³⁴, Ile⁵³⁶, Val⁵³⁸ and Gln⁸¹⁴ with the interaction fraction of 0.074, 0.021, 0.021, 0.157, 0.043, 0.019 and 0.020 respectively. While, the major hydrophobic interaction was formed by similar amino acids which were involved in H-bond i.e. Tyr⁵¹⁷, Tyr⁵³³ Ile⁵³⁶, Val⁵³⁸, Val⁵⁴⁰, Ala⁵⁵⁶, Ile⁷⁶², Pro⁸¹³, Ala⁸¹⁶ and Met⁸²⁰ with the respective interaction fraction of 0.034, 0.118, 0.068, 0.017, 0.043, 0.037, 0.041, 0.133, 0.066 and 0.031. Similarly, a large number of water bridges with high interaction fractions were also exhibited between those above-mentioned amino acids such as, Asp⁵¹⁶, Tyr⁵¹⁷, Tyr⁵³³, Met⁵³⁵, Ile⁵³⁶, Val⁵³⁸, Gln⁸¹⁴ etc. (Figure 6). The overall results exhibited that various parameters described above like RMSD, RMSF and protein–ligand contacts were analyzed by the Desmond simulation interaction diagram module and we depicted that the GR complex showed better stability than the CR complex during the simulation interval of 100 ns.

3.5. Carvacrol and geraniol potentially inhibit the HMG-R activity and enzyme kinetics study

The above in-silico data was further validated through the in-vitro analysis. In order to validate the cardioprotective effect of CR and GR, in-vitro HMG-R inhibitory assay was performed following the standard protocol [10]. Our results demonstrated that CR and GR have potent inhibitory potential of HMG-R with an IC₅₀ value of 78.23 ± 2.21 µM, and 72.91 ± 2.92 µM, respectively (Figure 7A,B). However, the standard drug pravastatin showed lower IC₅₀ (58.12 ± 3.01 nM) than these monoterpenes. The spectrophotometric time scans of CR and GR (100−300 µM), and pravastatin (25−100 nM) suggested the direct interaction of these monoterpene drug candidates with HMG-R (Figure 8). It is obvious from our results that the IC₅₀ value of our tested compounds against HMG-R is quite higher but nevertheless, they are from natural sources [1]. The toxicity analysis already depicted that pravastatin and atorvastatin are active against carcinogenicity, immunotoxicity and hepatotoxicity (Table 2).

Figure 7. The bar graphs (upper panels, 7A) represent the % inhibition of in-vitro HMG-R activity in the presence or absence of carvacrol (IC₅₀ = 78.23 µM), geraniol (IC₅₀ = 72.91 µM) and pravastatin (IC₅₀ = 58.12 nM). The values of % inhibition are expressed as mean ± SEM from triplicate assays for each concentration. The Significance was defined as ***p < 0.001. The lower panel 7B represents the Lineweaver–Burk plot for the inhibition of HMG-R in the absence and presence of inhibitor at varying concentrations (100, 200 and 300 µM) and substrate at varying concentrations (i.e. 66.66, 133.33 and 266.66 µM). Note: The concentration of pravastatin is in nM.

Figure 8. The spectrophotometric time scans showing the progress of the reaction catalyzed by HMG-R either in the presence or absence of test compounds, (panel A): time scan for carvacrol; (panel B): for geraniol and, (panel C): for pravastatin at the different concentration.

In addition, to depict the mode of inhibition of HMG-R activity via these compounds, enzyme kinetics studies were conducted. The Lineweaver–Burk double reciprocal plot that is 1/V vs 1/S plotted, where V represents the velocity of enzyme activity (absorbance change) and S stands for the concentration of substrate. The competitive inhibitor (all trendlines intersect at the same position on the Y axis) shows no change in the V_max and an increase in the K_m value, and for the non-competitive inhibitor (all the trendlines intersect at the same position on the X axis); V_max decreases and K_m remain same. Our results demonstrated that presence of different concentrations of inhibitor (CR & GR), the K_m was increased and there is no significant change was observed in the V_max (Figure 7B). Thus, we concluded that both the monoterpenes CR and GR are competitive inhibitors of HMG-R. However, pravastatin is a well-established competitive inhibitor of HMG-R (Figure 7B) [8,88]. Therefore, we can say that CR and GR occupy the active site of the HMG-R via hydrophobic as well as hydrophilic interactions.

4. Conclusion

Our pharmacokinetic profiling and drug-likeness analysis demonstrated that these monoterpenes (CR and GR) followed Lipinski’s rule of five and were inactive against the toxicity assessment. However, pravastatin and atorvastatin were active against carcinogenicity and immunotoxicity. In addition, the molecular docking studies demonstrated that CR and GR actively occupied the catalytic site of the HMG-R which is comparable to the pravastatin. The MDS exhibited high RMSD and RMSF for both the complex of CR and GR. Due to the high RMSD and RMSF, we can say that both the complexes undergo some structural modification, but GR showed lower RMSD and RMSF values as compared to CR, which in turn depicted that GR complex showed better stability than CR complex. Furthermore, the in-vitro inhibition and kinetic study revealed that CR and GR have potent inhibitory activity of HMG-R with IC₅₀ values of 78.23 ± 2.21 µM, and 72.9 ± 2.92 µM, respectively. Additionally, our enzyme kinetics analysis depicted that both CR and GR competitively blocked HMG-R activity. Based on combined in-vitro and in-silico data, our results suggested that both CR and GR exhibited significant anti-hypercholesterolemic activity while the molecular simulation studies depicted that GR complex showed better stability than CR complex. Further, a broader in-vivo study is needed for the better understanding of the cardioprotective effect of these compounds.

Acknowledgments

The authors would like to thank the Deanship of Scientific Research at Majmaah University for supporting this work under project no. R-2023-886. 

Disclosure statement

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

Data availability statement

All the data supporting the findings are included in the article.