Molecular modeling of the interaction of novel hydroxy-and aminobisphosphonates with hydroxyapatite

Bisphosphonates (BPs) are a class of drugs widely used in the treatment of several metabolic bone disorders associated with increased bone resorption. Although BPs can directly inhibit the cellular activity of osteoclasts, their ability to adsorb on bone mineral is also an important factor in determining their potency and duration of action. In this study, we performed a molecular mechanics analysis of the molecular structures of 18 novel hydroxy-and aminobisphosphonates and examined their interactions with hydroxyapatite. From the calculated interaction energies, general rules were extracted relating structural characteristics of BPs and their affinities to the mineral. The results obtained were shown to be in agreement with in vitro and in vivo studies performed for some of the studied BPs.


Introduction
Bisphosphonates (BPs) are a class of drugs widely used in the treatment of several metabolic bone disorders associated with increased bone resorption, including osteoporosis, Paget's disease and bone metastases.All osteopathies characterised by an excess (absolute or relative) of osteoclastic activity, comprising about 90% of all bone disorders, are now treated with bisphosphonates. 1 In addition, BPs are also useful as novel ligands for well-defined radioactive metal complexes that can be used both for bone scanning imaging and for therapeutic applications. 24][5][6] These compounds are characterized by two phosphonate groups linked to a central carbon atom (Scheme 1).They are analogues of pyrophosphate, in which the oxygen atom of the central P-O-P structure has been replaced by carbon, resulting in a P-C-P group.BPs have chemical and physical properties similar to pyrophosphate but, due to their P-C-P backbone, they are considerably more resistant to heat and enzymatic hydrolysis.Scheme 1.General structure of a bisphosphonate.
Besides the phosphonate groups, which are responsible by the general strong affinity and selectivity of BPs to the bone mineral, these compounds also possess two side groups (R 1 and R 2 in Scheme 1) which determine their structure-function specific profile.For example, the side groups fine tune binding to bone and determine in large extent the antiresorptive potency of the bisphosphonates. 7,8Thus, even though all bisphosphonates share many pharmacological features, subtle structural differences give rise to BPs that are biochemically distinct, particularly in the way in which they bind to bone mineral and on their effects on bone resorption, making it necessary to evaluate each compound individually.
Molecular modeling is an adequate and practical general methodology for evaluation of the relationship between molecular structure and biological function.Molecular modeling allows the determination of molecular geometries (molecular shape, bond lengths, bond angles), physicochemical properties, binding energies and the simulation of intermolecular interactions.Among the methods available to compute molecular structures, molecular mechanics (MM) is able to predict structures for large molecular systems with high accuracy and little computational expense.In molecular mechanics the molecules are considered as a collection of masses (corresponding to the atoms) held together by elastic forces (corresponding to the bonds).These forces can be represented by several terms describing the various structural aspects of a molecule.Giving these terms functional forms we obtain a mathematic expression (Equation 1) that, when combined with certain parameters (obtained from experimental data or high level ab initio theoretical calculations), defines a force field and allows calculation of the energy of a molecular system as a function of its molecular geometry.E = ∑Er + ∑Eθ + ∑EΦ + ∑Evdw + ∑Eel (+ ∑Eothers) The steric or potential energy (E), which represents the difference between the energies of the real molecule and the hypothetical molecule where all the structural parameters assume simultaneously their natural values, can be described as the sum of the contributions resulting from the valence interactions (bond stretching, Er, bond bending, Eθ, and torsions, EΦ) and from the non-bonded interactions (e.g., electrostatic interactions, Eel, van der Waals forces, Evdw). 9,10n this study, with the intent of evaluating specific binding abilities of BPs to HAP and find simple general rules relating structural characteristics of BPs to their affinities to the mineral that could be used to search for pharmacologically more promising BPs, we performed an extensive molecular modeling analysis of the structures of 18 novel hydroxy-and aminobisphosphonates 11 (Scheme 2; see also Table S1 in the Supplementary Information for the nomenclature of the studied compounds) and of their interaction with hydroxyapatite [Ca10(PO4)6(OH)2].The Universal Force Field (UFF) from Rappé and co-workers was selected for the calculations.12 The results obtained were compared with structural data of BPs obtained by X-ray crystallography and in vitro and in vivo hydroxyapatite affinity studies performed for some of the investigated BPs.[13][14][15] Scheme 2

