Ab initio investigation of the adsorption of zoledronic acid molecule on hydroxyapatite (001) surface: an atomistic insight of bone protection
aa r X i v : . [ c ond - m a t . m t r l - s c i ] S e p Ab initio investigation of the adsorption of zoledronic acid molecule onhydroxyapatite (001) surface: an atomistic insight of bone protection
Mun-Hyok Ri a , Chol-Jun Yu ∗ ,a , Yong-Man Jang b , Song-Un Kim b a Department of Computational Materials Design (CMD), Faculty of Materials Science, and b Department of Organic Chemistry, Faculty of Chemistry, Kim Il Sung University, Ryongnam-Dong, Taesong-District, Pyongyang, DPR Korea
Abstract
We report a computational study of the adsorption of zoledronic acid molecule on hydroxyapatite (001) surface within ab initio density functional theory. The systematic study has been performed, from hydroxyapatite bulk and surface,and zoledronic acid molecule to the adsorption of the molecule on the surface. The optimized bond lengths and bondangles were obtained and analyzed, giving an evidence of structural similarity between subjects under study. Theformation energies of hydroxyapatite (001) surfaces with two kinds of terminations were computed as about 1.2 and1.5 J / m with detailed atomistic structural information. We determined the adsorption energies of zoledronic acidmolecule on the surfaces, which are -260 kJ / mol at 0.25 ML and -400 kJ / mol at 0.5 ML. An atomistic insight ofstrong binding a ffi nity of zoledronic acid to the hydroxyapatite surface was given and discussed. Key words:
Zoledronic acid, Hydroxyapatite,
Ab initio method, Surface adsorption, Bone
1. Introduction
Bisphosphonates (BPs) are widely used as the mostpowerful drug for protecting bone and treating bone dis-ease such as osteoporosis and other metabolic bone dis-orders [1, 2, 3, 4, 5]. From the 1960s onward, a greatnumber of extensive studies on BPs with respect to thedevelopment and clinical treatment have been carriedout, but the research field on BPs is still very active,being continued to evolve towards a better understand-ing of treatment mechanisms and a finding of more po-tent BP on the basis of it. In this context, zoledronicacid (or zoledronate, ZOD), whose chemical nameis [1-hydroxy-2-(1H-imidazol-1) ethylidene] bisphos-phonic acid or 2-(imidazol-1-yl)-1-hydroxy-ethane-1,1-diphosphonic acid, was developed as a third-generationBP and was approved for the oral treatment of bone dis-ease in 2012 in DPR Korea after the success of its syn-thesis in our own way by two of the authors (Yong-ManJang and Song-Un Kim).BPs are characterized by two phosphonate groups,central carbon atom in between, and two side groups R1and R2, while human bone is composed of complex ar-ray of hydroxyapatite (HAP, Ca (PO ) OH) crystallites ∗ Corresponding author
Email address: [email protected] (Chol-Jun Yu) with nano sizes ranging from 30 to 200 nm embeddedwithin the collagen matrix [6, 7, 8]. It is well estab-lished that the function of BPs to inhibit bone resorp-tion is originated with their ability of binding stronglyto bone mineral, that is, HAP crystal, and sti ff resis-tance to chemical and enzymatic hydrolysis [4]. In moredetailed explanation, the P-C-P backbone of BPs madefrom two phosphonate groups and carbon atom has ahigh binding a ffi nity for the compositions of HAP –tetrahedral PO groups, OH groups and Ca ions – withan additional contribution from the R1 side group (inmost cases including ZOD, it is -OH) [4, 9]. Moreover,the P-C-P group of BPs is considerably more resistant tothe dissolution of HAP crystal than the P-O-P group ofpyrophosphate, of which BPs are stable structural ana-logues [10, 11, 12, 13, 14].Varying the another side group R2 of BPs can resultin di ff erences in antiresorptive potency with several or-ders of magnitude. It was observed that more potent BPsposses a primary, secondary or tertiary nitrogen atom inthe R2 side chain [15, 16]. At present, the most potentantiresorptive BPs include a heterocyclic R2 side chaincontaining a nitrogen atom like risedronate and zole-dronate [17]. The structure-activity relationship anal-ysis showed that the strong binding a ffi nity of zole-dronate for HAP is related to its 3D shape and atomicorientation, indicating an important role of 3D shape of Preprint submitted to J. Mater. Sci. October 16, 2018 itrogen-containing BP and the orientation of its nitro-gen in binding a ffi nity for HAP [4].Atomistic modeling and simulations are a powerfultool to describe the structural characteristics of BPs andbone mineral HAP and the interaction between them atatomic scale, as proved in materials science and molec-ular science through a vast number of applications. Forinstance, molecular dynamics (MD) simulations basedon the well-constructed classical force field have pro-vided the structural information and energetics of bonemineral, HAP crystal and its surfaces [18, 19, 20, 21].Bhowmik et al. [18] obtained the structural parame-ters of monoclinic HAP crystal and its surface ener-getics using the consistent valence force field (CVFF),presenting a valuable description of the interaction be-tween polyacrylic acid and HAP. The surface energet-ics of HAP crystalline surfaces using ab initio densityfunctional theory (DFT) calculations within the gener-alized gradient approximation (GGA) for the exchange-correlation functional has been studied by Zhu andWu [19], testing the e ff ects of slab thickness, vacuumwidth between slabs and surface relaxation on surfaceenergy. Barrios [21] has investigated the interactionbetween collagen protein and HAP surface by usinga combination of computational techniques, DFT andclassical MD methods. Duarte et al. [22] performedmolecular mechanics simulations for molecular struc-tures of 18 novel hydroxyl- and amino-bisphosphonatesto examine the interaction between BPs and hydroxya-patite and to extract relating structural characteristics ofBPs and their a ffi nities for the mineral, which are inagreement with in vitro and in vivo studies for some ofthe studied BPs. To the best our knowledge, however,investigation on surface adsorption of zoledronate onHAP surface with its detailed atomistic structure basedon quantum mechanics is still scarce, in spite of suchextensive theoretical studies of BPs and HAP surface,and we believe that ab initio study on this phenomenonshould definitely contribute to a better understandingof the interaction of zoledronate with bone mineral atatomic and electronic scale.In this paper, we carry out systematic ab initio DFTcalculations for zoledronate molecule, hydroxyapatitebulk and surface, and surface adsorption of zoledronatemolecule on hydroxyapatite surface. Our special focusis placed on the adsorption of zoledronate molecule onhydroxyapatite surface, providing the adsorption energyand atomistic structures of surface complexes, and aninsight how charge transferring is occurred in the eventof adsorption. This is aimed to get a reliable insightfor bone protection e ff ect of zoledronate at electronicscale. In the following, we first describe a computa- tional method in Sec. 2, present the results for structuralparameters and electronic properties of hydroxyapatitebulk and surface, zoledronate molecule, and for bindingof zoledronate to hydroxyapatite surface in Sec. 3, andfinally give our conclusions in Sec. 4.
2. Computational Method
For the DFT calculations in this work, we haveemployed SIESTA code[23] which solves numericallyKohn-Sham equation within DFT [24, 25] using a lo-calized numerical basis set, namely pseudo atomic or-bitals (PAO), and pseudopotentials for describing theinteraction between ionic core (nucleus plus core elec-trons) and valence electrons. The BLYP GGA func-tional (the Becke exchange functional [26] in conjunc-tion with the Lee-Yang-Parr correlation functional [27])was used for exchange–correlation interaction betweenelectrons. For all the atoms, Troullier-Martins [28]type norm-conserving pseudopotentials were generatedwithin local density approximation (LDA) [29], andchecked carefully. The basis sets used in this work werethe DZP type (double ζ plus polarization). The meshsize of grid, which is controlled by energy cuto ff to setthe wavelength of the shortest plane wave representedon the grid, has taken a value of 200 Ry. Non-fixedatoms were allowed to relax until the forces convergeless than 0.02 eV / Å.We first determine the crystal lattice parameters ofbulk HAP by performing structural relaxation withatomic coordinates optimization using the conjugategradient scheme and appropriate Monkhorst-Pack k-points. Through a comparison of the lattice constantswith the experimental values, we have a confidence ofthe above mentioned selection for calculation parame-ters. Structural optimization for isolated ZOD moleculeis performed as well, and this for several conformationsto select the most stable one, of which total energy isthe lowest among studied conformations.We then select the interesting surface index as (001)on the basis of experimental evidence that this faceprovides binding sites for many ionic species [30,31]. In order to simulate the surface within a three-dimensional simulation code, we use two-dimensionalperiodic slabs, which have a thickness of atomic layersof 2 crystal unit cells plus 30 Å vacuum layer, allowingthe atoms in one half of a crystal unit cell at each sideof the slab to relax. These structural parameters in thesupercell can give well-converged result: increasing thethickness of slab to 3 and 4 crystal unit cells and thevacuum width up to 50 Å result in variations of surfaceenergy smaller than ± / m . To check the stability2f the surface, the surface formation energies are calcu-lated approximately as the di ff erence between the totalenergies of the slab and the corresponding bulk crystal: γ = E slab − N slab N bulk E bulk ! / A , (1)where N slab and N bulk are the numbers of atoms in thesurface slab and in the bulk unit cell, A is the slab sur-face area, and E slab and E bulk are the total energies ofthe slab and bulk structures, respectively [32]. The signof surface formation energy is a test of surface stabil-ity: positive (negative) means energy should be pro-vided (released) in order to create a surface. We notethat there might be several possible cutting planes for aparticular Miller index while guaranteeing neutrality ofthe surface charge perpendicular to the surface[33].Final step of the work is to simulate the adsorptionof ZOD molecule on the relevant relaxed HAP (001)surfaces. In order to have the rough, initial surface ad-sorption structure, we use classical molecular mechan-ics (MM) method. In this MM simulation, GeneralUtility Lattice Program (GULP) [34] is utilized. Fordescribing the interaction between atoms, we use theDreiding force field, where the potential energy is de-scribed as the sum of the contributions resulting fromthe bonded interactions (bond stretching, bond bend-ing, and torsions) and from the non-bonded interactions(e.g., electrostatic interaction, and van der Waals inter-action) [35]. After making a guess for the initial stateof the adsorption complex, we also carry out atomic re-laxation for the surface complex – ZOD molecule andsurface atoms in the slab – to obtain a stable final state.Then, the adsorption energy can be calculated as fol-lows: E ads = N mol [ E mol-slab − ( E slab + N mol E mol )] , (2)where N mol is the number of adsorbed molecules, and E mol-slab and E mol are the total energies of the adsorbedsurface and of the isolated molecule, respectively. Theadsorption energy can be either negative or positive:negative (positive) means energy should be released(provided) during the adsorption, indicating that the ad-sorption is (not) spontaneous exothermic (endothermic)process.
3. Results and discussion
The crystal structure of hydroxyapatite is hexagonalwith space group P / m , which contains a formula unit Ca (PO ) (OH) . In fact, there must be four hydroxylgroups in the unit cell with the P / m space group,each oxygen atom of hydroxyl group with 1 / P / m to P . How-ever, the modified model has a net electric polarizationcontrary to the experiment which shows zero polariza-tion, because all the hydroxyl groups in the model areoriented in the same direction. To mimic the real struc-ture where exists disorder in the relative orientation ofthe parallel OH channels, so that electric polarizationsare compensated each other, we have created a super-cell containing two unit cells in the [100] direction, andassigning opposite orientations to the two OH channelsin the supercell. The crystal structure with antiparal-lel hydroxyl groups in a double unit cell compared tothe original hexagonal structure is monoclinic with P space group. We have used this structure in this work,thus allowing a direct comparison of the simulation re-sults with experimentally determined surfaces. Figure 1: (Color online) Top view (a) and side view (b) of supercell ofbulk hydroxyapatite crystal with monoclinic P space group, whichis doubled the unit cell in [100] direction so as to compensate theelectric polarization by assigning opposite orientations to the two OHchannels in hexagon surrounded by triangle of Ca1 ions. (a) top viewand (b) side view. (Ca: Green, P: purple, O: red, and H: white) The full optimization of the structure, allowing thecell shape and volume, and ionic positions to relax, wasperformed. We have checked the convergence with re-spect to the special k-points; increasing the k-points setfrom (1 × ×
3) to (2 × ×
6) leads to the change be-3ween the total energies as about 5 meV, thus indicat-ing the safe use of (1 × ×
3) set. Figure 1 depicts thefully optimized supercell of bulk HAP in this work. Thecalculated structural parameters and chemical bondingproperties of bulk HAP are given in Table 1. The lat-tice constants ( a = b = c = a = b = c = a is the half of latticeconstant of the supercell. Table 1: Structural parameters of bulk hydroxyapatite crystal struc-ture with monoclinic P space group, compared with experimentalresults. Lattice parametersThis work Experiment a a , b , c (Å) 9.348, 9.352, 6.955 9.430, 9.430, 6.891 α, β, γ ( ◦ ) 90, 90, 119.86 90, 90, 120Average of bond lengths (Å)P − O 1.578 1.540Ca1 − O 2.423 2.478Ca1 − O H − O 2.553 2.556Average of bond angles ( ◦ )O − P − O 109.44 109.45O − Ca1 − O 98.62O − Ca2 − O 99.40 a Ref. [36]
As shown in Figure 1, there are two di ff erent Ca sites(denoted as Ca1 and Ca2), four oxygen sites (O1, O2,and O3 of phosphate tetrahedron, and O H of hydroxylgroup), one P site, and one H site. The Ca1 atom is coor-dinated to seven oxygen ions, which are six oxygens ofdi ff erent five PO groups and one oxygen of OH group,while the Ca2 atom has a ninefold coordination withoxygen ions situated in six di ff erent PO tetrahedra.The three Ca1 atoms also form a triangle at the sameplane normal to the c axis, whose center is occupiedby an OH group, and the two triangles form a hexag-onal screw configuration. The calculations yield aver-age bond lengths as 1.578 Å for P-O bond, 2.423 Å forCa1 − O, 2.339 Å for Ca1 − O H , and 2.553 Å for Ca2 − O,compared to the corresponding experimental values of1.54, 2.478, 2.354, and 2.556 Å, respectively. The sim-ilar agreement is obtained in the bond angles, as shownin Table 1. These comparisons confirm that our com-putational parameters are reasonably satisfactory for agood assessment of surface calculations.
In this work, we have selected the Miller indices ofHAP surface as (001), since there exists experimentalevidence that the (001) surface is the most stable amongseveral surfaces with low indices and provides bindingsites for many ionic species [40]. As we use the super-cell doubled in the [100] direction as a unit cell for bulkHAP, the (001) surface unit cell is an oblique (2 ×
1) sur-face cell and there might be two kinds of terminations;(1) 2Ca2 − (6PO · · − · · OH) − − (3PO · · OH),denoted as Type (II). They are repeated in the [001] di-rection, assigned to building layer, and belong to theTasker (II)-type surface [33] with no electric dipole mo-ment perpendicular to the surface.
Table 2: Surface energies (J / m ) according to vacuum thickness andbuilding layers of hydroxyapatite (001) surfaces with (2 ×
1) surfaceunit cell. The values in parenthesis represent the number of atoms.
Type (I); 2Ca2 − (6PO · · − a Vacuum (Å) 20 1.460 1.20830 1.461 1.20540 1.448 1.19750 1.449 1.205Layer 4 (176) 1.461 1.2056 (264) 1.474 1.2188 (352) 1.496Type (II); (3PO · · OH) − − (3PO · · OH)Layer 4 (176) 2.084 1.514 1.016 (264) 2.057 1.5068 (352) 2.078 a Ref. [21]
The HAP (001) surface has been modelled usingthree-dimensional periodic supercell (slab) with twoequivalent surfaces at bottom and top side. To ensure nointeraction between bottom surface of the above imageslab and top surface of the present slab across the vac-uum region, the vacuum region must be wide enough.And the atomic layer, which is consisted of surface layerallowed to relax and crystal layer fixed at its crystallineposition, should be thick enough so that the two sur-faces of each slab do not interact through the crystallayer. To check the convergence according to the vac-uum region, we tested 20, 30, 40, and 50 Å thickness,and confirmed that there is no distinct change between20 and 50 Å thick vacuums (Table 2). Therefore, thevacuum region of 30 Å thickness will be used in thefollowing calculations. The thickness of the slab is usu-ally expressed in terms of a number of building layers,4 igure 2: (Color online) Top view (upper panel) and side view (lower panel) of optimized structure of hydroxyapatite (001) surface with(2 ×
1) surface unit cell, vacuum thickness of 30 Å, 4 building layers, and termination 2Ca2 − (6PO · · − · · OH) − − (3PO · · OH) (b). The top and bottom layers are allowed to relax, while the layers between dashed lines are fixed attheir crystalline positions. (Ca: green, P: purple, O: red, H: white) where one layer contains 44 atoms. The four layers (twobulk unit cells, 176 atoms), six layers (three unit cells,264 atoms), and eight layers (four unit cells, 352 atoms)were tested with the vacuum width of 30 Å, and thechange of surface energies between 4 layers slab and8 layers slab is only 0.05 J / m . Thus we will use a slabwith 4 building layers in the study of surface adsorption.As listed in Table 2, the surface energy of Type (II)surface is slightly higher than that of Type (I) surface,indicating the Type (I) terminated surface is more fa-vorable than the Type (II) terminated surface. We seethat the surface relaxation in the Type (I) surface is notreally much compared to the Type (II) surface, since thedi ff erence between unrelaxed and relaxed surface ener-gies in the former case is not remarkable contrary to thelatter case. Since there is no data available in the litera-ture for Type (I), we only compared the calculated Type(II) surface energy with the previous SIESTA work [21].Figure 2 shows the fully relaxed atomistic structureof the HAP (001) surfaces modelled by a slab with vac-uum thickness of 30 Å and four building layers. It is ob- served that the coordination number (CN) of Ca1 atomchanges from 7 in bulk to 6, and CN of Ca2 from 9 to 6in the Type (I) termination, while in the Type (II) termi-nation the CN of Ca1 atom changes from 7 to 6, and Ca2atom from 9 to 5, and 6. The coordination environmentsaround the surface oxygen atoms (O1, O2, O3, and O H )are also changed from the their bulk environment, whileP atoms are still fully surrounded by four oxygen atomslike in bulk.Table 3 shows the bond lengths of cation-oxygen andbond angles of oxygen-cation-oxygen at the top surfacelayer. In the case of Type (II) surface, we consider thetwo kinds of Ca1 and three kinds of Ca2 surface atoms.In both cases, the averages (1.583, 1.585 Å) of bondlengths of P − O are a bit expanded compared with thebulk (1.578 Å), and the averages of bond angles getsmaller than in the bulk. The bond lengths of Ca2 − Oare clearly contracted at the surface; 2.434 in Type (I),and 2.325, 2.380, and 2.414 Å in Type (II), which are allsmaller than the bulk value 2.553 Å. From these obser-vations, it can be concluded that the undercoordinated5 able 3: Bond lengths of cation-oxygen and bond angles of oxygen-cation-oxygen at hydroxyapatite (001) surface with coordinationnumbers (CN) of cations.
Type (I); 2Ca2 − (6PO · · − ◦ )Length (Å) Range Average CNP 1.583 103.47 ∼ ∼ ∼ · · OH) − − (3PO · · OH)P 1.585 102.07 ∼ ∼ ∼ ∼ ∼ ∼ surface Ca atoms and O atoms can be favorable adsorp-tion sites. To simulate an isolated molecule, we have used a cu-bic supercell with lattice constants of a = b = c =
50 Å,which has been proved to be enough to prevent the ar-tificial interaction between adjacent molecules. Theremight be four di ff erent conformations, as presented inFigure 3. We can make a distinction between the confor-mations according to the directions of two pairs of OHgroups attached to the two P atoms; (1) ZOD Bout for thecase where the directions of both OH pairs are outwardagainst the nitrogen heterocyclic ring, (2) ZOD
Bin forthe case of inward directions of both pairs, (3) ZOD
Nin for the case of inward direction of one pair in the sameside of nitrogen atom, and (4) ZOD
Nout for the case ofoutward direction of one pair in the N side. After thetotal energy minimizations to get the optimized struc-ture of molecule, we have compared the total energiesof four conformations. The calculations tell us the moststable structure is just the third case, ZOD
Nin , which willbe used in the following adsorption study.The optimized bond lengths and bond angles relatedwith P atoms in ZOD
Nin conformation are listed in Ta-ble 4. The averages of P − O bond lengths are 1.597Å in P1 side and 1.586 Å in P2 side, which are simi-lar to those of bulk HAP (1.578 Å) and of HAP (001)surface (1.583 Å). The averages of O − P − O bond an-gles (112.11 and 111.42 ◦ ) are also close to those of bulkHAP (109.44 ◦ ) and HAP surface (109.39 ◦ ). It is foundthat those values of P2 side are slightly closer to thebulk and surface values than P1 side, so that P2 side are PHO HOOHO HOON N HOP
Figure 3: (Color online) Optimized structures of zoledronic acid withfour di ff erent conformations, (1) ZOD Bout , (2) ZOD
Bin , (3) ZOD
Nin ,and (4) ZOD
Nout . (C: grey, P: purple, N: blue, O: red, H: white) more favorable to binding with HAP due to the closerstructural similarity.Figure 4 shows the computed isodensity surfacesfor the highest occupied molecular orbitals (HOMO)and lowest unoccupied molecular orbitals (LUMO) ofZOD (conformation ZOD
Nin ). While the HOMO andHOMO-1 are localized on the heterocyclic ring con-taining nitrogen atom, LUMO and LUMO + able 4: Optimized bond lengths and angles related with phosphorusatoms in zoledronic acid conformation, ZOD Nin . Ave. means average. P relatedBond length (Å) Bond angle ( ◦ )P1 = O1 1.495 O1 = P1 − O2 117.26P1 − O2 1.654 O1 = P1 − O3 115.10P1 − O3 1.641 O1 = P1 − C1 115.22P1 − C1 1.902 O2 − P1 − O3 103.97O2 − P1 − C1 102.96O3 − P1 − C1 100.12Ave. O − P1 − O 112.11Ave. O − P1 − C1 106.10P relatedP2 = O4 1.495 O4 = P2 − O5 120.01P2 − O5 1.622 O4 = P2 − O6 117.35P2 − O6 1.641 O4 = P2 − C1 99.93P2 − C1 1.934 O5 − P2 − O6 96.91O5 − P2 − C1 112.73O6 − P2 − C1 110.40Ave. O − P2 − O 111.42Ave. O − P2 − C1 107.69P1 − C1 − P2 98.65
Figure 4: (Color online) Molecular orbitals of zoledronic acid confor-mation, ZOD
Nin . heterocyclic ring containing nitrogen atom could be adonor. To begin with adsorption, we have utilized GULP toobtain rough structure of ZOD absorbed on HAP (001)surface, where Dreiding forcefield was adopted. Simu-lated annealing was performed, increasing the tempera-ture from 300 K to 10000 K and subsequently decreas-ing back to 300 K with the temperature interval of 50 K.Using the obtained rough structure as the starting struc-ture, we then performed atomic relaxation with SIESTAcode to get the final optimized structure of zoledronicacid adsorbed HAP (001) surface.In this work we have tried to simulate the adsorptionof ZOD molecule at 0.5 ML (one molecule on (2 × ×
2) surface cell). As mentioned above, there aretwo possible terminations in the HAP (001) surface, andtherefore four kinds of simulation tasks were carriedout. We have used the slab model with 4 building layersand vacuum thickness of 30 Å, which guarantee to givea reliable adsorption energy. To make a systematic er-ror small, the pristine surface energies for (2 ×
2) surfaceslabs were also calculated. In Table 5, the adsorption
Table 5: Calculated adsorption energies (kJ / mol) of zoledronic acidon hydroxyapatite (001) surface. Coverage0.5 ML (2 ×
1) 0.25 ML (2 × − − − − energies are listed. All the adsorption energies are neg-ative, indicating that all the adsorptions are exothermicreactions. We see that the adsorption on Type (II) sur-face is slightly more favorable than on Type (I) at both0.5 ML and 0.25 ML, although the Type (I) surface for-mation from the bulk needs more energy than the Type(II) surface.Figure 5 shows the optimized atomistic structure ofZOD adsorption complexes on HAP (001) surfaces at0.25 ML coverage. It is observed that in the Type(I) surface the hydrogen atom of ZOD’s phosphonategroup moved to oxygen atom of HAP surface’s phos-phate group, making formation of additional OH groupon the surface, and thus indicating that the adsorptionis kind of chemisorption caused by proton exchange.There are several hydrogen bonds between ZOD’s phos-phonate OH group and oxygen atoms of the surface’sphosphate groups, and vice versa. It is also found thatCa1 and Ca2 atoms with 6 coordination number (CN)7 igure 5: (Color online) Top view (upper panel) and side view (lower panel) of optimized atomistic structure of zoledronic acid adsorptioncomplexes on hydroxyapatite (001) surfaces with Type (I) (a) and Type (II) (b) terminations at 0.25 ML coverage. (Ca: green, P: purple, O: red, N:blue, H: white) make bonds with oxygen and nitrogen atoms of ZOD,respectively. Such bonding was expected based on thesurface relaxation analysis and HOMO-LUMO distri-bution of isolated ZOD molecule. In the case of Type(II) surface, there are also hydrogen bonds between thesurface and ZOD, but the move of hydrogen atom is notobserved, which indicates that the adsorption is kind ofphysisorption. Nevertheless, the magnitude of adsorp-tion energy on Type (II) surface is larger than that onType (I) surface.To make it clear the charge transfer occurred duringthe adsorption, we show the electron density di ff erence, ρ surf + mol − ( ρ surf + ρ mol ), plots in Figure 6, where (a) is forthe isosurface figure with the value of ± / Å , and(b) and (c) are the contour plots on the planes around hy-drogen bond and Ca cation, respectively. These figuresclearly give the evidence of charge transferring at theevent of adsorption, from hydrogen of phosphate groupof ZOD molecule to oxygen of PH group of HAP sur-face making formation of hydrogen bond, and from Caof HAP surface to oxygens of phosphonate group form-ing the coordinate bond. Therefore the adsorption ofZOD on the HAP (100) surface is surely chemisorption. We calculated the density of states (DOS) of elec-trons, total and partial DOS shown in Figure 7 and atomprojected and partial DOS in Figure 8. In Figure 7 wesee the energy spectrum of isolated ZOD molecule (a),the DOS of ZOD molecule that was adsorbed on HAPsurface (b), the DOS of HAP (100) surface (c), and thewhole complex system comprised of the adsorbed ZODmolecule and HAP (100) surface. Some hybridizationbetween ZOD and HAP surface electrons is shown inthe figures. Which atoms cause the charge transferringat the adsorption? To make an answer to this question,we carefully investigate the atom projected partial DOS.The 1s peak of hydrogen of isolated ZOD moleculeplaced over 0 eV (red line in Figure 7 (a) panel) is re-markably weakened after the adorption (blue line in (a)panel), indicating the loss of electrons from hydrogenand the gain by oxygen, and thus the formation of hy-drogen bond between hydrogen atom of ZOD and oxy-gen of HAP surface. We also see the hybridizations be-tween hydrogen 1s state and oxygen 2p state in (a) and(b) panels, and between oxygen 2s state and calcium 3pstate in (c) and (d) panels.8 igure 6: (Color online) Electronic charge density di ff erence of zoledronic acid adsorption complex on hydroxyapatite (001) surface. (a) isosurfacefigure with the value of ± / Å (yellow for + and blue for - value), (b) contour figure on the plane containing hydrogen bond between O ofphosphonate group of HAP surface and hydrogen of phosphate group, and (c) contour figure around Ca cation. (Ca: green, P: purple, O: red, N:blue, H: white)
4. Conclusion
In conclusion, we have attempted to make a model-ing of adsorption of zoledronic acid on hydroxyapatite(001) surface to get an atomistic insight of bone protec-tion. The systematic study has been performed, fromhydroxyapatite bulk and surface, and zoledronic acid toadsorption of the molecule on the (001) surface. Wehave carried out the structural optimizations and atomicrelaxations of the bulk, molecule, surface, and adsorp-tion complexes on the surfaces, and obtained the struc-tural informations and energetics.Using the three-dimensional periodic supercellmodel, we determined the stable conformation of themolecule and calculated the molecular orbitals. It wasconcluded that two phosphonate groups can play a roleof electron acceptor in the chemical reaction, while theheterocyclic ring containing nitrogen atom can be aelectron donor. After verifying the validity of compu-tational parameters of SIESTA work through the bulk hydroxyapatite calculation, surface modeling and relax-ations were performed, and surface formation energieswere calculated for two kinds of (001) surface termina-tions, which are about 1.2 and 1.5 J / m . Subsequently,the adsorption of zoledronic acid on the relaxed sur-face was studied, obtaining the surface binding structureand calculating the adsorption energies. We found thatthe surface Ca atoms and oxygen atoms of phosphategroups can form surface bond including hydrogen bondand coordinate bond with nitrogen, hydrogen, and oxy-gen atoms of zoledronic acid. The calculated adsorp-tion energies are about -260 kJ / mol at 0.25 ML cover-age and -400 kJ / mol at 0.5 ML coverage, indicating thestrong binding a ffi nity of zoledronic acid to hydroxyap-atite surface.We have made an interpretation of such strong bind-ing a ffi nity through the analysis of atomistic structures,electron density di ff erence, the density of states andHirshfeld charges of atoms relevant to surface binding.The results showed that bond lengths and bond angles9 spdTot D e n s it y o f S t a t e s ( e V - ) -20 -10 0 10Energy (eV)020004000 ZODZOD on HAPHAPZOD + HAP(a)(b)(c)(d)
Figure 7: (Color online) Total and partial density of states in iso-lated zoledronic acid (a), adsorbed zoledronic acid (b), hydroxyap-atite (100) surface (c) and zoledronic acid adsorption complex on thesurface (d). of phosphonate group of zoledronic acid are similar tothose of phosphate group of hydroxyapatite bulk andsurface, indicating the strong binding a ffi nity relatedto the structural similarity. It was also found that thecharge transferring is occurred mainly on the side ofphosphonate group in one kind surface, while in anothersurface it is occurred on both imidazol ring contain-ing nitrogen atom and phosphonate group of zoledronicacid. Acknowledgments
The simulations have been carried out on the HPBlade System c7000 (HP BL460c) that is owned andmanaged by Faculty of Materials Science, Kim Il SungUniversity. This work was supported from the Com-mette of Education (grant number 02-2014), DPR Ko-rea. -20 -10 0 10Energy (eV)0100200
Ca 4sCa 3p D e n s it y o f s t a t e s ( e V - ) O 2sO 2p
ZOD OHZOD OH on HAPHAP O on HAP
Isolated ZODZOD on HAP (a)(b)(c)(d)
H 1sO 2s O 2p
Figure 8: (Color online) Atom projected partial density of states; (a)for 1s electrons of hydrogen in di ff erent situations, (b) for 2s and 2pelectrons of oxygen in di ff erent situations, (c) for 2s and 2p states ofoxygens surrounding Ca cation in ZOD adsorbed HAP surface, and(d) for 4s and 3p of calcium forming the coordinate bond on HAPsurface. Notes
The authors declare that they have no conflict of in-terest.
References [1] F. H. Ebetino, A.-M. L. Hogan, S. Sun, M. K. Tsoumpra,X. Duan, J. T. Triftt, A. A. Kwaasi, J. E. Dunford, B. L. Barnett,U. Oppermann, M. W. Lundy, A. Boyde, B. A. Kashemirov,C. E. McKenna, R. G. G. Russell, The relationship between thechemistry and biological activity of the bisphosphonates, Bone49 (2011) 20–33.[2] R. Bartl, B. Frisch, E. von Tresckow, C. Bartl, Bisphosphonatesin Medical Practice, Springer-Verlag, Berlin, Germany, 2007.[3] F. P. Coxon, K. Thompson, M. J. Rogers, Recent advances in un-derstanding the mechanism of action of bisphosphonates, Curr.Opin. Pharmacol. 6 (2006) 307–312.
4] R. G. G. Russell, N. B. Watts, F. H. Ebetino, M. J. Rogers,Mechanisms of action of bisphosphonates: similarities and dif-ferences and their potential influence on clinical e ffi cacy, Osteo-poros Int. 19 (2008) 733–759.[5] R. G. G. Russell, Bisphosphonates: The first 40 years, Bone 49(2011) 2–19.[6] T. Hassenkam, G. E. Fantner, J. A. Cutroni, J. C. Weaver, D. E.Morse, P. K. Hansma, High-resolution AFM imaging of intactand fractured trabecular bone, Bone 35 (2004) 4.[7] P. Alberius-Henning, E. Adolfsson, J. Grins, Triclinic oxy-hydroxyapatite, J. Mater. Sci. 36 (2001) 663.[8] J. C. Elliot, Structure and Chemistry of the Apatites and OtherCalcium Orthophosphates, Elsevier Science, Amsterdam, 1994.[9] R. G. G. Russell, R. C. Muhlbauer, S. Bisaz, D. A. Williams,H. Fleisch, The influence of pyrophosphate, condensed phos-phates, phosphonates and other phosphate compounds on thedissolution of hydroxyapatite in vitro and on bone resorptioninduced by parathyroid hormone in tissue culture and in thy-roparathyroidectomised rats, Calcif Tissue Res. 6 (1970) 183.[10] E. G. Kontecka, J. Jezierska, M. Lecouvey, Y. Leroux, H. Ko-zlowski, Bisphosphonate chelating agents. Coordination abil-ity of 1-phenyl-1-hydroxymethylene bisphosphonate towardsCu(2 + ) ions, J. Inorg. Biochem. 89 (2002) 13.[11] H. Fleich, Bisphosphonates: mechanisms of action, Endocr.Rev. 19 (1998) 80.[12] M. J. Rogers, S. Gordon, H. L. Benford, F. P. Coxon, S. P. Luck-man, J. Monkkonen, J. C. Frith, Cellular and molecular mech-anisms of action of bisphosphonates, Cancer Suppl. 8 (2000)2961.[13] A. J. Roelofs, K. Thompson, S. Gordon, M. J. Rogers, Molec-ular mechanisms of action of bisphosphonates: current status,Clin. Cancer Res. 12 (2006) 6222.[14] J. E. Dunford, K. Thompson, F. P. Coxon, S. P. Luckman, F. M.Hahn, C. D. Poulter, F. H. Ebetino, M. J. Rogers, Structure-Activity Relationships for Inhibition of Farnesyl DiphosphateSynthase in Vitro and Inhibition of Bone Resorption in Vivo byNitrogen-Containing Bisphosphonates, J. Pharm. Exp. Therap.296 (2001) 235–242.[15] G. H. Nancollas, R. Tang, R. J. Phipps, Z. Henneman, S. Gulde,W. Wu, A. Mangood, R. G. G. Russell, F. H. Ebetino, Novelinsights into actions of bisphosphonates on bone: Di ff erences ininteractions with hydroxyapatite, Bone 38 (2006) 617.[16] M. A. Lawson, Z. Xia, B. L. Barnett, J. T. Tri ffi tt, R. J. Phipps,J. E. Dunford, R. M. Locklin, F. H. Ebetino, R. G. G. Russell,Di ff erences Between Bisphosphonates in Binding A ffi nities forHydroxyapatite, J. Biomed. Mater. Res. B 92B (2010) 149.[17] J. R. Green, K. Muller, K. A. Jaeggi, Preclinical pharmacologyof CGP 42 446, a new, potent, heterocyclic bisphosphonate com-pound, J. Bone Miner. Res. 9 (1994) 745–751.[18] R. Bhowmik, K. S. Katti, D. Katti., Molecular dynamics simu-lation of hydroxyapatite-polyacrylic acid interfaces, Polymer 48(2007) 664.[19] W. Zhu, P. Wu, Surface energetics of hydroxyapatite: a DFTstudy, Chem. Phys. Lett. 396 (2004) 38.[20] K. Matsunaga, A. Kuwabara, First-principles study of vacancyformation in hydroxyapatite, Phys. Rev. B 75 (2007) 014102.[21] N. A. Barrios, A Computational Investigation of the Interactionof the Collagen Molecule with Hydroxyapatite, Ph.D. thesis,Department of Chemistry, University College London, 2010.[22] L. F. Duarte, F. C. Teixeira, R. Fausto, Molecular modeling ofthe interaction of novel hydroxyl- and aminobisphosphonateswith hydroxyapatite, ARKIVOC (v) (2010) 117.[23] J. M. Soler, E. Artacho, J. D. Gale, A. Garc´ıa, J. Junquera,P. Ordej´on, D. S´anchez-Portal, The Siesta method for ab ini-tio order-N materials simulation, J. Phys:Condens. Matter 14 (2002) 2745.[24] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Phys.Rev. 136 (1964) B864–B871.[25] W. Kohn, L. J. Sham, Self-consistent equations including ex-change and correlation e ff ects, Phys. Rev. 140 (1965) A1133–A1138.[26] A. D. Becke, Density-functional exchange-energy approxima-tion with correct asymptotic behavior, Phys. Rev. A 38 (1988)3098–3100.[27] C. Lee, W. Yang, R. G. Parr, Development of the Colle-Salvetticorrelation-energy formula into a functional of the electron den-sity, Phys. Rev. B 37 (1988) 786.[28] N. Troullier, J. L. Martins, E ffi cient pseudopotentials for plane-wave calculations, Phys. Rev. B 43 (1991) 1993–2006.[29] J. P. Perdew, A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys.Rev. B 23 (1981) 5048.[30] N. Almora-Barrios, K. Austen, N. de Leeuw, A Density Func-tional Theory study of the binding of glycine, proline and hy-droxyproline to the hydroxyapatite (0001) and (0110) surfaces,Langmuir 25 (2009) 5018–5025.[31] A. Wierzbicki, H. S. Cheung, Molecular modeling of inhibitionof hydroxyapatite by phosphocitrate, J. Mol. Struc-Theochem.529 (2000) 73–82.[32] C.-J. Yu, J. Kundin, S. Cottenier, H. Emmerich, Ab initio mod-eling of glass corrosion: hydroxylation and chemisorption ofoxalic acid at diopside and åkermanite surfaces, Acta Mater. 57(2009) 5303–5313.[33] P. W. Tasker, The stability of ionic crystal surfaces, J. Phys. C12 (1979) 4977–4984.[34] J. D. Gale, A. L. Rohl, The general utility lattice program(GULP), Mol. Sim. 29 (2003) 291–341.[35] A. Leach, Molecular Modelling: Principles and Applications,Prentice Hall, 2001.[36] J. Y. Kim, R. R. Fenton, B. A. Hunter, B. J. Kennedy, Powderdi ff raction studies of synthetic calcium and lead apatites, Aust.J. Chem. 53 (2000) 679.[37] A. S. Posner, A. Perlo ff , A. F. Diorio, Refinement of the hydrox-yapatite structure, Acta Cryst. 11 (1958) 308.[38] N. H. de Leeuw, Density functional theory calculations of localordering of hydroxy groups and fluoride ions in hydroxyapatite,Phys. Chem. Chem. Phys. 4 (2002) 3865–3871.[39] D. E. Ellis, J. Terra, O. Warschkow, M. Jiang, G. B. Gonzalez,J. S. Okasinski, M. J. Bedzyk, A. M. Rossi, J. G. Eon, A the-oretical and experimental study of lead substitution in calciumhydroxyapatite, Phys. Chem. Chem. Phys. 8 (2006) 967–976.[40] N. H. de Leeuw, A computer modelling study of the uptakeand segregation of fluoride ions at the hydrated hydroxyap-atite (0001) surface: introducing a Ca (PO4)6(OH) potentialmodel, Phys. Chem. Chem. Phys. 6 (2004) 1860–1866.potentialmodel, Phys. Chem. Chem. Phys. 6 (2004) 1860–1866.