Cyclocurcumin, a curcumin derivative, exhibits immune-modulating ability and is a potential compound for the treatment of rheumatoid arthritis as predicted by the MM-PBSA method.

The control and treatment of rheumatoid arthritis is a challenge in today's world. Therefore, the pursuit of natural disease-modifying antirheumatic drugs (DMRDs) remains a top priority in rheumatology. The present study focused on curcumin and its derivatives in the search for new DMRDs. We focused on prominent p38 mitogen-activated protein (MAP) kinase p38α which is a prime regulator of tumor necrosis factor-α (TNF-α), a key mediator of rheumatoid arthritis. In the present study, we used the X-ray crystallographic structure of p38α for molecular docking simulations and molecular dynamic simulations to study the binding modes of curcumin and its derivatives with the active site of p38α. The ATP-binding domain was used for evaluating curcumin and its derivatives. Molecular docking simulation results were used to select 4 out of 8 compounds. These 4 compounds were simulated using GROMACS molecular simulation platform; the results generated were subjected to molecular mechanics-Poisson Boltzmann surface area (MM-PBSA) calculations. The results showed cyclocurcumin as a potential natural compound for development of a potent DMRD. These data were further supported by inhibition of TNF-α release from lipopolysaccharide (LPS)-stimulated human macrophages following cyclocurcumin treatment.

Introduction

Rheumatoid arthritis (RA), an autoimmune disorder, evolves over time causing progressive damage of synovial-lined joints and variable extra-articular manifestations (1). RA is the most common autoimmune disease, with an incidence higher than psoriasis, Crohn's disease, type I (insulin-dependent) diabetes, lupus and multiple sclerosis (2). The progression of RA is characterized by symmetric polyarticular inflammation of the synovium, resulting in pain and stiffness and over a period of time leads to joint damage resulting in loss of function (3). The etiology of the autoimmune disorder leading to RA remains unknown, yet the mechanisms of how it evolves have been decoded to a great extent (4–7). RA affects 1% of the world population; in 2010 the global prevalence of this disease was 0.24% and was found to be more common in females (8–10). In Australian, Western European and North American regions the modeled age-standardized prevalence was found to be highest and the lowest was found to be in Asia and North African regions (10). In the USA, RA affects 1.3 million people, making it the third common type of arthritis among Americans (11).

The pathophysiology of this disease has been linked to genetic and environmental factors. These factors lead to a clinical syndrome spanning several disease subsets that require inflammatory cascades (12–14). The overexpression of tumor necrosis factor-α (TNF-α) is a prime inflammatory cascade involved in RA. This cascade leads to the overexpression of cytokines such as interleukin-1 (IL-1), IL-6 and IL-10, which drives persistent inflammation and joint destruction (15–18). These components of the chronic inflammatory response have been used for designing RA management therapies. To date, the management of RA has been achieved using disease-modifying antirheumatic drugs (DMARDs). These drugs are either of biological origin or are synthetic chemical compounds and based on this classification are termed as biological (bDMARDs) and synthetic (sDMARDs), respectively. Both bDMARDs and sDMARDs are characterized by their capacity to reduce or reverse signs and symptoms of RA. They interfere with the pathophysiology of this disease, which leads to improvements in disability, quality of life, ability to work, and suppression of joint damage (19–21).

DMARDs include etanercept, infliximab, methotrexate (MTX), sulfasalazine, leflunomide, rituximab, abatacept and others. However, these DMARDs are associated with numerous side effects (20). In order to overcome this factor more and more research is being conducted on natural comp ounds (21). Curcumin (diferuloylmethane) extracted from Curcuma longa is one of the many natural products that have been reported to be potential bDMARDs (22–25). Here in this study, we employed in silico techniques to investigate the natural compound curcumin and its derivatives as possible bDMARDs. Molecular docking simulations of curcumin and its derivatives with p38α, one of the prominent p38 mitogen-activated protein (MAP) kinase that regulates the expression of TNF-α and thus plays a role in RA was performed (26). p38α is involved in the p38 MAP kinase pathway that regulates a number of cellular processes critical to the etiology of RA. The activation and infiltration of leukocytes as well as the production of inflammatory cytokines are p38-dependent processes. Targeting of p38α has demonstrated efficacy in animal models of RA and various p38α inhibitors are currently in phase II and I clinical trials for RA. Thus using p38α as a therapeutic target for the development of a novel compound is a viable option. The 4 top complexes generated by targeting p38α with curcumin and its derivatives, were subjected to molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculations. The top compound, cyclocurcumin, was further studied, and inhibition of the release of TNF-α from LPS-stimulated macrophages by cyclocurcumin treatment confirmed its role as a potent p38α inhibitor.

