Molecular simulations of doxorubicin complexed with native and modified cyclodextrins in water

In this study, all-atom molecular dynamics simulations have been performed to investigate from a nanoscopic point of view the interactions in supramolecular complexes formed between doxorubicin and cyclodextrins in an aqueous phase at a temperature of 310 K. All the simulations are on the time scales of microseconds and are carried out on High Performance Computing resources provided by the Greek Supercomputer ARIS. Doxorubicin, which is a well-known anticancer drug belonging to the class of anthracyclines, is simulated with two cyclodextrins of broad utilization in drug delivery, i.e.: γ-cyclodextrin and 2-hydroxypropyl-β-cyclodextrin. In both cases, the complexation process takes place spontaneously during the course of the simulations and a detailed analysis of the formed complexes is conducted in terms of structural, dynamical and thermodynamic properties. In addition, we present a quantitative analysis of how doxorubicin affects, upon complexation, the geometry, the hydrogen bonding network and the hydration of the selected cyclodextrins. Special attention is paid to the thermodynamic description of the simulated systems. More specifically, the standard binding free energy is calculated with an extended version of the linear interaction energy method. The presented in silico approach reveals that doxorubicin forms stable complexes with the chosen cyclodextrin molecules, there is a water solubility enhancement of the drug with respect to its unbound state and also cyclodextrins act protectively to the drug preventing undesired effects, such as the glycosidic degradation. Excellent agreement is found between the presented results and published in silico and experimental findings.

