Molecular modeling of β -cyclodextrin inclusion complexes with pharmaceutical compounds

Ten molecularly diverse compounds, selected on the basis of the standard free energy of complexation with β -cyclodextrin, have been used as training set of molecular docking experiments. After the conformational search of each compound the Molecular Interaction Evaluation (MOLINE) program has been adopted to generate the inclusion geometries. The recognition process has been thermodynamically evaluated by the computational method, theoretical and experimental free energies of complexation have been compared obtaining an high correlation ( r 2 > 0.9). A statistical and graphical analysis of the binding inclusion geometries completed the work.


Introduction
Molecular recognition (MR) is a general term designating non-covalent interactions between two or more compounds belonging to host-guest, enzyme-inhibitor and/or drug-receptor complexes.A rigorous approach to the MR study should involve the adoption of a computational method independent to the chemical intuition of the researcher.Drug design purposes prompt another challenging feature of such an ideal computational method, the ability to make sufficiently accurate thermodynamic predictions about the recognition process.In this study we report an application of the computational method MOLINE (Molecular Interaction Evaluation) created for studying complexes with an unbiased fashion. 1 With respect to recently reported computational studies on cyclodextrin inclusion recognition mostly by molecular dynamics techniques, 2,3 MOLINE is based on a systematic, automatic and quasi-flexible docking approach that prevents the influence of the chemist's intuition in the configuration generation.The inclusion molecular system focused in this work is based on β-cyclodextrin as host (Figure 1) and a set of pharmaceutical organic molecules as guests.The aim of this work is to find the best computational condition able to reproduce the same trend of free energy interaction in this set of host-guest inclusion complexes.The selection of ten ligands (Figure 2) was carried out selecting quite strong interaction experimentally characterized by standard free energy of complexation (∆G c ) lower than -4.90 kcal/mol. 4This level of binding configures interaction comparable to drug-receptor complexes.All data were obtained in water environment with specific pH condition and 1:1 stoichiometry.Chemically 1 and 2 are acid compounds with ∆G c reported for pH equal to 7.20 and 9.94. 5,63, 4 and 5 are tricyclic antidepressants (nortriptylin, amitriptylin and doxepin) experimentally measured with CYD at neutral pH conditions. 76 is the chlorpromazine, one of most known antipsychotic agent, that was reported as chlorhydrate. 8Aliphatic neutral molecules are the compounds 7-9, characterized by ∆G c with CYD in the order of magnitude of -4.92±0.05kcal/mol. 9,10,11In the same range, but positively charged due the quaternary ammonium nitrogen, it is compound 10.

Computational procedures
The MOLINE docking approach required the conformational search of the interacting molecules.Following the experience of two previously reported works 13,14 the β-cyclodextrin was considered in "open conformations" with three structures symmetrically generated considering the gauche +, gaucheand trans rotamers of the χ 5 torsion angle.This dihedral angle is responsible of the high flexibility of primary hydroxyl sidechain.Compounds 1-9 were submitted to a 1000 iterations Monte Carlo (MC) search with MacroModel, 15 force field AMBER* united atoms and GB/SA water implicit model of solvatation. 16Compound 10 showed no rotatable bonds due to the symmetry of the trimethylammonium moiety.Most stable conformers obtained in this search with Boltzmann population higher than 90% at 300 °K were considered in the automatic docking procedure of the MOLINE program.According to our previous experience 13 the grid resolution G R and the Van der Waals compression factor χ were respectively fixed to 6 and 0.6 for a total of 65.712 configurations for each compound/cyclodextrin conformer combination.The force field used in the docking procedure was the same AMBER* united atom, a dielectric constant equal to 80 was considered for mimicking the water environment.After rigid optimization performed by three cycles of 100 Simplex iterations, the most stable configurations within 3 kcal/mol were submitted to a full energy multiminimization procedure with MacroModel in the same conformational search conditions.Duplicated structures were identified and removed by considering conformers within an energy window of 1 kcal/mol and RMS lower than 0.25 Å computed onto the entire complex atomic coordinates.Thermodynamic calculations were performed according to the MOLINE theory 1 in the full energy minimized ensemble considering the host-guest interaction energy only as reported in another molecular recognition study. 17A statistical analysis of the inclusion geometries was performed following the procedure of a previously reported work. 13Two dummy atoms averaging the fourteen CYD oxygen atoms in C2/C3 positions (center of the large entrance) and seven methylene C6 carbons (center of the small entrance) were fixed and used for determining two distances with respect to the centroid of the guest compounds 1-10.The classification of the binding mode was carried out considering the closest of the two distances.All the calculations were performed using a DS20 HP-Compaq workstation 666 MHz.

