Cyclodextrins are widely used for the solubilisation of poorly soluble drugs in the formulations. However, current cyclodextrin formulation development strongly depends on trial-and-error in the laboratory, which is time-consuming and high cost. The aim of this research was to compare three modeling approaches (Docking, molecular dynamics (MD), and quantum mechanics (QM)) for cyclodextrin/drug complexation. Ibuprofen was used as a model drug. Binding free energy from three simulation methods was calculated as an important parameter to compare with the experimental results. Docking results from AutoDock Vina program showed that the scoring of complexation capability between ibuprofen and cyclodextrins is alpha (α), gamma (γ), beta (β), and HP-beta-cyclodextrins, which indicated similar ranking with the results from phase, solubility diagram experiments. MD simulation indicated that ibuprofen could form the stable complexes with β-, γ-, and HP-β-cyclodextrins, but not for alpha cyclodextrin. Binding free energies from the MD simulation for β-, γ-, and HP-β-cyclodextrins were −3.67, −0.67, and −3.87 kcal/mol, individually. The enthalpies of QM simulation for β-, γ-, and HP-β-cyclodextrins were −17.22, −14.75, and −20.28 kcal/mol, respectively. Results from all three modeling approaches showed similar ranking between ibuprofen and four cyclodextrin molecules as the experimental data. However, MD simulation with entropy calculation had the closest value to experimental data for β and HP-beta-cyclodextrins. Thus, MD simulation with MM-PBSA method may be fit to in silico screen for cyclodextrin formulations.

1. Introduction

Cyclodextrins (CDs) are natural pharmaceutical ingredients for the solubility-enhancement of poorly soluble drugs for over half century. Chemical composition of cyclodextrin consists of D-glucopyranoses, which are linked by alpha 1,4-glycosidic bonds to form the macrocyclic ring conformations. Cyclodextrins are classified into alpha, beta, and gamma cyclodextrins, comprised of six, seven, and eight D-glucopyranose residues, respectively [1]. These macrocyclic structures have hydrophilic outer with hydroxyl groups, while the nonpolar inner cavity is hydrophobic [1]. The unique structure enables CDs to accommodate poorly water-soluble drugs to form drug-inclusion complexes. However, cyclodextrins are poor across most of the biological membranes due to their high molecular weight and inherent hydrophilicity [2]. There are over 35 commercially available pharmaceutical formulations including various natural or chemically modified cyclodextrins [1]. Some experimental techniques are able to characterize drug-cyclodextrins inclusion complexes, such as thermoanalytical [3], spectroscopic [4], microscopic [5], and chromatographic techniques [6]. However, current cyclodextrin formulation development mainly depends on trial-and-error by formulation scientists in the laboratory, which is time-consuming and high cost.

Molecule modeling method is able to investigate the structure, dynamics, and energetics of cyclodextrin systems. From the search of Web of Science, there are over 700 research papers about theoretical calculations on cyclodextrins complexation with various drugs until 2014. These theoretical calculations could be classified into four different types, such as quantum mechanics (QM) [7], molecular dynamics (MD) [8], docking [9], and quantitative structure activity relationships (QSARs) [10]. However, the results from different modeling and experimental approaches were conflicting and contradictory. For example, ibuprofen and cyclodextrin complexes have been studied using different experimental techniques [11, 12]. There are also a few theoretical research papers on the complexation of ibuprofen and cyclodextrins [13, 14]. However, these experimental and theoretical studies obtained various results. Some recent research found the binding free energy between ibuprofen and β-CD was stronger than complex of ibuprofen and HP-β-CD [11, 12]. However, Mura (1998) found the HP-β-CD solubilizing effectiveness for ibuprofen was stronger than β-CD in aqueous environment [12]. Núñez-Agüero (2006) used MM-GBSA method to calculate the enthalpy between ibuprofen and β-CD [15]. The results from recent researches were summarized in Table 1.

Thus, it is necessary to establish a standard method to investigate the complexation between cyclodextrin and drug. The purpose of this study was to compare different simulation methods of ibuprofen and cyclodextrin complexation. Ibuprofen was chosen as a model drug.

2. Simulation Methods

2.1. Model Building of Drug and Cyclodextrins

The crystal structures of alpha (α), beta (β), gamma (γ) cyclodextrins were taken from The Cambridge Crystallographic Data Centre (CCDC) (α-CD for “BAJJAX,” β-CD for “ARUXIU,” and γ-CD for “CYDXPL”). Hydroxypropyl-beta-CD (HP-β-CD) was modified from the structure of β-CD using Discovery Studio Visualizer. Ibuprofen was built by Discovery Studio Visualizer. The structure of all molecules was optimized with a fast, Dreiding-like force field by Discovery Studio Visualizer.

2.2. Docking Simulation Setting

AutoDock Vina was one of widely used docking programs, which can predict free energy of ligand-receptor interaction [20]. In our research, the AutoDock Tools package (version 1.5.6) and AutoDock Vina (version 1.1.2) were used to perform docking studies [9].