Materials and methods

Data set

The X-ray crystallographic structure of p38α MAP kinase was downloaded from the RCSB Protein Data Bank (http://www.rcsb.org/pdb/home/home.do). The structure was cocrystallized with the pyrimidine-urea analogue 15 (PDB ID 2GHL, 2.10 Å X-ray resolution). The structure was stripped of the pyrimidine-urea analogue, where after removing water molecules, the charge and polar hydrogen atoms were added. The incomplete residues were corrected and the energy minimization of the modified structure was set to the default cut-off root mean square deviation (RMSD) value of 0.3 Å using OPLS 2001 force field by SPDB viewer (27). The structure of curcumin and its derivatives were retrieved from the PubChem database (Fig. 1).

Figure 1

Two dimensional structures of curcumin and its derivatives used in the present study.

Molecular docking simulations

Molecular docking simulation was carried out using AutoDock 4.2.6 tool (28). This tool uses binding free energy evaluation for mining the top compound-protein pose. The tool uses binding site definition to create a rectangular grid. The execution of the docking runs are performed within the grid. The energy values calculated by this suite are the culmination of the intermolecular energy of the protein, the internal energy of the ligand and the torsional free energy. The default parameters were set, and the maximum evaluations were performed to elucidate the top 10 binding modes for each run. The best compound among curcumin and its derivatives was selected for further evaluation at the atomic level using molecular dynamics simulation. The results were analyzed using Discovery Studio Visualizer (29).

Molecular dynamics simulation

GROningen MAchine for Chemical Simulations GROMACS 4.5.6 atomic simulation tool was used to evaluate the ligand-p38α complexes over 100 nanoseconds (nsec) (30). GROMOS 43a1 force field was used for p38α and the force field for the top ligand was calculated using PRODRG webserver (31,32). The complexes were kept in a water bath, where water molecules were represented by simple point charge (SPC216) model. The system was ensured to be overall charge neutral before being subjected to energy minimization process, where position restrain procedure was performed with NVT and NPT ensembles, both of which follow Berendsen coupling scheme. Time duration of 1 nsec with coupling constant of 0.1 psec at a constant temperature of 300 K was used for NVT ensemble. This was followed by NPT ensemble with a time duration of 1 nsec; here a constant pressure of 1 bar was subjected to a coupling constant of 5 psec. The complexes were then subjected to simulation for 100 nsec at 300 K. For each run, the SHAKE algorithm for bond lengths and particle mesh Ewald method for long-range electrostatic interactions was employed, where a 14 Å cut-off for van der Waals interactions, a 12 Å cut-off for Coulomb interaction with updates every 2 steps (30,33). The trajectories generated were subjected to GROMACS inbuild tools for evaluation. The results generated were visualized using Pymol, Ligplus and VMD (34-36) and the graphs were plotted by using Grace program (37).

MM-PBSA calculations

The trajectories generated were used for MM-PBSA calculations (38,39). For each trajectory 1001 complexes were retrieved, each after 100 psec. Binding modes were analyzed using Grace program (37). The tool gave us the binding free energy of the four complexes. The binding free energy (ΔGbind) calculated by the tool includes the following species:

In the equations Gcomplex, Gprotein and Gligand represent the free energy of the respective species, ΔEMM, ΔGsol and TΔS represent gas phase energy, solvation energy and an entropy term. The gas phase energy is the summation of the internal energy of the bonds, angle and torsion (ΔEval), electrostatic interaction energy (ΔEele), and van der Waals interaction energy (ΔEvdw). The salvation free energy comprises the polar salvation free energy (ΔGp), and the nonpolar salvation free energy (ΔGnp).

Principal component analysis

g_covar, g_anaeig in-build tools of GROMACS suite were used for indentifying the dominant motions of the selected complexes. PCA technique was used to trace these high-amplitude concerted motions in the p38α-curcumin-based complexes (40). Eigenvectors were generated from the covariance matrix generated by simple linear transformation. The first two eigenvectors (EV1 and EV2) having cosine content <0.2 were used to study free energy landscape (FEL). g_sham in-build tool was used to calculate FEL.

Cell culture

THP-1 monocytes were procured from American Type Culture Collection (ATCC; Manassas, VA, USA) and cultured in RPMI-medium (Sigma-Aldrich, St. Louis, MO, USA) supplemented with 10% FBS and 100 U/ml PenStrep (Gibco, Grand Island, NY, USA). The cells were grown under standard culture conditions at 37°C under 5% CO2 in a humidified incubator. For macrophage differentiation, monocytes were cultured with macrophage CSF (M-CSF) at 100 ng/ml. Stock solution of 100 mM cyclocurcumin (Sigma-Aldrich) was prepared in RPMI-medium and diluted to different concentrations. Cells were stimulated with 10 ng/ml lipopoly-saccharide (LPS) (Sigma-Aldrich) for TNF-α production.

Cell proliferation assay

3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay is a standard test to measure cell viability and is based on the conversion of MTT to formazan crystals by active mitochondrial dehydrogenases. Briefly, a cell suspension containing ~3×104 cells was plated into each well of a 96-well plate and allowed to attach for 24 h. Cyclocurcumin prepared in DMEM was added to the wells at a concentration of 10, 20, 40, and 60, 80 and 100 µM for 24 h. After treatment, 100 µl of MTT (Sigma-Aldrich) solution was added to each well at a concentration of 0.1 mg/ml (dissolved in PBS). After 4 h of incubation at 37°C in dark, the solution was removed and 100 µl of dimethyl sulfoxide (DMSO) was added to solubilize the produced formazan. Absorbance was measured at 570 nm using an automated microplate reader (BioTek Instruments, Winooski, VT, USA). The results are expressed as percentages relative to the controls. The percentage of proliferation inhibition was calculated as = (1 - ODsample/ODcontrol) ×100%.

TNF-α, enzyme-linked immunosorbent assay (ELISA)

Following treatment with cyclocurcumin and stimulation of human macrophages with LPS, the supernatants were harvested. Concentrations of TNF-α (BD Pharmingen, San Diego, CA, USA) were determined by ELISA, following the manufacturer's instructions. Absorbance was read at 450 nm on an ELISA plate reader (BioTek Instruments) using an in-built software program.

Statistical analysis

Statistical analysis was carried out by one-way ANOVA and the level of significance was tested at a P-value ranging from 0.01 to 0.001. Data are expressed as the mean ± SD of three independent experiments.

Results and Discussion

Molecular docking analysis

Molecular docking analysis was carried out using AutoDock 4.2 tool; 100 conformations for each compound were generated and the binding energy (kcal/mol) was set as the criteria for the cut-off for the molecular dynamics simulation analysis. Table I shows the results generated by docking curcumin and its selected derivatives against the active site of p38α. The top 4 compounds with the least binding energy were curcumin, bisdemethoxycurcumin, cyclocurcumin and demethoxycurcumin. Protein/ligand complexes from the molecular docking were subjected to post docking analysis, where weak interactions involved in the complex were studied (Fig. 2). Four complexes shortlisted for this analysis showed the favorable binding energy with a range -5 to -7 kcal/mol (Table II). The effective interactions studied here were supported by energies including final intermolecular energy, electrostatic energy and van der Waals energy. The effective interactions of the curcumin-p38α complex were evaluated using Discovery Studio Visualizer 4.5. The binding pocket of the complex involved GLY110, ASN115, ASP112, ASP168, VAL38, LEU104, GLU71, LYS53, VAL105, LEU75, LEU86, VAL30, THR106 and ALA111. The binding energy of curcumin-p38α was calculated to be the highest among the studied complexes, with no non-covalent interaction observed (Fig. 3A). The binding energy of the bisdemethoxycurcumin-p38α complex at -6.38 kcal/mol was the second highest calculated among the studied complexes. The bisdemethoxycurcumin-p38α complex formed one hydrogen bond (Fig. 4A). The interacting residue PHE169 of p38α is involved in active p38α inhibition. The oxygen atom of PHE169 formed a physical non-covalent interaction with bisdemethoxycurcumin, and the distance of the hydrogen bond formed was 2.36 A°. The binding pocket (Fig. 3B) of the bisdemethoxycurcumin-p38α complex involved LEU75, LEU74, THR106, ASP168, GLY170, GLU71, ILE84, PHE169, LEU167, MET109, GLY110, ASP112, ALA111 and ALA157. The third complex was that of cyclocurcumin with p38α protein kinase. The binding energy of the complex was calculated to be −6.12 kcal/mol. In the cyclocurcumin-p38α complex the results showed two hydrogen bonds (Fig. 3C) between the protein and ligand atoms. The interaction residues included LEU104 and VAL105 involved in active p38α inhibition. The binding pose of hit CID 69879809 clearly suggested effective inhibition. The atom 'O2' of the key residue LEU104 interacted with cyclocurcumin through hydrogen bond interactions (Fig. 4B). The key residue VAL105 of the active binding site cavity also interacted with cyclocurcumin oxygen atoms. Overall, the number of hydrogen bonds confirmed the high affinity of cyclocurcumin towards inhibition of p38α. In the case of the demethoxycurcumin-p38α complex, the results showed 3 hydrogen bonds between p38α and demethoxycurcumin (Fig. 3D). The interaction residues included SER154 and ASN155 (Fig. 4C). The interaction of an 'O' atom of SER154 and two 'O' atoms in the demethoxycurcumin aromatic ring showed a strong interaction pattern. The 'O' atom in both sides of the demethoxycurcumin aromatic ring was more potent in p38α inhibition.

Table I

Binding energy score of curcumin and its derivatives with p38α.

Name	CID	ΔG kcal/mol
Curcumin	969516	−6.67
Bisdemethoxycurcumin	5315472	−6.38
Curcuminglucuronide	71315012	−3.15
Cyclocurcumin	69879809	−6.12
Demethoxycurcumin	5469424	−5.88
Dihydrocurcumin	10429233	−4.48
Hexahydocurcumin	5318039	−4.81
Tetrahydrocurcumin	124072	−3.33

CID, PubChem compound identifier.

Figure 2

Binding modes of the top four curcumin-based compounds with the active site of p38α.

Table II

Autodock analysis of the top four curcumin-based compounds.

Complexes	Ligand binding pocket	Interaction
Curcumin-p38α	GLY110, ASN115, ASP112, ASP168, VAL38, LEU104, GLU71, LYS53, VAL105, LEU75, LEU86, VAL30, THR106, ALA111	Nil
Bisdemethoxycurcumin-p38α	LEU75, LEU74, THR106, ASP168, GLY170, GLU71, ILE84, PHE169, LEU167, MET109, GLY110, ASP112, ALA111, ALA157	OneBisdemethoxycurcumin: H39 - A:PHE169:O 2.36 A°
Cyclocurcumin-p38α	ILE84, LEU75, LEU86, LYS53, VAL105, GLY85, THR106, LEU104, VAL52, GLU71, ASP168, ARG67, TYR35, PHE169	TwoCyclocurcumin:H40 - A:LEU104:O 2.13 A°A:VAL105:N - Cyclocurcumin:O3 3.02 A°
Demethoxycurcumin-p38α	THR106, LYS53, LEU108, ILE84, LEU104, MET109, LEU167, VAL38, ASP112, HIS107, ASP168, TYR35, ASN155, SER154, GLU71, LEU75	ThreeDemethoxycurcumin:H39 - A:SER154:OG 1.92 A°ASN155:ND2 - :Demethoxycurcumin:O1 3.04 A°A:SER154:OG - :Demethoxycurcumin:O2 2.67 A°

Figure 3

Two-dimensional docked conformation of (A) curcumin (B) bisdemethoxycurcumin (C) cyclocurcumin and (D) demethoxycurcumin with the ATP binding domain of p38α.

Figure 4

Mapping of three compounds with their respective hits; green dotted lines represent hydrogen bonds. (A) Bisdemethoxycurcumin interacting with single amino acid PHE169. (B) The two hydrogen bond interactions between cyclocurcumin and LEU104 and VAL105. (C) Demethoxycurcumin forming three hydrogen bonds with SER154 and ASN155.

Molecular dynamic simulation analysis

After identification of the top 4 curcumin-based compounds from the molecular docking, molecular dynamic simulation (MDS) was employed to refine the top 4 complexes. MDS is a powerful technique to monitor protein-ligand complexes over a period of time. MDS run at a 100-nsec production was carried out under GROMOS force field. The analysis was evaluated by group properties, RMSD, hydrogen bond interactions and advance properties including PCA and MM/PBSA free energy calculations using trajectory files. All the 4 complexes were stable throughout the run. For each complex, RMSD was evaluated for the convergence of the protein structure towards an equilibrium state after lead binding. Fig. 5 shows the RMSD plot based on backbone atoms. Curcumin and bisdemethoxycurcumin in complex with p38α exhibited a jump at 30 nsec and equilibrated well and both of the complexes converged with an RMSD range almost >0.2 nm. The RMSD value of the bisdemethoxycurcumin-p38α complex was least among the complexes. With an RMSD of 0.3 nm, cyclocurcumin showed the highest deviation. From the RMSD plot of the demethoxycurcumin-p38α complex, the structure equilibrated at 30 nsec and then started to converge near 0.25 nm. The change in RMSD in all 4 complexes up to 30 nsec was due to ligand binding with the active site and after 30 nsec all of the structure obtained stability with less structural deviation at the end of the 100-nsec simulation. After ligand binding, the structural deviation of the bisdemethoxycurcumin-p38α complex was the most stable among the 4 under consideration. Inter-hydrogen bond interactions between protein and ligands were evaluated for the p38α-curcumin-based complexes (Fig. 6). In the case of the curcumin complex (Fig. 6A), NH plot results showed a range of 4–6 hydrogen bond interactions observed throughout the 100-nsec simulation and a maximum of 8 hydrogen bonds. NH analysis confirmed strong inhibition of p38α by curcumin in dynamic system similar to the docking results inferred with 6 hydrogen bonds. In the case of the bisdemethoxycurcumin complex (Fig. 6B), the results showed a range of 3–4 hydrogen bond interaction found throughout the 100-nsec simulation and a maximum of 6 hydrogen bonds similar to the docking results. In the case of the cyclocurcumin complex (Fig. 6C), the results showed a range of 1–5 hydrogen bond interactions found throughout the 100-nsec simulation and a maximum of 5 hydrogen bonds similar to the docking results. The inter-hydrogen bond interaction pattern suggested the plausible mode of the strong binding of the lead candidates with p38α which favored the inhibition mechanism.

Figure 5

Root mean square deviation (RMSD) of the top four p38α-cur-cumin-based complexes calculated over the backbone atoms as a function of time with respect to the starting conformation.

Figure 6

Hydrogen bond pattern over time for the three physically interacting curcumin-based compounds with p38α. (A) curcumin, (B) bisdemethoxycurcumin and (C) cyclocurcumin.

MM/PBSA calculation

The best 4 ligands in complex with p38α from MDS were used to calculate the binding free energy using MM-PBSA method. Snapshots were extracted at every 10-psec of stable intervals from 100-nsec MD trajectory and served as input for calculation. The binding free energy and its corresponding components obtained from the MM/PBSA calculation of the p38α curcumin-based complexes are listed. The results indicated that 4 curcumin-based compounds possessed a negative binding free energy of −141.727, −132.758, −151.455 and −133.56 kcal/mol (Table III). Moreover, van der Waals, electrostatic interactions and non-polar solvation energy negatively contribute to the total interaction energy while only polar solvation energy positively contributes to total free binding energy. The relative binding free energies of the 4 p38α-curcumin-based complexes supported the strong binding in the dynamic system. The interaction with the active site of p-p38α is shown in Fig. 7.

Table III

MM-PBSA score of the top four curcumin-based compounds.

Complex	ΔGbind (kcal/mol)
Curcumin-p38α	−141.727
Bisdemethoxycurcumin-p38α	−132.758
Cyclocurcumin-p38α	−151.455
Demethoxycurcumin-p38α	−133.56

MM-PBSA, molecular mechanics Poisson-Boltzmann surface area.

Figure 7

Residue wise calculation of the binding mode interaction between the top 4 curcumin-based compounds and the active site of p38α.

Principal component analysis or essential dynamics was performed on the cyclocurcumin complex with p38α. The analysis was used to monitor the overall strenuous motions of the complex. The results showed that the first 6 eigenvectors formed 90% of the total motion of the complex. The first 2 components of the eigenvectors were projected at 300 K and analysis of this revealed the clusters of stable states for the complex (41). The FEL was plotted over this projection, which provided a clear description of the stability of the complex (Fig. 8).

Figure 8

Projection of the motion along the first and second principal eigenvectors in phase space plotted against each other for the cyclocurcumin-p38α complex.

Cyclocurcumin inhibits the release of TNF-α from human macrophages

TNF-α is a key factor in a variety of inflammatory diseases. Keeping in view the role of p38α in regulating TNF-α expression, we evaluated its expression in the culture supernatants of human macrophages (which constitute the major TNF-α-producing cells in highly relevant inflammatory disorders) treated with non-cytotoxic doses of cyclocurcumin. In order to determine the suitable dosage of cyclocurcumin that result in less cytotoxicity, macrophages (obtained with treatment of THP-1 monocytes with CSF), were treated with different concentrations of cyclocurcumin (10–100 µM) for 24 h, and the viability was measured by MTT assay. The results demonstrated that treatment with cyclocurcumin at concentrations such as 10, 20 and 40 µM exhibited the least cytotoxic effects (Fig. 9). Therefore, these concentrations of cyclocurcumin, without cytotoxic effects on macrophages were used. To test the effect of cyclocurcumin, macrophages were treated with the indicated concentrations of cyclocurcumin or 10 µm SB203580, a selective inhibitor of p38 MAPK. This was followed by stimulation with 10 ng/ml LPS for 18 h to induce the production of TNF-α. It was observed that treatment of the LPS-stimulated macrophages with cyclocurcumin led to significant inhibition in the release of TNF-α in a dose-dependent manner, similar to the inhibition caused by SB203580 (Fig. 10). Thus, these data confirm the role of cyclocurcumin in overcoming p38α-induced production of TNF-α and hence can be used as a therapeutic agent to target RA.

Figure 9

Evaluation of the least cytotoxic doses of cyclocurcumin on macrophages by MTT assay. Results are expressed in terms of percentage inhibition of cell proliferation. Data are expressed as mean ± SD of results obtained from at least three independent experiments. Statistical analysis by one-way ANOVA, **P≤0.01 and ***P≤0.001.

Figure 10

Effect of cyclocurcumin on tumor necrosis factor-α (TNF-α) production in macrophages. Cells were treated with the indicated concentrations of cyclocurcumin or 10 µm SB203580 or left untreated and then stimulated with 10 ng/ml lipopolysaccharide (LPS) for 18 h. The culture supernatants were then removed and assayed for TNF-α by enzyme-linked immunosorbent assay (ELISA). Data are expressed as mean ± SD of results obtained from at least three independent experiments. **P<0.01 and ***P<0.001 compared to the control.

In conclusion, natural products have been described as a gold mine for arthritis treatment. These products can suppress the expression of TNF-α, IL-1β, cyclooxygenase-2, lipooxygenase, matrix metalloproteinases or adhesion molecules, or suppress the activation of nuclear factor-κB (NF-κB). All are factors linked to the treatment of RA (42). Curcumin has long been reported to have an antirheumatoid effect. Initial studies have linked it with the downregulation of NF-κB (43). A plethora of research on curcumin has given us an opportunity to investigate it and its derivatives as potential kinase inhibitors for p38α. To the best of our knowledge, our study is the first to use an advanced computational approach to identify the possible DMRDs. MM-PBSA is a potent tool in the field of drug design (44) and this approach was employed in this study for short listing of compounds. From the molecular docking and molecular dynamic simulation studies, the potential binding mode of cyclocurcumin with p38α with stability was revealed as a top compound for the treatment of RA. Finally, inhibition of the release of TNF-α from LPS-stimulated macrophages by cyclocurcumin treatment confirms its role as a potent p38α inhibitor.