Cyclodextrins (CDs) are cyclic oligosaccharides that are composed of D-glucopyranoside units [17].Their shape is a truncated cone displaying an inner hollow hydrophobic core, whereas the outer surface is hydrophilic [17].This special geometry makes CDs excellent hosts for the formation of supramolecular complexes with either hydrophobic or amphiphilic drug molecules and improving their pharmacological abilities by enclosing them partly or fully into their hydrophobic cavity.There exists a rich literature concerning drug delivery of DOX or DOX derivatives with native or modified CDs [2][3][4][5]7,[10][11][12][13][18][19][20][21][22].Among the native CDs, γ-CD forms the strongest supramolecular complexes with DOX due to its larger cavity size [7,13,18,19,23], while the very low toxicity of this cyclodextrin [13] is another significant advantage over other molecules of this class.In the framework of this article, two cyclodextrins are considered, i.e.: γ-CD and 2-hydroxy-propyl-β-cyclodextrin (HP-β-CD); the latter is utilized widely for drug delivery purposes.To be more specific, HP-β-CD, which is of very low toxicity as well [24] and may form stable complexes with drug molecules [5,12,[20][21][22].

Computational Methods and Simulations Details
In this article, atomistic molecular dynamics simulations in the isothermal-isobaric statistical ensemble are carried out to investigate from a nanoscopic point of view the interactions in supramolecular complexes formed between DOX and the selected CDs (γ-CD, HP-β-CD) in water phases at a temperature of 310 K and at a pressure of 1 bar, so as to mimic human body conditions.All the simulations are on the time scales of microseconds and are carried out by making use of the GROMACS [25] package on High Performance Computing resources provided by the Greek Supercomputer ARIS.The selected force fields and the simulation protocol adopted herein are exactly the same as those of reference [26] (see Section 2 therein) which concerns noncovalent complexes formed between a wellknown antidepressant drug and β-CD.The initial structures of DOX and γ-CD are obtained from the RSCB protein data bank (PDB IDs: 1I1E [27] and 2ZYK [28], respectively).As far as the HP-β-CD molecule is concerned, the initial configuration is generated from the β-CD by replacing manually the primary hydroxyl groups with 2-hydroxy-propyl substituents via the Maestro Schrödinger software [29].
The chemical structure of DOX is depicted in Fig. 1 which is created with the ChemDraw online platform [30].As seen from the aforementioned figure, DOX is simulated in the protonated as well as in the neutral state, which henceforward are referred to as pDOX and nDOX, respectively; it is pointed out that both pDOX and nDOX are present in the normal blood pH range.In addition, a specific group of atoms that are necessary for the proper discussion of the results in the next section are marked in red color in Fig. 1.The simulation boxes of these complex systems are constructed via the Packmol software [31].In all simulations, the constituent parts of the complexes are initially placed in the water phase in the unbound state so as to investigate whether the complexation is a spontaneous process.In the case of pDOX, a chloride ion is added to neutralize the simulated system.The molecular visualization images are generated employing VMD [32].Regarding the Solvent Accessible Surface Area (SASA) calculations [33], they are realized by making use of the sasa utility of GROMACS [34,35].
The binding free energy is computed via an extended version of the Linear Interaction Energy (LIE) method [36], as proposed by Montalvo-Acosta et al. [37], which calculates binding free energies with statistical errors less than 6.4 kJ/mol with respect to the corresponding experimental values.The binding Gibbs energy in the context of this version of the LIE method is calculated by the equation ( 1): (1) the first two terms of equation ( 1) are explained in detail elsewhere [26,37].The third term stands for the host's strain energy, upon complexation, and the procedure for computing it is analytically explained in the original work of Montalvo-Acosta et al. [37].In order to calculate this quantity, a minimization protocol and a simple molecular mechanics energetic calculation are needed.To this end, we employed the sander and cpptraj modules of the AMBER14 package [38,39].

Results and Discussion
The analysis of results begins with a study of complex formation by monitoring the center of mass distances (COM-to-COM distances) between the constituent molecules of the complexes as a function of the simulation time.This approach was also used successfully elsewhere [26,40].The COM-to-COM distance plots are depicted in Fig. 2, from which it is concluded that there are three distinct regions in each plot.The first one is the region in which the two molecules are in their unbound state and the COMto-COM distances are big (>2 nm).The second region has a duration on the order of 10-200 ns and corresponds to the external complexes in which the drug is bound to the outer surface of the CD (the COM-to-COM distance is nearly 1 nm); such complexes are considered to be in a metastable state.In the third region the two molecules have formed a stable inclusion complex and their COM-to-COM distances take values less than 0.5 nm.Representative visualizations of the three regions are provided in Fig. 2 as well.These structures were obtained via cluster analysis over nearly 25,000 time-frames in GROMACS with the RMSD-cutoff value being 0.15 nm.In all the complexes, we observe that the main hydrophobic parts of DOX have entered the hydrophobic cavity of CDs in order to avoid water and to interact favorably with the CDs' walls mainly via van der Waals forces, whereas the hydrophilic parts lie outside the cavity, interacting favorably with the hydrophilic outer CD's surface and water, such as the hydrogen bonds.As far as the latter interactions are concerned, some results are presented in Table 1 in which the number of hydrogen bonds (HBs) formed between various pairs of molecules are provided.We can clearly see from Table 1 that the unbound pDOX forms more HBs with water than nDOX, revealing the more hydrophilic character of pDOX, due to the net charge of its amine group (this group contributes 0.680 ± 0.007 more HBs with water in the protonated case).In all the cases, we can easily see that DOX loses HBs with water, upon complexation and this is due to the fact that when it enters into the cavity, water molecules are stereochemically hindered, by CD, to reach near DOX's HBcontributing groups.CDs also lose HBs with water, upon complexation and the reason is that DOX displaces the water molecules which were initially entrapped into CD's cavity and were interacting via HBs with cyclodextrin.Due to the size of DOX molecule, it perturbs the external water shell too (mainly in the HP-β-CD complexes), albeit to a much smaller extent.The above HB analysis shows that the HB network between DOX and CD, as well as water is one of the driving forces of the complexation, a finding which is in accordance with the literature [17,26,41].For further validation of our calculations, we mention that the unbound γ-CD formed 36.45 HBs with water which is in excellent agreement with other in-silico studies (36)(37) [42,43].Very good agreement is also obtained in the case of the intramolecular HBs in the unbound γ-CD molecule; the calculated value from the simulation 5.97 ± 0.06 is close enough to the value 6.4 reported in reference [42].It is known that the CDs entrap water molecules into their cavities (named cavity waters), the release of which is an important driving force for the CD complexations [17,26,41].This is analyzed in this paragraph with the selected Radial Distribution Function (RDF) plots which describe the distribution of specific atoms or molecules around a reference point.Fig. 3 displays the RDFs of water molecules around γ-CD's COM in both the unbound and bound states.Fig. 4 displays the RDFs around the COM of the cavity in a HP-β-CD in both the unbound and bound states.The number of cavity waters may be extracted by integrating the RDF, (), over a sphere with its center being the RDF's reference point and having an equivalent radius,   , as shown in equation ( 2): (2) where  ∞ is the bulk numerical density of water at the temperature and the pressure of the simulation.The equivalent radius,   , is taken here as the half of the CD's internal diameter [44], which has been used elsewhere too [26,40].The RDF around γ-CD in the unbound state (Fig. 3a) is similar with the prediction from other computational works [42,45,46], confirming the validity of our results.In both CD molecules, studied here, we observe similar RDFs qualitatively in their unbound states in water (Fig. 3a for γ-CD and Fig. 4a for HP-β-CD).An important point on these plots is that both of them have a peak at small distances (0.26 nm for γ-CD and 0.21 nm for HP-β-CD approximately), indicating that there are indeed cavity waters.Then, we find a minimum (0.62 nm for γ-CD and 0.50 nm for HP-β-CD approximately), due to the CD's walls.Lastly, we have the external hydration shell, giving a maximum and a minimum to the plots.Using equation (2), we find 11.9 and 3.0 water molecules inside the cavities (γ-CD and HP-β-CD, respectively), which are close to results from other computational works [42,46,47].In all cases, we observe that the first maxima and minima have nearly vanished in the bound states (Fig. 3b, 4b), because the water molecules have been excluded from the CDs, due to the drug's insertions.The corresponding cavity waters are now 0.3, 1.0, 0.3, 0.0 for the nDOX/γ-CD, pDOX/γ-CD, nDOX/HP-β-CD and pDOX/HP-β-CD cases, respectively.These findings confirm the above-mentioned bibliographic fact that the cavity waters are extracted, upon complexation.Next, we perform a SASA analysis which may reveal whether there is a solubility enhancement of the selected drug upon complexation [26,48].These results are summarized in Table 2.It is worthmentioning that the SASA for the unbound γ-CD is close to the one calculated in reference [49], where similar SASA calculated are conducted.The DOX's SASA increased dramatically from approximately 7.7 nm 2 in its free state to 13.6 nm 2 and 16.6 nm 2 when it is complexed with γ-CD and HP-β-CD, respectively.The SASA calculations may provide insight about the parts of each molecule buried, upon complexation.It has been mentioned in the Introduction section that DOX have the serious problem of being unstable in a water solution and having cardiotoxic side-effects.In the literature, these undesired effects seem to come from the glycosidic bond, because of the glycosidic degradation and the anthraquinone group, the reactivity of which causes instability in water and cardiotoxicity [6,9,14].It is known bibliographically that CD complexation is able to resolve these problems to some extent and this is confirmed here too by an additional SASA analysis.We can see that the glycosidic oxygen atom has been buried into the cavity of CDs by (93 ± 1)%, (98.9 ± 0.2)%, (98.1 ± 0.5)% and (98.4 ± 0.9)% and the anthraquinone group has been buried by (55.4 ± 0.2)%, (45.8 ± 0.2)%, (86 ± 1)% and (84.7 ± 0.8)% in the nDOX/γ-CD, pDOX/γ-CD, nDOX/HP-β-CD and pDOX/HP-β-CD respectively.These mean that the above-mentioned parts of DOX are covered and protected from the CDs.As already mentioned, we pay attention to the thermodynamic description of the systems under consideration herein.The binding Gibbs energies, as calculated by equation ( 1), are −19.8 ± 0.4, −21.2 ± 0.8, −20.2 ± 3.8 and −24.9 ± 0.8 kJ/mol for the nDOX/γ-CD, pDOX/γ-CD, nDOX/HP-β-CD and pDOX/HP-β-CD complexes respectively.For validation purposes, we mention that UV-vis spectroscopy and isothermal titration calorimetry experiments report binding Gibbs energies −22.0 and −22.8 kJ/mol, respectively (310 K) for the pDOX/γ-CD complex [13], while fluorescence and Lineweaver-Burk analysis place this value at −25.0 and −27.4 kJ/mol (297 K) for the pDOX/HP-β-CD complex [4].Both of these experimental works are in agreement with our results.One last note here is that in both cases, the protonated complexes give more negative binding free energies, compared to the corresponding neutral ones.This finding is in agreement with experimental work [18], in which the association constant of DOX/γ-CD has been measured at different pHs and the maximum affinity was found to be at relatively neutral pHs, rather than extremely acidic or basic ones.This was explained by the mono-cationic state of DOX at these pHs, which is something that we can clearly see too.
It is bibliographically known that a guest molecule is able to induce conformational changes to the CD molecule, upon complexation [41].A good measure for the geometric effect of the guest molecule to the CD's structure is the gyration radius [50].The results indeed indicate some minor conformational changes to CDs, upon the complexation.In the γ-CD complexes, we saw that the gyration radius decreased by (1.2 ± 0.1)% and (0.3 ± 0.1)% in the neutral and the protonated cases respectively, which means that the γ-CD's structure converges a little in both cases, in order to maximize the interactions with the drug.The conformational changes were more important in the HP-β-CD case, in which we saw an increase of (1.8 ± 0.4)% and (2.2 ± 0.4)% in the cavity in the neutral and protonated cases respectively and an increase of (11 ± 3)% and (8 ± 1)% in the substituents in the neutral and protonated cases respectively.The cause of this is that HP-β-CD's internal spaces are smaller than the γ-CD's ones and its sizes are comparable with DOX's size.So, DOX's insertion creates repulsions between itself and the CD's walls, so the walls move away from the drug, resulting in having a wider distribution, thus a greater gyration radius too.Also, substituents show significantly greater changes and corresponding errors and this is due to their greater flexibility, compared with the main cavity.These findings can easily be related to the fact that DOX barely forms complexes with α-and β-CD or if it does, the affinity is smaller, due to the smaller host cavity [7,18,19,23].Lastly, we confirm our results on the gyration radii by seeing the agreement with the computational works [2,51], in which the gyration radii of γ-CD and HP-β-CD, respectively have been calculated.

Conclusions
DOX, in the neutral as well as in the protonated state, forms stable inclusion complexes with γ-CD and HP-β-CD.In all cases, the complexation process takes place spontaneously, during the unbiased molecular dynamics simulations and this is confirmed by the negative value of the binding Gibbs energies.The analysis of the previous section shed light into several aspects of complexation mechanism and the underlying driving forces in which are included the hydrophobic interactions, the formation of an HB-network between the molecules and water too, the exclusion of cavity waters and of course the nonbonded interactions, such as the electrostatic and Lennard-Jones ones, since these interactions gave the negative value to the binding free energy, according to the LIE method.The complexations improve some important pharmaceutical characteristics of DOX, such as water-solubility.They also, protect parts of the DOX molecule which are responsible for some serious side-effects of this drug in its action, namely the glycosidic bond that causes degradation in water and the anthraquinone which causes photolytic degradation and cardiotoxicity.A last conclusion is that the complexation process affects CD's geometrical conformation as well, compared with the one that it adopts in the unbound states in water solution.We observed that γ-CD's walls came closer to each other, resulting in a more compact conformation, while HP-β-CD exhibited a wider conformation, due to the similarity between the sizes of the drug and the cavity.Our results are in agreement with experimental findings, making the computational models and algorithms used here, a valuable tool for the simulation of similar systems.

Figure 1 .
Figure 1.2-D representation sketches of DOX in its (a) neutral and (b) protonated states.

Figure 2 .
Figure 2. COM-to-COM distance between the molecules at the (a) DOX+γ-CD and (b) DOX+HP-β-CD cases.Among the inclusion complexes, the left ones are the neutral complexes, while the right ones are the protonated complexes.

Figure 3 .
Figure 3. RDFs of water molecules around the γ-CD's COM in the (a) unbound state and (b) bound state.A water molecule is represented by its oxygen atom.

Figure 4 .
Figure 4. RDFs of water molecules around the HP-β-CD molecule's cavity in the (a) unbound state and (b) bound state.A water molecule is represented by its oxygen atom.

Table 1 .
Number of Hydrogen bonds (HBs) in the unbound and bound states for the pairs watermolecule and CD-DOX.

Table 2 .
SASA calculations of DOX and CDs in their unbound and bound states.