Results and Discussion
Before simulating the interactions of the bisphosphonates with hydroxyapatite, a detailed study of the preferred conformations assumed by the isolated BPs molecules was undertaken.This study provided the necessary basic knowledge on the molecular structure of the BPs, thus allowing for the subsequent identification of the main structural changes suffered by the BPs due to the interaction with hydroxyapatite.For BP1 and BP2 (see Scheme 2) an extended systematic conformational search was performed.The conformationally relevant torsional angles were incrementally rotated (see "Computational Methods" section for details) to generate the initial structures to be submitted to the geometry optimization procedures, allowing, in this way, a more thorough search of the conformational space.For the conformationally more flexible bisphosphonates BP3 to BP18, we choose to use a probabilistic methodsimulated annealing.The simulated annealing is a molecular dynamics technique applied very often to conformational searches that basically consists of the 'heating' of a system to elevated temperatures, so that different energy states are sampled, followed by the slow 'cooling' of the system.In theory, if the temperature is sufficiently high and the cooling process slow, the global energy minimum can be found.Table 1 shows the number of conformers obtained for each BP within an energy limit of 5 kcal mol -1 and the energy differences between the lowest and second lowest energy conformers.These energy differences suggest that for many of the compounds, the most stable conformer is considerably more stable than the other forms and will then assume a particular practical importance, thus justifying its use as reference structure in the simulation of the interactions between the BPs and hydroxyapatite (see Figure S1 of the Supplementary Information for the lowest energy conformations obtained for all studied bisphosphonates).
In order to test the performance of the selected force field, the geometries of BP2, BP6 and BP12 calculated by molecular modeling were compared with those obtained experimentally by X-ray crystallography 13 (Figure 1; see also Tables S2-S4 of the Supplementary Information for numerical comparison of bond length values).Even though the structures obtained experimentally correspond to the crystalline statewhere crystal packing and the dielectric effect may have a significant impact on the molecular structureswhereas the calculated structures refer to the isolated molecules, the theoretical and experimental data show a fairly good agreement, validating the use of the selected theoretical method.In particular, calculated bond lengths and bond angles fit well the experimental values, the most significant differences being observed for bonds that involve hydrogen atoms, where the experimental uncertainty is also larger due to the non-isotropic electron distribution (Tables S2-S4).It is also worth mentioning that the calculated and experimental conformations around the central carbon atom in the case of the more flexible BP6 and, specially, BP12 molecules were found to be somewhat different in the crystal and for the isolated molecule situation (Figure 1), as could be expected taking into account the considerable number of low energy conformers they possess and the low energies associated with internal rotations about the C-R 2 and C-P bonds.When interaction of BPs with hydroxyapatite is considered, the calculations clearly reveal that BPs conformationally relax onto the crystal surface during the energy minimization process, adopting different conformations from those corresponding to the conformational ground state of the isolated molecule.In most of the cases, the BPs positioned themselves in a surface trigonal cavity defined by the HAP calcium ions.As described in previous reports, 16,17 the most important interactions between the ligands and the surface are electrostatic in nature and occurred between the phosphonate groups and nitrogen atoms of the side groups of the BPs and the calcium ions of hydroxyapatite.The oxygen atoms of the phosphonate groups and the nitrogen atoms have relatively high partial negative charges (-0.2 to -0.8e) and interact more strongly with the positively charged calcium atoms (Figure 2).S2 of the Supplementary Information for side views of the lowest energy complexes formed between others bisphosphonates and HAP).Cgrey; Nblue; Ored; Pviolet; Hwhite; Ca 2+green.In the top view representation, hydroxyapatite is showed with a light grey coloration, except for the calcium atoms involved in the coordination.
The calculated interaction energies (Table 2) indicate that, among all the studied bisphosphonates, BP5 binds strongly to hydroxyapatite, with an interaction energy almost two fold larger than those of the other bisphosphonates.This larger interaction can be justified by the fact that BP5 possesses two additional phosphonate groups when compared with the other BPs.As shown in Figure 2, these supplementary phosphonate groups are involved in additional interactions with the calcium ions of HAP.On the contrary, BP12 showed a weaker interaction with HAP (see Table 2).Compared to the other BPs, BP12 has ethyl groups connected to the oxygen atoms of the phosphonate groups instead of hydrogen atoms.This substitution then precludes effective phosphonate/calcium ions interaction, reducing the binding ability of BP12.The interaction energy of BP18 is also low, when compared with those of the other bisphosphonates.BP18 has the same number of phosphonate groups as all other bisphosphonates studied (except BP5) but does not possess the typical geminal bisphosphonate P-C-P structure that is then shown to be of fundamental importance in determining the binding capability of the compounds through calcium coordination.a Average of 50 independent minimizations.Reported errors are standard deviations expressed as a percentage of the interaction energy.
In general, all the remaining bisphosphonates were predicted to have similar interaction energies with hydroxyapatite.Nevertheless, the calculated values for BP1, BP2, BP3 and BP4 are in agreement with results of studies of in vitro and in vivo binding affinities to hydroxyapatite performed for these bisphosphonates marked with 153 Sm. 14,15 The interaction energies calculated for BP2 and BP4 are very similar (difference of 0.1%) and larger than those obtained for BP1 and BP3.In the in vitro studies, 14 the BP2 and BP4 exhibit similar binding percentages to different concentrations of hydroxyapatite and they also bind more strongly than BP1 and BP3.In addition, in vivo biodistribution studies performed on mice for BP1 and BP2 14,15 showed that the two species exhibit different pharmacokinetic parameters, with the 153 Sm-BP2 complex presenting initial higher levels of bone uptake than the 153 Sm-BP1 complex, which concentration is increasing during the 24 h period studied.
The present study demonstrates the good performance of the UFF force field and of the consequent molecular modeling methodology in predicting relative affinities of a series of novel BPs to the bone mineral hydroxyapatite.General rules relating structural characteristics of BPs and their affinities to HAP and information on the nature of the dominant interaction forces could be established: i) the interactions are essentially of electrostatic nature, involving the phosphonate and the nitrogen atom of the BP side groups and the calcium ions of HAP; ii) the geminal P-C-P moiety is essential to ensure a favourable coordination of the BPs to the calcium; iii) increasing the number of P-C-P fragments in the BP molecule makes the interaction energy proportionally larger, so that production of BPs with a larger number of geminal bisphosphonate groups might lead to substances presenting better pharmacological properties; iv) among the series of the studied BPs, subtle differences in the structure of the substituents bound to the bisphosphonate group were predicted to affect only slightly the binding affinity to HAP, though the calculated interaction energies were found to be in consonance with the available in vitro and in vivo studies even in this case.
In the present study, the preferred conformers of the investigated BP molecules were also determined and the relative energies of other less stable conformers estimated.It was shown that for many of the compounds the most stable conformer is considerably more stable than the other forms and will then assume a particular practical importance.Nevertheless, most of the compounds were also shown to have several conformers with energies not higher than 5 kcal mol -1 in relation to the most stable form, testifying their appreciable conformational flexibility.This feature may also be considered relevant in terms of the ability of BPs to bind efficiently HAP, since it makes it possible to attain better geometrical arrangements of the BP molecule on the mineral surface.