Results and Discussion
With the aim to reproduce the experimental trend of the free energy of complexation in 1:1 inclusion complexes of CYD and the guest compounds 1-10 the MOLINE docking approach was carried out.The MC search of the guest molecules revealed different conformational distribution related to the flexibility of each compounds.In table 1  The conformational distribution observed for compounds 1-10 is function of their flexibility, i.e. the number of rotatable bonds.The statistical descriptor used as reference for the conformational search quality is the average number of duplicates.In all cases this number is higher than 12, attesting that the conformational space of all compounds has been adequately explored by the MC protocol.In order to reduce the computational cost of the MOLINE procedure, only those conformers with Boltzmann population higher than 90 % at 300 °K were submitted to the docking computations.
The three symmetric "open" cyclodextrin conformations were docked using the protocol described in computational procedure paragraph with the selected conformers of the ten guest compounds.The interaction energies of most stable docking configurations obtained within 12.5 kcal/mol above the global minimum bimolecular complex were used to estimate the thermodynamics of the recognition process as summarized in table 2.  The correlation with the experimental free energy of complexation is reported in Figure 3.The computational method is not able to reproduce absolute free energy of complexation probably due to the limit of the molecular mechanics force field AMBER*.Quantitatively the miscalculation seems to be addressed to a systematic overestimation of the interaction energies as shown by the deep slope trend of the correlation equation y = 11.85 x + 37.59.The computational method does not predict the complexation free energy trend of compound 4 and 10, respectively over and underestimating their CYD interaction with respect to compounds 3 and 9. Despite these observations, the high correlation factor demonstrates the ability of the computational method to predict well the stability trend of CYD•1-10 complexes.The analysis of the inclusion complex geometries was carried out by two distance descriptors related with the preferential recognition of each compound into the large and small CYD rims.

The statistical analysis carried out onto the ensembles obtained after full minimization reveals some cases ([CYD•2], [CYD•5], [CYD•10]
) with relevant Boltzmann populations not related to their number of occurred docking configurations.This effect was due to the relative probability weight of each host-guest configuration that can dramatically bias the Boltzmann analysis.Anionic compounds 1 and 2 clearly prefer to bind into the large cavity orienting the acidic moieties in the internal side of the primary oxygen atoms.Several hydrogen bonds with such hydroxyl stabilize both guest compounds in the most stable configurations (Figure 4).Compounds 3 and 4, only differing in one N substitution, show an interaction pattern similar in terms of rim recognition with a week preference for the large one.The most stable CYD inclusion complexes reveal two hydrogen bonds between CYD primary hydroxyls and the secondary amine moiety of 3 and no intermolecular hydrogen bonds, i.e. pure hydrophobic interaction, for compound 4. Compound 5 shows a preference for the small rim inclusion with a best configuration binding mode comparable to the closely related compound 3.The compound 6 brings a positive charge onto the sidechain tertiary nitrogen.The inclusion preference is statistically toward the large rim.In the global minimum energy configuration the ammonium group interacts with a secondary hydroxyl of the large rim with one hydrogen bond interaction, on the other hand the net positive charge of such a moiety induces a local conformational change of another secondary hydroxyl due electrostatic repulsion (Figure 4).Compound 7 interacts with CYD mainly in the small rim.Its extended conformation interacts with a consistent hydrogen network (Figure 4).Similar considerations can be done for the complex with the other alcohol 8, interacting mainly in the small rim accepting in its most stable configuration two hydrogen bonds from primary CYD hydroxyls.The carboxylic moiety of compound 9 has been modeled in neutral form according to the acid pH conditions.The lack of the negative net charge probably is the cause of the reduced free energy of complexation with respect to anionic form of compound 2. The preferred rim recognition is for the small one.Four hydrogen bonds stabilize the global minimum configuration (Figure 4).The ammonium cation 10 interacts in one conformation only with CYD.Also in this case the preferred rim recognition is for the small one.The global minimum configuration shows the positively charged nitrogen exposed in the small rim cavity toward the solvent.It is worth noting the electrostatic effect of such charged atom onto all primary hydroxyl hydrogen atoms forced to loss the intramolecular hydrogen bond network of the small ring (Figure 4).The hydrophobic effect not explicitly considered in our simulations is evaluated by the GB/SA implicit model of solvation mimicking the presence of water.As shown in the energy minimum configurations (figure 4) most of hydrophobic moieties of compounds 1-10 fit into the cyclodextrin cavity minimizing their exposition to the solvent.Other effects, such as the volume and the dipole moments of the guests respectively taken into account by Van der Waals and electrostatic terms of the force field, are responsible of the different recognition of similar compounds, as in the case of acids 2 and 9.

Conclusions
A set of ten molecularly diverse compounds was considered for a docking study with the inhouse MOLINE program. 1 After the conformational search of the isolated guest molecules the docking with three "open" conformations of β-cyclodextrin was performed following the "quasiflexible" approach of the methodology.The ensemble obtained after the full minimization of the complex configurations was thermodynamically evaluated with computational estimation of state functions related to the recognition process.Theoretical and experimental free energies of complexation were compared observing an high level of correlation with r 2 > 0.9.The statistical analysis of the inclusion geometries and the visual inspection of global minimum configurations revealed the trend of anionic compounds to interact with CYD into the large cavity and, in most cases, intermolecular hydrogen bond stabilization with hydroxyl moieties of the small ring.

Figure 2 .
Figure 2. Chemical structures of the ten ligands sorted by the experimentally determined free energy of complexation with β-cyclodextrin.Anionic and cationic forms are related to the free energy experimental conditions.4

Figure 3 .
Figure 3.Correlation plot between experimental and theoretical free energy of complexation in kcal/mol computed at 300°K.The linear correlation is reported as solid line (r 2 = 0.938).

in Honor of Prof. Vincenzo Tortorella ARKIVOC 2004 (v) 107-117 ISSN 1424-6376 Page 111
the results of the conformational search are reported.

Table 1 .
Summary of Monte Carlo search results a no conformational search was carried out because no rotatable bonds exist.

Table 2 .
Thermodynamics of the cyclodextrin 1-10 complexes computed with MOLINE at 300 K a Number of occurred docking configurations found within 12.5 kcal/mol above the global minimum after full energy minimization.

Table 3 .
Statistical details of the CYD complexes