In docking simulation, ibuprofen was used as a ligand, while cyclodextrins were used as receptor. Before the docking, AutoDock Tools was used to optimize the structure of ibuprofen from the PDB files. Atom C2 of ibuprofen was chosen to be the root atom because all its bonds to other atoms were rotatable. AutoDock Tools was used to determine the possible rotations of bonds. Current total five rotatable bonds were the bonds between atom C8 and C7, C7 and C5, C2 and C10, C10 and C11, and C11 and O12 as shown in Figure 1. The pdbqt format files (required as input to AutoDock Vina) were generated using AutoDock Tools. During the docking procedure, ibuprofen was docked to cyclodextrin molecules using the following Cartesian coordinates as the center of the search space:  Å,  Å, and  Å. The dimensions of docking grid were 60 Å × 60 Å × 60 Å. All the other parameters were the default values by AutoDock Vina. The docking pose with the highest score was treated as the stable binding model. The docking score from AutoDock Vina output was used as the binding free energy.

2.3. MD Simulation

The stable binding model from docking was used as the starting structure of MD simulation. The MD simulation was performed by the Amber14 and Amber Tools 14 software package. Currently, there is no specific force field for cyclodextrin molecules. Unlike the traditional Amber force fields, general Amber force field (GAFF) is more general for organic chemical [21]. Some recent articles chose GAFF force fields in MD simulation of CDs [2224]. So the GAFF force field was selected in Antechamber module with AM1-BCC charge method for all molecules. A water cube of 10 Å thickness was added to solvate the system with the TIP3P water model. All MD setting details were shown in Table 2. After two stages of minimization, 30 ns MD simulations were performed. The detailed simulation protocol was similar to our previous simulation researches [2527].

2.4. QM/MM Simulation

The stable binding model from docking was used as the starting structure of QM simulation. Coupled potential QM/MM simulation has been performed in Amber Tools 14. The complex was masked as quantum mechanics (QM) modeling by the semiempirical PM3 Hamiltonian method, while all water molecules were treated as classical MM. The periodic boundaries and PME approaches were used to calculate the long-range QM-MM electrostatic energies and forces and the long-range QM-QM forces with 8 Å cutoff distance in both two stages. The SHAKE algorithm was applied to QM atoms in order to constrain all bonds involving hydrogen.

In case of ibuprofen and cyclodextrin simulation in water, the system was optimized with the 200-steps steepest descent method and then 300-steps conjugate gradient method. Then the system was heated up to 300 K. In QM/MM stage, the system carried out 1 ns simulation at 300 K with an average pressure of 1 atm by using Langevin dynamics with the collision frequency 1 ps−1.

2.5. MM-PBSA Method for Binding Free Energy Calculation

MM-PBSA method was chosen to calculate the binding free energy [28]. In this study, the binding free energy of the ibuprofen and CD inclusion complex () was calculated by the free energy of complex () and the isolated ibuprofen () and CDs () as the following equation:

The Gibbs free energy () was calculated by the enthalpy () and entropy with invariable temperature ():

2.6. Phase-Solubility Diagram

Over amount of ibuprofen was added in the 10 mL of cyclodextrin solution with the concentration of 0, 0.0025 M, 0.005 M, 0.01 M, 0.02 M, and 0.03 M, respectively. The solutions were sonicated for 60 minutes to reach the equilibration. After the samples were filtered, the filtrate was measured by UV-visual Spectrometer at the wavelength of 264 nm. The experimental result was shown in Figure 2.

The stability constants of ibuprofen and CDs complexes in water were calculated by the following equation:

3. Result and Discussion

3.1. Docking Result

The molecular docking is an attractive strategy for ligand-receptor interaction. In our study docking was conducted to calculate the binding energies between ibuprofen and cyclodextrins complexes by AutoDock Vina. Figure 3 showed the optimal models for ibuprofen and CDs. Unlike the interaction of α-CD/ibuprofen complex, the carboxylic part of ibuprofen inserted deeply into the cavity of β- and γ-CD molecules due to their large cavity. Moreover, ibuprofen molecule had a slight inclination in the cavity of γ-CD because γ-CD has the largest cavity.

According to the docking results, the binding energies between ibuprofen and cyclodextrins were −3.6, −4.4, −4.2, and −5.0 kcal/mol, respectively. The binding energy of ibuprofen and α-CD was the lowest value among four complexes, which confirmed the shallow structure among them. Ibuprofen and HP-β-CD complex had the strongest binding energy. The possible reason is that the interactions between hydroxypropyl group of HP-β-CD and isopropyl group of ibuprofen make the complex more stable. The negative binding affinity energy values suggested that the formation of all these inclusion complexes was a spontaneous reaction.

3.2. MD Trajectory Analysis

Figure 4 showed the snapshots during the 30 ns MD simulations of ibuprofen and α-CD. These snapshots indicated the unstable state of the complex. At the first 6 ns, the structure of ibuprofen and α-CD complex did not change much. After that, the ibuprofen started to reverse its structure. After 24 ns, the ibuprofen separated from α-CD. However, the complex of ibuprofen and α-CD was found to be unstable during another three MD simulations. The possible reason may be the weak binding energy of ibuprofen/α-CD complex, which is also confirmed by the low solubilisation of α-CD to ibuprofen in the experiments.

