Virtual Screening of Sesquiterpenoid Pogostemon herba as Predicted Cyclooxygenase Inhibitor

Abstract

To analyze the structural features that dictate the selectivity of the two isoforms of the cyclooxygenase (COX), the three-dimensional structure of COX-1/COX-2 was assessed by means of binding energy calculation by way of virtual molecular dynamic simulations using ligand sesquiterpenoid Pogostemon herba. This study was conducted to investigate the molecular interaction between ligand alpha-bulnesene (CID94275), alpha-guaiene (CID197152), seychellene (CID519743), and alpha-patchouli alcohol isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) to COX-1 and COX-2. Molecular docking tools proposed by Hex 8.0 were employed in this research. Discovery Studio Client 3.5 software tool and virtual molecular dynamic 1.9.1 software were also used to visualize the molecular interactions identified in this research. In order to calculate the binding energy of the molecular dynamic interaction, AMBER12 software was utilized. Results of the analysis on all sesquiterpenoid indicate that those compounds were the inhibitors of COX-1 and COX-2. Overall, the binding energy calculations (using PBSA Model Solvent) of alpha-patchouli alcohol (CID521903) and seychellene (CID519743) have been identified as the candidates of non-selective inhibitor; alpha-bulnesene (CID94275), alpha-guaiene (CID107152), and alpha-patchouli alcohol isomers (CID6432585, CID3080622, CID10955174, CID56928117) have been suggested as the candidates for a selective COX-1 inhibitor; whereas alpha-patchouli alcohol (CID442384) was the candidate for a selective COX-2 inhibitor.

1. Introduction

The rapid development of high-performance computing intensifies the competition to invent faster supercomputers which invention is considered as a national pride. High-performance computing machines are highly valued for having adequate ability to solve complex problems related to national interests in several sectors including national defense, energy, financial sectors and science. Within the global economic growth competition and the advancement in science and technology including the advancement in biology, chemistry, pharmacy and medicine, supercomputers play key roles. In-silico analysis has been developed as the computational approach [1, 2].

Molecular interactions including protein-nucleic acid, drug-protein, protein-protein, enzyme-substrate, and drug-nucleic acid play important roles in many essential biological processes, such as enzyme inhibition, signal transduction, antibody-antigen recognition, transport, gene expression control, cell regulation, up to multi-domain proteins assembly. Stable protein-protein or protein-ligand complexes are often produced by the interaction which complexes are considered essential in performing their biological functions. To determine the binding mode and affinity between interacting molecules, tertiary structure of proteins should be first identified. Unfortunately, conducting experiments to obtain complex structures has been considered challenging and expensive because the experiments would require X-ray crystallography or NMR. Docking computation has been considered a feasible and important approach to comprehend the protein–protein or protein-ligand interactions [3]. Experimental technique has been frequently used to determine the three-dimensional protein structures—and structure of the databases such as Protein Data Bank (PDB) and Worldwide Protein Data Bank. A total of 88,000 protein structures have been identified and most of them are significant in critical metabolic pathways which might be regarded as potential therapeutic target. Therefore, specific databases that contain binary complexes structures would be available, along with information about binding affinities, such as in PDBBIND, PLD, AffinDB and BindDB, molecular docking procedures improvement [3, 4].

In silico virtual screening is a popular identification technique used in in drug discovery projects which distinguishes true actives from inactive or decoy molecules effectively. To have better comprehension on the dynamic behavior of protein drug targets, compound databases can be screened against an ensemble of protein conformations through experiments or generated-computation [5]. Screening states include ligand preparation, protein targets, molecular docking, visualization, binding energy calculation, and scoring [6]. A computer simulation procedure in the form of molecular docking is commonly used to predict the conformation of certain receptor-ligand complex, which receptor is usually a protein or a nucleic acid molecule or the ligands in the form of a small molecule or other protein, or sesquiterpenoid/sesquiterpenoid alcohol interaction to protein cyclooxygenase, as shown Figure 1. In modern structure-based drug design, accurate prediction is necessary to determine the binding modes between the ligand and protein. Virtual screening is the most popular docking application that selects molecules from an existing database for further research. As the demand on this computational method keeps increasing, people expected a fast and reliable method. Another application used in this study was molecular complexes investigation [3, 6, 7, 8, 9, 10, 11]. Previous studies have shown that dynamic molecular-generated conformations play considerable role in the identification of novel hit compounds because structural rearrangements obtained from molecular dynamic show novel-targetable areas. However, predicting the priori is still considered difficult, especially when a molecular dynamic conformation outperforms a virtual screening against the crystal structure.

Figure 1.

Molecular docking-Molecular dynamic ligand to protein.

This study evaluated whether molecular dynamic conformations lead to better virtual screening performance for nine ligand sesquiterpenoid/sesquiterpenoid Pogostemon herba to protein cyclooxygenase (COX-1/COX-2). The results of in silico analysis data will be completed with IC₅₀ value determination and in-vitro and in-vivo evaluation of the biological activity.

2. Theory of docking and virtual molecular dynamic