Computational methods
Before simulating the interactions with hydroxyapatite, the preferred conformations assumed by the isolated BP molecules were studied.For BP1 and BP2 a systematic conformational search on the potential energy surface of the BPs was undertaken, in which all conformationally relevant torsional angles were rotated by increments of 60 degrees, generating for each BP a total of 1296 different starting conformations.For bisphosphonates BP3 to BP18 a probabilistic method was used (simulated annealing).In this latter case, ten judiciously selected conformations of each BP were submitted to simulating annealing, generating a total of 2000 different conformations.All conformations were energy-minimized using a smart algorithm (cascade of the steepest descent, adjusted basis set Newton-Raphson and quasi-Newton methods) with convergence criteria of gradient less than 0.5 kcal Å -1 mol -1 or 5000 iterations.For each BP, the lowest energy conformer was selected and used as starting geometry in the studies of interaction of the BP with hydroxyapatite.The interaction energies were calculated by initially minimizing the energy of the selected bisphosphonate far (20 Å) from hydroxyapatite.The ligand was then brought closer to the surface (ca. 10 Å) and energy minimized again.50 different orientations of the BP placed at 10 Å from the (001) Miller face of the HAP crystal, in vacuo, were subjected to minimization.The difference in energies between the two systems is the interaction energy.For BP2, BP6 and BP12 the geometry of the low energy conformation was also compared with those obtained experimentally by X-ray crystallography. 13The crystalline structure of hydroxyapatite reported by Kay and Young 18 was used in the simulations and was considered fixed during the energy minimization process.The Universal Force Field (UFF) of Rappé et al. 12 was used in all calculations.Charges were assigned by the charge equilibration method (Qeq). 19except for the calcium ions of HAP for which a full ionic charge of +2 was used.All simulations were made using the Cerius-2 system of programs, running on Silicon Graphics O2 workstations. 20

Figure 2 .
Figure 2. Side view (a) and top view (b) of the lowest energy complex formed between BP5 and hydroxyapatite (see FigureS2of the Supplementary Information for side views of the lowest energy complexes formed between others bisphosphonates and HAP).Cgrey; Nblue; Ored; Pviolet; Hwhite; Ca 2+green.In the top view representation, hydroxyapatite is showed with a light grey coloration, except for the calcium atoms involved in the coordination.

Table 1 .
Number of conformers found for each bisphosphonate, within an energy limit of 5 kcal mol -1 relative to the lowest energy conformer and energy difference between the lowest and the second lowest energy conformers (∆E)

Table 2 .
Interaction energies between the bisphosphonate ligands and hydroxyapatite, as determined by molecular mechanics using the UFF force field