Figures 57 showed the stable structures between ibuprofen and other three cyclodextrin molecules during the 30 ns simulation. In the simulation, the isopropyl group of ibuprofen totally inserted into the cavity of β-CD and the carboxylic part will stick out from the primary side of β-CD, which was in agreement with previous study about crystal structure of ibuprofen/β-CD complex [29]. However, our result was different from previous simulation papers [15]. The possible reason was improper force field and too short simulation time (4.5 ns) in previous simulations. For γ-CD, ibuprofen lied in the cavity with about 45° inclination angle. In addition, some portion of the γ-CD molecule became twisted to fit the inclusion ibuprofen. This phenomenon was not found in docking studies because all CD molecules were treated as the rigid structure in the docking. Figure 7 displayed that ibuprofen binding with HP-β-CD on the opposite direction to β- and γ-CD. The possible reason was nonpolar interaction between hydroxypropyl group of HP-β-CD and isopropyl part of ibuprofen. The MD simulation predicted model was in agreement with previous research about the spectroscopic characterization of ibuprofen and 2-hydroxypropyl-β-CD inclusion complex [16].

Figure 8 showed the root mean square displacement (RMSD) for backbone heavy atoms of the complexes in 30 ns MD simulation. RMSD value of α-CD/ibuprofen complex increased sharply after 20 ns MD simulation, which agreed with the structural change in Figure 4. However, for the complexes ibuprofen and β-, γ-, and HP-β-CD, the RMSD values were relatively stable, which indicated the complexes had reached the equilibrated state after 30 ns simulation. Therefore, only the inclusion complex of β-, γ-, and HP-β-CD was further calculated for binding free energy and QM modeling.

The MM-PBSA approach is a useful tool to calculate the binding free energies. The 100 MD snapshots from the last 1 ns of each system were used for binding free energy calculations. The binding free energies () and the other energy contributions were shown in Table 3. The van der Waals force () made the key contribution to the complexation between ibuprofen and cyclodextrin molecules. The contribution of the HP-β-CD complex was stronger than those of β- and γ-CD complexes. The results of were not accurately equal to the binding free energy because the entropy contribution was not estimated. Thus, () were calculated for the Gibbs free energy of the complexes. The binding energies of β- and HP-β-CD complex are −3.67 and −3.87 kcal/mol, respectively, which were close to the experimental values (−3.58 and −3.78 kcal/mol) in Table 4. These results suggested that the entropy and solvation contribution were paramount for the binding free energy between ibuprofen and cyclodextrin molecules. MD simulation had very good prediction for ibuprofen with β-CD and HP-β-CD complex, but not for complex of γ-CD. The possible reason may be the binding mode of γ-CD and ibuprofen was not 1 : 1 because of the large cavity of γ-CD.

3.3. QM/MM Result

Due to the separation of α-CD and ibuprofen in MD simulation, we only calculated the QM/MM simulation for ibuprofen with β-, γ-, and HP-β-CD complexes. In QM/MM simulation, the binding modes of ibuprofen and CDs were shown in Figure 9, which were very similar to the models in MD simulation. The PM3 theory was used for quantum calculation. The QM/MM-GBSA method predicted the enthalpy of the ibuprofen with β-, γ-, and HP-β-CD complexes being comparatively close to the enthalpy values in MD simulation, as shown in Tables 3 and 4. However, QM calculation took much longer computing time than MD simulations and the obtained energy magnitudes were far off from the experimental values because the entropy was not included. Thus, QM simulation is not recommended for in silico formulation screening.

In this paper, three different kinds of computer modeling methods were used to investigate the complexation of ibuprofen with α-, β-, γ-, and HP-β-CD. The binding energies of AutoDock Vina results show stronger binding of HP-β-CD than other complexes. However, the values from docking were a little overestimated compared to experimental data. MM-PBSA approaches estimated that the VDW mainly contributed to the complexation formation. In addition, the entropy and solvation contribution were very important for the binding energy calculation, which was ignored in previous researches [15]. The binding free energy from MD simulation was very close to experimental results. QM/MM simulation predicted the similar enthalpy of ibuprofen with β-, γ-, and HP-β-CD with MD simulation, but QM simulation took significantly longer simulation time. In comparison to the experimental value, the values from docking results for all complexes were found a little overestimated, but the ranking between them were correct. MD simulation had very good prediction for ibuprofen with β-CD and HP-β-CD complex, but not for complex of γ-CD.

4. Conclusion

In this study, three modeling methods were used to investigate the complexation of ibuprofen and CDs in water. All modeling and experimental results showed that HP-β-CD inclusion was stronger than other complexes. Docking result has similar ranking with experiments. The Gibbs free energy for β-CD and HP-β-CD complex from MD simulations and MM-PBSA method were the most close to the experimental value. QM method took too expensive computing resource to use in the wet laboratory. MD simulation with MM-PBSA method may be fit to in silico screen for cyclodextrin formulations.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This study was supported by Science and Technology Development Fund (FDCT), Macao SAR, (074-2012-A3) and Research Fund of University of Macau (MRG013-WYT-2013-ICMS and CPG2014-00012-ICMS).