Within the process of living system, protein-ligand interactions have been known to play central roles. It has been considered interesting to obtain more comprehensive understanding of protein interactions with small molecules because it leads to better understanding of various functions and therapeutic intervention. As a matter of fact, molecular recognition is a complex interplay of several factors including inter-molecular interactions of protein, ligand and the surrounding solvent, conformational variations of binding partners and the thermodynamics of molecular association. The non-covalent reversible binding of small-molecules to proteins also plays a central role in the field of biology. Several processes are known crucial in living systems that involve specific recognition of small molecule ligands by proteins. For instance, certain enzymes affect their substrates and catalyze chemical reactions inside the cells, where transporters recognize specific molecules based on the movement across membrane barriers, receptors that are specifically bind to hormones or other chemical messengers for inter- and intracellular communication. Finally, antibodies uniquely can bind to other chemical agents to mount vital defense mechanisms against infections and diseases. In general, protein-ligand binding in an aqueous environment is described as follows.

$$\text{Protein}{\left(\mathrm{P}\right)}_{\left(\mathrm{aq}\right)}+\text{Ligand}{\left(\mathrm{L}\right)}_{\left(\mathrm{aq}\right)}➔\text{Protein}-\text{Ligand}{\left(\mathrm{PL}\right)}_{\left(\mathrm{aq}\right)}$$ E1

A change in the free energy (ΔG) is always followed by chemical reactions and change in two other important quantities; enthalpy (∆H)—the heat content and entropy (∆S) that showed disorder of temperature-independent degree. The relationship between these quantities is shown as follows:

$${\mathrm{\Delta G}}^{°}=\mathrm{\Delta H}–\mathrm{T\Delta S}$$ E2

Some factors including electrostatic and van der Waals interactions, ionization effects, conformational changes and the role of solvent affect the changes in the binding of free energy. Those factors are manifested as either favorable or unfavorable changes in entropy and enthalpy. In order to create a spontaneous reaction, the free energy change should be negative at equilibrium, ∆G° which relates to the equilibrium constant (K) in this following expression:

$${\mathrm{\Delta G}}^{°}=-\text{RTlnK}$$ E3

where R is the gas constant and T is the absolute temperature. Using this relationship, free energy changes can be derived from experimentally measurable quantity, K. Biological K values exhibit a wide range from weak to very strong binding.

Scoring function in ligand-protein docking is expected to identify the preferred binding poses of ligands. However, considering the computational efficiency, approximations were usually introduced into the scoring functions, which unfortunately, often impair the prediction accuracy [12]. The scoring functions of ligand–protein docking can be roughly categorized into two classes: force field-based scoring functions and knowledge-based scoring functions. A force field-based scoring consists of a few potential terms, such as van der Waals interactions, electrostatic interactions, hydrophobic effects, desolvation energies, and entropic effects, and the total energy of a conformation is calculated by summing up the contributions of all energy terms [13, 14].

Inter-atomic interactions mediate the non-covalent binding of small-molecule ligand to proteins. The interactions usually include electrostatic and van der Waals interactions (Figure 2). The affinity of receptor-ligand binding is strongly determined by other factors such as entropy, desolvation, flexibility of receptor structure and the structural water molecules in the binding site [15, 16]. A brief literature review of the importance of protein-ligand interactions and other factors contributing to binding affinity is described as follows.

Figure 2.

Major type of non-bonded interactions in protein-ligand complexes [17].

Protein-ligand electrostatic complementarity and the ligand at the binding interface are both vital for the formation of complex. The predominant types of electrostatic interactions appear in the form of hydrogen bonding, salt bridges, and metal interactions.

As the most important directional interaction in biological macromolecules, hydrogen bonding is known for conferring stability to protein structure and selectivity to protein-ligand interactions [17]. Hydrogen bonding normally occurs between two electronegative atoms, which donor is covalently bound to hydrogen atom, while the acceptor contains a lone pair of electrons. The attractive interaction between partial positive charge in the hydrogen atom and partial negative charge on the acceptor atom forms strong electrostatic attraction. Several theoretical and experimental studies have successfully confirmed an additional covalent component to hydrogen bonds based empty σ^* anti- bonding orbital of the hydrogen atom and highest occupied orbital of the acceptor interaction [17, 18]. In hydrophobic interactions, non-polar parts of the molecule interact (Figure 2). The non-polar parts of protein-ligand complexes at the interacting surfaces are covered by the binding causing water molecules displacement which eventually increases the entropy. Hydrophobic interactions are entropy-driven and they play crucial roles in ligand binding [16, 18].

Finding an accurate modeling of protein-ligand binding is an extremely challenging task due to its complexity. Usually, thermodynamics and statistical mechanical principles are employed to develop relatively accurate, but computationally demanding treatment of protein-ligand interactions. In this method, full-scale molecular dynamic simulation using explicit solvent and flexible protein and ligand molecules is employed [15, 17]. Both, absolute and relative binding energy can be measured using free energy approach. The absolute binding free energy method which has been considered an accurate method; it involves separate simulation treatment for solvated protein, ligand and the complex. Prior information is not quite necessary regarding the structure and binding affinity of the complex. A well-known structure for the complex is used within the context of relative free energy calculation as a reference, while the gaps in binding free energy are measured for the ligand of interest. Measurement can be carried out in the form of alchemical transformation of reference ligand into target ligand. Molecular dynamics is used in exhaustive sampling of the configuration space. The accuracy of these methods is determined by the underlying atomic force field and proper selection of protocol to address certain problem at hand [17, 18].

Scoring functions are based on knowledge focus on the optimization of specified terms and they use other optimization methods to set the best weighing for each scoring term regarding the training sets. Hence, these methods are often called as informatics-driven methods including IFACE [19, 20], DARS [21], SPA-PP [22], DrugScorePPI [23], TS [24], etc. In the past decades, free energy perturbation (FEP) [25] as well as energy representation (ER) [21] that are more theoretically rigorous free energy calculation methods, has emerged. Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) [26] and Molecular mechanics/generalized Born surface area (MM/GBSA) [27] are commonly employed to estimate ligand-protein binding free energies. Unfortunately, these methods are rather time-consuming compared to other scoring functions. Yet, rapid advancement of computer hardware technology seems promising, and it is expected to allow these methods to be used in protein-protein docking in the near future.

The MM/PBSA and MM/GBSA approaches are more computationally efficient when they were compared to thermodynamic integration (TI) and free energy perturbation (FEP) approaches. Besides, they allow decomposition into different interaction terms to occur [13, 26, 28]. MM/PBSA and MM/GBSA are more efficient in tem of computation. Another similar approach is the linear interaction energy (LIE) method, which calculates the average energy interaction in MD simulations to estimate the absolute binding free energy. Similarly, LIE restricts the simulations only to two end points of ligand binding. Different from most molecular docking empirical scoring functions, MM/PBSA and MM/GBSA do not demand a large training set to fit different parameters for each energy term [25]. In addition, MM/PBSA and MM/GBSA allow rigorous free energy decomposition into contributions to occur, originating from different groups of atoms or types of interaction [29, 30]. The binding free energy (∆G_bind) between a ligand and a receptor protein offered in these methods to form a complex Receptor protein - Ligand is calculated as follows.

$$\Delta {\text{G}}_{\text{bind}}={\text{G}}_{\text{complex}}\left({\text{G}}_{\text{protein}}+{\text{G}}_{\text{ligand}}\right)$$ E4

$$\Delta {\text{G}}_{\text{bind}}=\mathrm{\Delta H}-\mathrm{T\Delta S}˜\Delta {\text{E}}_{\text{MM}}+\Delta {\text{G}}_{\text{sol}}–\mathrm{T\Delta S}$$ E5

$$\Delta {\text{E}}_{\text{MM}}=\Delta {\text{E}}_{\text{internal}}+\Delta {\text{E}}_{\text{electrostatic}}+\Delta {\text{E}}_{\text{vdW}}$$ E6

$$\Delta {\text{G}}_{\text{sol}}=\Delta {\text{G}}_{\text{PB/GB}}+\Delta {\text{G}}_{\text{SA}}$$ E7

where ΔG_bind shows the free energy of ligand-protein total binding; ΔE_MM reflects the total gas phase energy (sum of ΔE_internal, ΔE_{electrostatic}, and ΔE_vdw); ΔG_sol is the sum of polar (ΔG_PB/GB) and non-polar (ΔG_SA) contributions to solvation; and -TΔS refers to the conformational binding entropy (commonly calculated by normal-mode analysis). ΔE_internal shows the internal energy that arises from different bond, angle, and dihedral in molecular mechanics (MM) force field (in the MM/PBSA and MM/GBSA, this always amounts to zero as shown in single trajectory of a complex calculation). ΔE_{electrostatic} and ΔE_vdw are the electrostatic and van der Waals energies resulted from the calculation of MM, while ΔG_PB/GB shows the polar contribution to the solvation free energy (calculated using Poisson–Boltzmann (PB) or generalized Born (GB) method). ΔG_SA is the nonpolar solvation free energy that is usually computed using a linear function of the solvent-accessible surface area (SASA). E_MM is molecular mechanical energy calculated from CHARMM force field, G_elec and G_NP are electrostatic polar components and non-solvation free energy. TS term refers to the entropy of the solute which is assumed to be constant between one set of poses for the same ligand on the active side. E_MM is a gas phase forcefield energy and consists of internal energy (E_int), electrostatic energy (E_elec) and van der Waals energy components. E_int is further divided into E_bond, E_angle, E_torsion and E_oop to calculate account energy related to bonds, angles, torque and outside as shown in Figure 3.

Figure 3.

Calculation of binding energy using the MM-GB/SA approach [36].

MM/GBSA and MM/PBSA have been successfully applied to predict the binding free energies for various ligand sesquiterpenoid/sesquiterpenoid alcohol to protein COX-1/COX-2, but the previous studies mostly focused on certain specific systems and the prediction results cannot afford the overall accuracy of MM/PBSA and MM/GBSA for ligand-protein systems.

3. Ligand sesquiterpenoid

Ligand sequiterpenoid was obtained from pubchem.ncbi.nlm.nih.gov, such as (alpha-bulnesene (CID94275), alpha-guaiene (CID107152), and seychellene (CID519743). Also, sesquiterpenoid alcohols, including alpha-Patchouli alcohol isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) as 3D-SDF format. Then, its energy was minimized which files were converted to 3D-PDB format by Open Babel 2.3.1 in Hex 8.0 as the ligands prepared for virtual screening [6, 8, 11, 31]. The structures of studied ligands are shown in Figure 4.

Figure 4.

Structure of sequiterpenoid and sesquiterpenoid alcohols Pogostemon herba.

Sesquiterpenoid, such as alpha-bulnesene (CID94275), alpha-guaiene (CID107152), and seychellene (CID519743) has molecular weight 204.35 g/mol, molecular formula C₁₅H₂₄, and XLogP3-AA: 4.60; 4.6; and 5.10 respectively (Table 1). And alpha-patchouli alcohol isomers has molecular weight: 222.36634 g/mol; molecular formula:

Description	sequiterpenoid			sesquiterpenoid alcohol
	CID94275	CID107152	CID519743	CID442384, CID521903, CID6432585, CID3080622,CID1095517,CID56928117
Molecular Weight (g mol^−1)	204.35	204.35	204.35	222.36
Kinase inhibitor	−1.33	−1.33	−1.30	−0.88
Nuclear receptor inhibitor	0.19	0.19	0.27	0.55
Protease inhibitor	−0.60	−0.60	−0.50	−0.32
Enzyme inhibitor	0.07	0.07	0.28	0.40
xlogP3-AA	4.6	4.60	5.10	4.10
H-Bond donor	0	0	0	1
H-Bond aseptor	0	0	0	1

Table 1.

Physical–chemical properties and predicted activity of sesquiterpenoids and sequiterpenoid alcohols Pogostemon herba.

C₁₅H₂₆O; XLogP3-AA: 4.1; H-Bond Donor: 1; and H-Bond Acceptor: 1. The number of isomers of alpha-Patchouli alcohol is six. In Figure 4 2D-sesquiterpernoid/sesquiterpenoid alcohol, such as alpha-bulnesene (CID94275), alpha-guaiene, and seychellene (CID519743), and alpha-Patchouli isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) show the different position of hydroxyl group and hydrogen atom. The 3D structure of sesquiterpenoid/sesquiterpenoid alcohol was retrieved in 3D-SDF format from http://pubchem.ncbi.nlm.nih.gov/. For the preparation of docking, 3D-SDF format of isomers was converted to 3D-PDB using open babel software. This program helps to search, convert, analyze, or store data which has a wide range of applications in the different fields of molecular modeling, computational chemistry, and so forth. For a common user, it helps to apply chemistry aspects without worrying about the low level details of chemical information. It also converts crystallographic file formats (CIF, ShelX), reaction formats (MDLRXN), molecular dynamics and docking (AutoDock, Amber), 3D viewers (Chem3D, Molden), and chemical kinetics and thermodynamics (ChemKin, Termo) [6, 8].

4. Cyclooxygenase protein receptor (COX-1 and COX-2)

3D model from PDB ID: 1PTH was obtained from SWISS-MODEL repository for cyclooxygenase-1 (COX-1) (http://www.rcsb.org/pdb/explore/explore.do?structureId=1pth) and 3D model from PDB ID: 6COX for cyclooxygenase-2 (COX-2) [6, 8]. We used Ramachandran plot analysis for validation protein receptor [32].

In Figure 5 shows the Ramachandran plot analysis of COX-1 and COX-2 protein receptor before rigid docking. It showed that COX-1 protein receptor had 97.5% favored regions, 2.4% allowed regions, and, 0.2% outlier regions. Whereas, COX-2 protein receptor had 81.9% favored regions, 15.4% allowed regions and 2.7% outlier regions. Ramachandran plot displays the main chain torsion angles phi, psi (φ, Ψ) (Ramachandran angles) in a protein of known structure. The model was verified to guarantee the validity of programming and algorithms implemented. Results of the validity test showed that amino acid residues were distributed at the most favorable region in the Ramachandran plot. This is an indication of the stereochemical quality of the model taken for the structural analysis and also validated the target-ligand binding efficacy of the structure. The Ramachandran plot presents the angle of phi-psi torsion of all residues in the structure (except those at the chain termini) which were classified according to their regions in the quadrangle. The most favored regions are colored yellow, additional allowed/generously allowed region, and outlier regions are indicated as blue and pink fields, respectively [6, 8, 32].

Figure 5.

Ramachandran plot analysis of COX-1 and COX-2.

5. Molecular docking ligand and binding energy interaction sesquiterpenoid/sesquiterpenoid alcohols to protein COX-1 and COX-2

We used Hex8.0 software (http://hex.loria.fr) for rigid docking to compute possible interaction COX-1 and COX-2 with (alpha-bulnesene (CID94275), alpha-guaiene (CID107152), seychellene (CID519743) and sesquiterpenoid alcohols such as alpha-Patchouli alcohol isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) on the interaction site. Output of the docking was refined using Discovery Studio Client 3.5 software. We used Discovery Studio Client 3.5 to perform interactions, ligand binds to COX-1/COX-2 and Ramachandran plot analysis.

The repeat rigid docking used Hex 8.0 software to compute possible interaction COX-1 and COX-2 with sesquiterpenoid/sesquiterpenoid alcohols such as alpha-bulnesene (CID94275), alpha-guaiene, seychellene (CID519743), and alpha-patchouli alcohol isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) on its interaction site and the data are represented by Discovery Studio 3.5 software in (Figure 5(a1–l1)). The interaction site position of COX-1/COX-2-sesquiterpenoid/sesquiterpenoid alcohol complexes were analyzed using Discovery Studio-3.5 Client software to get the receptor-ligand interaction and Ramachandran plot, as shown in Figure 5; some of them are alpha-patchouli alcohol isomers-COX-1/COX-2 complexes.

In Table 2 and Figure 6, the interactions active site of ligand sesquiterpenoid/sesquiterpenoid alcohol with COX-1 and COX-2 protein receptor showed the differences in the position active site. The different positions were analyzed and presented in the Ramachandran plot analysis and its amino acid residues in the receptor active site of COX-1 and COX-2 in which hydrogen atoms and hydroxyl groups on each of the 3D-isomers of alpha-patchouli alcohol structure were in different position (Figure 4). The results of docking and analysis of the active site also show that all ligands sesquiterpenoid/sesquiterpenoid alcohol are in the catalytic domain. Thus, all the compounds have the capability of blocking oxygenated reaction and reaction peroxides; currently substrate arachidonic acid becomes PGH2.

No.	Virtual modeling	Amino acid residues in the active site (by Hex 8.0 software and then; Discovery Studio 3.5 software)
		COX-1	COX-2
1.	alpha-Patchouli alcohol CID442384	TRP141A, GLU142A, SER145A, ASN146A, LEU226B, GLY227B, ASP231B, GLN243B, GLY237B, ASN239B, LEU240B, ASP238B, ARG335B	TRP139A, GLU140A, SER143A, ASN144A, LEU145A, GLY235B, GLU236B, THR237B, LEU238B, GLN241B, GLN330B
2.	alpha-Patchouli alcohol CID521903	SER123A, ASN124A, LEU125A, ILE126A, PRO127A, SER128A, PRO129A, GLN372A, PHE373A, GLN274A, LYS534A, PRO544B, GLU545B	LYS211B, THR212B, ASP213B, HIS214B, LYS215B, ARG222B, ILE274B, GLN298B, GLU290B, VAL291B, HEM682B
3.	alpha-Patchouli alcohol CID643285	TRP141A, GLU142A, SER145A, ASN146A, LEU226B, ASP231B, GLY237B, ASP238B, ASN239B, LEU240B, GLN243B, ARG335B	TRP139A, GLU140A, SER143A, ASN144A, LEU145A, THR237B, LEU238B, GLY235B, GLU236B, GLN241B, GLN330B, LYS333B
4.	alpha-Patchouli alcohol CID3080622	SER123A, ASN124A, ILE126A, PRO127A, SER128A, PRO129A, GLN372A, PHE373A, GLN374A, LYS534A, PRO544B, GLU545B	TRP139A, SER143A, ASN144A, LEU145A, GLY235B, GLU236B, THR237B, LEU238B, GLN241B, LYS333B
5.	alpha-Patchouli alcohol CID10955174	TRP141A, GLU142A, SER145A, ASN146A, LEU226B, ASP231B, GLY237B, ASP238B, ASN239B, LEU240B, GLN243B, ARG335B	TRP139A, GLU140A, SER143A, ASN144A, LEU145A, GLY235B, GLU236B, THR237B, LEU238B, GLN241B, GLN330B, LYS33B
6.	alpha-Patchouli alcohol CID56928117	Electrostatic: ASN146B.; Van der Walls: LEU226A, GLY237A, ASP238A, ASN239A, LEU240A, GLU241A, GLN243A, ARG335A, TRP141B, GLU142B, SER145B	Electrostatic: SER143B.; Van der Walls: GLY235A, GLU236A, THR237A, LEU238A, ASP239A, GLN241A, LYS333A, TRP139B, GLU140B, ASN144B, LEU145B
7.	alpha-bulnesene CID94275	Van der Walls: VAL147A, LYS224A, ALA225A, LEU226A, GLY227A, ASP231A, GLY233A, GLY237A, ASP238A, ASN239A, LEU240A, ARG335A, TRP141B, GLU142B, SER145B, ASN146B, VAL147B	Van der Walls: GLY225A, ASP229A, GLY235A, GLU236A, LEU238A, GLN241A, GLN330A, THR237A, LYS333A, SER143B, TRP139B, GLU140B, ASN144B, LEU145B
8.	alpha-guaiene CID107152	Van der Walls: TRP141A, GLU142A, SER145A, ASN146A, LEU226B, GLY227B, ASP231B, GLY237B, ASN239B, ASP238B, LEU240B, GLU241B, GLN243B, ARG335B	Van der Walls: GLY225A, ASP229A, ASN231A, GLY235A, GLU236A, THR237A, GLN241A, GLN330A, LYS333A, TRP139B, GLU140B, SER143B, ASN144B, LEU145B, LEU238A
9.	Seychellene CID519743	Van der Walls: PRO544A, GLU545A, SER123B, ASN124B, LEU125B, ILE126B, PRO127B, SER128B, PHE373B, GLN372B, GLN374B, LYS534B	Van der Walls: ASP213A, HIS214A, LYS215A,LYS211A, THR212A, ARG222A, ILE274A, GLN289A, GLU290A, VAL291A, HEM682A

Table 2.

Analysis of virtual modeling of COX-1/COX-2-sesquiterpenoid/sesquiterpenoid alcohol complexes.

Figure 6.

Modeling analysis alpha-patchouli alcohol isomer in complex with COX-1 and COX-2. (a1) – (l1) 3D active site structure of COX-1/COX-2-alpha-patchouli alcohol isomers complexes; (a2) to (l2) Ramachandran plot analysis of COX-1/COX-2-alpha-patchouli alcohol complexes using discovery studio 3.5 viewer Software.

Each ligand, CID521903, was seen interacting with HEM682B group in COX-2-CID521903 complexes. This result proved that it would lead to inhibition of enzymatic reactions occurring COX-1 and COX-2. The analysis of active site showed that there are any difference and similarities of the active site of all ligand alpha-patchouli alcohol isomers which is interact with receptor proteins COX-1 and COX-2. This difference is caused by different stereoisomers of hydrogen atoms and hydroxyl group in alpha-patchouli alcohol isomers. The different position active site the complexes have led to interaction types, such as hydrogen bond, van der Waals, electrostatic and covalent bond. The different types of interactions in this complex will certainly affect its binding free energy.

6. Molecular dynamic screening of sesquiterpenoid/sesquiterpenoid alcohol Pogostemon herba as predicted cyclooxygenase inhibitor selective

After the results of the rigid docking to compute possible interaction COX-1 and COX-2 with (alpha-bulnesene (CID94275), alpha-guaiene (CID107152), and seychellene (CID519743). And also, sesquiterpenoid alcohol, such as alpha-Patchouli alcohol isomers (CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117) to performed active visualization-interaction 2D and 3D, and binding energy using Discovery Studio 3.5 software. The output of the docking, visualization, and binding energy calculation using AMBER12 software and Virtual Molecular Dynamics 1.9.1 obtained the most possible native complex structure of sesquiterpenoid/sesquiterpenoid alcohol of CID94275, CID107152, CID519743, CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117, respectively, that bind with COX-1 and COX-2 in molecular dynamic with Model Solvent of MM-PBSA (Molecular Mechanics Poisson-Boltzmann Surface Area), which included both backbone and side-chains movements. Therefore, we used AMBER12 to refine the candidate models according to a binding energy calculation for scoring of virtual screening sesquiterpenoid/sesquiterpenoid alcohol compounds as selective inhibitor for COX-1 and/or COX-2. Molecular dynamics (MD) were carried out using AMBER12 and the AMBER-99 force field. The initial structure of the sesquiterpenoid/sesquiterpenoid alcohol inhibitor complex was taken for each compound from the Hex 8.0 docking study. The ligand force fields parameters were taken from the General Amber force Field (GAFF), whereas AM1 ESP atomic partial charges were assigned to the inhibitors. Prior to the free MD simulations, two steps of relaxation were carried out; in the first step, we kept the protein fixed with a constraint of 500 Kcal·mol⁻¹ · °A⁻¹. In the second step, the inhibitor structures were relaxed for 0.5 pico second, during which the protein atoms were restrained to the X-ray coordinates with a force constant of 500 Kcal·mol⁻¹ · °A⁻¹. In the final step, all restraints were removed and the complexes were relaxed for 1 pico second. The temperature of the relaxed system was then equilibrated at 300 Kelvin through 20 pico second of MD using 2 fs time steps. A constant volume periodic boundary was set to equilibrate the temperature of the system by the Langevin dynamics using a collision frequency of 10 ps⁻¹ and a velocity limit of five temperature units. During the temperature equilibration routine, the complex in the solvent box was restrained to the initial coordinates with a weak force constant of 10 Kcal·mol⁻¹ · °A⁻¹. The final coordinates of the temperature equilibration routine (after 20 ps) were then used to complete a 1 ns molecular dynamics routine using 2 fs time steps, during which the temperature was kept at 300 Kelvin. For the Langevin dynamics a collision frequency of 1 ps⁻¹ and a velocity limit of 20 temperature units were used. The pressure of the solvated system was equilibrated at 1 bar at a certain density in a constant pressure periodic boundary by an isotropic pressure scaling method employing a pressure relaxation time of 2 ps. The time step of the free MD simulations was 2 fs with a cut-off of 9°A for the non-bonded interaction, and SHAKE was employed to keep all bonds involving hydrogen atoms rigid. Calculation of binding energy was administered using this equation:

$$\Delta {\text{G}}_{\text{bind}}={\text{G}}_{\text{complex}}-\left[{\text{G}}_{\text{protein}}+{\text{G}}_{\text{ligand}}\right]$$

[6, 8, 33, 34, 35, 36, 37].

We were using AMBER12 software and Virtual Molecular Dynamics 1.9.1 to simulate the most possible native complex structure of sesquiterpenoid/sesquiterpenoid alcohol (CID94275, CID107152, CID519743, CID442384, CID521903, CID6432585, CID3080622, CID10955174, and CID56928117), respectively, that binds with COX-1 and COX-2 in molecular dynamic with MM-PBSA Model Solvent. The MD simulations of the sesquiterpenoid/sesquiterpenoid alcohol-inhibitor, some of them are alpha-patchouli alcohol-COX-1/COX-2 complexes. The structure of the complexes is shown in Figure 7(a–f) and (j–o). We also acquire the results of the analysis of 200 poses: the complex energy, energy ligand protein, and energy. The binding energy was calculated use the following equation:

Figure 7.

Binding energy calculation of alpha-patchouli alcohol isomers binds to COX-1/COX-2. (a) and (f) and (j)–(o) virtual molecule dynamic complexes of COX-1/COX-2-alpha-patchouli alcohol isomers. (g), (h), (i), (p),(q), and (r) comparison of binding energy calculation of alpha-patchouli alcohol isomer-COX-1 (blue) and COX-2 (red) complexes. (s) Histogram of binding energy calculation of COX-1 (blue)/COX-2 (red) sesquiterpenoid/sesquiterpenoid alcohol complexes by discovery studio 3.5 (s-1) and Amber 12 (s-2).

$$\Delta {\text{G}}_{\text{bind}}={\text{G}}_{\text{complex}}-\left[{\text{G}}_{\text{protein}}+{\text{G}}_{\text{ligand}}\right]$$

as shown in Figure 7(g–i), (p–r) and (o).

Analysis of the active site and the binding energies COX-1/COX-2-sesquiterpenoid/sesquiterpenoid alcohol are by Discovery Studio 3.5 and Amber 12, summarized and presented in Figure 7 (s).

The different position active site the complexes have led to interaction types, such as hydrogen bond, van der Waals, electrostatic and covalent bond. The different types of interactions in this complex will certainly affect its binding free energy. The use of Poisson-Boltzmann (PB) and Generalized Born (GB) characterized the binding free energy calculation model solvent MMPB/SA and MM-GB/SA in computing the electrostatic component of the solvation free energy. The following equation was employed in binding free energy of the protein-ligand complex.

$$\mathrm{\Delta G}=\mathrm{\Delta H}–\mathrm{T\Delta S}$$ E8

T is the temperature of the system at 300 Kelvin. The free binding energy (ΔGbinds) of the protein-ligand-complex were evaluated using MMPBSA (Molecular Mechanics Poison Blotzmann Surface Area) method as implemented in Discovery Studio 3.5 and AMBER12. MMPBSA has always been considered as a proper method to compare binding energies of similar ligands. MMPBSA measures the binding free energy based on thermodynamic cycle in which molecular mechanical energy and the continuum solvent approaches are simultaneously used [6, 8, 33, 38]. The calculation of binding free energy is computed as:

$$\Delta {\text{G}}_{\text{bind}}={\text{G}}_{\text{complex}}-\left[{\text{G}}_{\text{protein}}+{\text{G}}_{\text{ligand}}\right]$$ E9

In (5.2), G_complex is the absolute free energy of the complex, G_protein is the absolute free energy of the protein, and G_ligand is the absolute free energy of the ligand [6, 8, 33, 38]. The free energy of each term was estimated as a sum of the three terms:

$$\left[\mathrm{G}\right]=\left[{\text{E}}_{\text{MM}}\right]+\left[{\text{G}}_{\text{sol}}\right]-\mathrm{T}·\left[\mathrm{S}\right]$$ E10

[G_MM] is the molecular mechanics energy of the molecule expressed as the sum of the internal energy (bond, angle, and dihedral) (E_int), electrostatic energy (E_ele), and van der Waals term (E_vdw):

$$\left[{\text{E}}_{\text{MM}}\right]=\left[{\text{E}}_{\text{int}}\right]+\left[{\text{E}}_{\text{ele}}\right]+\left[{\text{E}}_{\text{vdw}}\right]$$ E11

[E_ele] solvation energy can be categorized as polar and nonpolar part. Polar part gives electrostatic contribution to solvation by solving the linear Poisson Boltzmann equation within the solvent’s continuum model [33]. The binding energy calculation in AMBER12 includes preparation, minimization, heating, and energy calculations (complex, protein, and ligand). We extracted 200 snapshots (at time intervals of 2 ps) for each species (complex, protein, and ligand). Furthermore, the visualization using virtual model dynamic (VMD 1.9.1 software) is shown in Figure 7(a–f) and (j–o), and then the binding energy calculation can be obtained from the data ligand energy, protein energy, and energy complex by AMBER12, 200 times/poses, respectively; next, the binding free energy calculation is calculated by Eq. (7.2) and shown in Figure 7(g), (h), (i), (p), (q), and (r) and summarized in Figure 5(s). Figure 5(s) shows that the binding energy calculation (PBSA Model Solvent) of COX-1 CID442384 complexes (−28.386 ± 1.102 Kcal/mol) was smaller than the COX-2 CID442384 complexes (−16.215 ± 0.985 Kcal/mol) and also ligands CID6432585, CID3080622, CID10955174, and CID56928117. The similar research, docking studies ligand salicin compound from D. gangeticum to COX-1 and COX-2 protein receptor, showed high binding affinity COX-2 protein (−5 Kcal/mol) and lesser interaction with COX-1 (−3.79 Kcal/mol). Therefore, salicin could predict as COX-2 inhibitor selective and anti-cancerous compound [6].

E_binds (ΔG) was determined on the basis of calculation of the Eq. (5.2). Gligand value is influenced by the type of ligand. Gligand will affect the value Ebinds and ratio of E_binds COX-1 and E_binds COX-2. Hence, in-silico analysis can be used as an approach to determine the selectivity of the ligand as an inhibitor of COX-1/COX-2. Ebinds (binding energy calculations) seychellene (CID519743) (Figure 7(s)) showed as candidate non-selective COX inhibitor and it’s similar to value of selective IC₅₀, as shown in Figure 8.

Figure 8.

Regression linier analyses of IC₅₀ fraction-5 to COX-1 and COX-2 [39].

The relationship binding energy, K_i and IC₅₀ is defined by Eq. (5.5) and (5.6) [39, 40].

$$\Delta {\text{G}}_{\text{bind}}=2.303\mathrm{R}·\mathrm{T}\text{log}{\mathrm{K}}_{\mathrm{i}}$$ E12

$$\text{For competitive inhibition}:{\mathrm{K}}_{\mathrm{i}}=\left({\mathrm{IC}}_{50}–\mathrm{E}/2\right)/\left(\mathrm{S}/\mathrm{Km}+1\right)$$

$$\text{For uncompetitive inhibition}:{\mathrm{K}}_{\mathrm{i}}=\left({\mathrm{IC}}_{50}–\mathrm{E}/2\right)/\left(\mathrm{Km}/\mathrm{S}+1\right)$$

$$\text{if}\mathrm{S}=\mathrm{Km},{\mathrm{K}}_{\mathrm{i}}={\mathrm{IC}}_{50}/2\text{;}$$

$$\text{if}\mathrm{S}>>{\mathrm{K}}_{\mathrm{m}},{\mathrm{K}}_{\mathrm{i}}<<{\mathrm{IC}}_{50}\text{;}$$

$$\text{if}\mathrm{S}<<{\mathrm{K}}_{\mathrm{m}}{\mathrm{K}}_{\mathrm{i}} \approx {\mathrm{IC}}_{50}.$$

For non-competitive inhibition: K_i = IC₅₀ when S = K_m or S << K_m and for tightly bound inhibitor:

$${\mathrm{K}}_{\mathrm{i}}={\mathrm{IC}}_{50}–\mathrm{E}/2$$ E13

where: E = enzyme, S = Substrate, P=Product.

The latest development is more selective selective COX-2 drugs, such as valdecoxib (Bextra™) and etoricoxib (Arcoxi™) and lumiracoxib (Prexige). Several COX-2-selective drugs in NSAIDs are presented in Figure 9. The classification of COX inhibition is based on the potential inhibition of COX isoforms and specifically the IC₅₀ ratio of COX-1 and COX-2 (or selectivity index) [20].

Figure 9.

The relative selectivity of COX-1 and COX-2 inhibitors based on the IC₈₀ ratio is declared logarithmic, so 0 is the baseline, that is, the compound in the line is equiactive to COX-1 and COX-2. Compounds above the COX-1-selective line and below are COX-2 selective [34].

Eq. (5.5) can be used as the COX-1/COX-2 selectivity approach in in-silico analysis, which without calculating for competitive, un-competitive and non-competitive, shows that ΔG_bind are directly proportional to IC50 values.

While the selectivity of COX-1/COX-2 is expressed in the equation:

$${\mathrm{IC}}_{50}\text{selectivity},\mathrm{COX}-1/\mathrm{COX}-2=\text{log}\left({\mathrm{IC}}_{50},\text{ratio}\left(\mathrm{COX}-2/\mathrm{COX}-1\right)\right)$$ E14

Therefore selectivity in in-silico analysis can be expressed as:

$${\mathrm{E}}_{\mathrm{bind}} \text{selectivity},\mathrm{COX}-1/\mathrm{COX}-2=\text{log}\left({\mathrm{E}}_{\mathrm{bind}},\text{ratio}\left(\mathrm{COX}-2/\mathrm{COX}-1\right)\right)$$ E15

According Eq. (5.5), selectivities Ebinds and selectivities IC50 some of them are complexes of CID442384, CID519743, CID3060622, CID107152, and CID94275 with COX-1/COX-2, as shown in Figure 10.

Figure 10.

Selectivities of IC₅₀ versus E_binds sesquiterpenoid/sesquiterpenoid alcohol Pogostemon herba to COX-1/COX-2.

Collectively, our results suggest that alpha-Patchouli alcohol (CID442384) as candidate COX-2 inhibitor selective, alpha-guaiene (CID107152), alpha bulnesene (CID94275), alpha patchouli alcohol isomers (CID3060622, CID6432585, CID10955174, and CID56928117) as candidate COX-1 inhibitor selective, and alpha-patchouli alcohol CID521903, seychellene as candidate COX non selective. These in silico analysis data will be completed with the biological activity analysis.

7. Conclusion

Exploration of sesquiterpenoid/sesquiterpenoid alcohol compounds as inhibitors of COX isoenzymes as development of group NSAIDs, was carried out by means of in silico tools. The binding energy calculation (using PBSA Model Solvent) of sesquiterpenoid/sesquiterpenoid alcohol compounds: alpha patchouli alcohol (CID521903) and seychellene (CID519743) were identified as the candidates of non-selective inhibitor; alpha bulnesene (CID94275), alpha guaiene (CID107152), and alpha-patchouli alcohol isomers (CID6432585, CID3080622, CID10955174, while CID56928117) had been suggested as the candidate for a selective COX-1 inhibitor. Whereas, alpha-patchouli alcohol (CID442384) was the candidate for a selective COX-2 inhibitor.

Acknowledgments

This study was supported by Doctoral program scholarship of “Sandwich-Like Program 2013,” DGHE, Ministry of Education and Culture, RI. The authors acknowledge all facilities of Bioinformatics Laboratory Department of Creation Core, Ritsumeikan University, for providing the in silico analyses. Special thanks are due to Prof. Takhesi Kikuchi (my sensei), Masanari Matsuoka, Antonius Christianto, Michirou Kabata, Sayaka Ohara, Yousuke Kawai, and working group Bioinformatic Laboratory, Biwako Kutsasu Campus, Ritsumeikan University, for helpful support analysis and discussions. And also my wife (Sri Wahjoeni), my son (Damara Setyowijayanto Raharjo), Prof. Fatchiyah, Prof. Chanif Mahdi, Prof. Sutiman, Dr. Nurdiana, Lectures of Biology Dept. Brawijaya University, Lectures Academic of Pharmacy and Food Analysis “Putra Indonesia Malang”, Yayasan Putera Indonesia Malang, and Pharmacy Dept. Faculty Medicine Brawijaya University.

Conflict of interests

The authors declare no conflict of interests.