Surface Structure of Hydroxyapatite from Simulated Annealing

Article pubs.acs.org/Langmuir

Surface Structure of Hydroxyapatite from Simulated Annealing Molecular Dynamics Simulations Hong Wu,† Dingguo Xu,*,‡ Mingli Yang,*,† and Xingdong Zhang§ †

Institute of Atomic and Molecular Physics, Key Laboratory of High Energy Density Physics and Technology of Ministry of Education, Sichuan University, Chengdu, Sichuan 610065, China ‡ College of Chemistry, and §National Engineering Research Center for Biomaterials, Sichuan University, Chengdu, Sichuan 610064, China S Supporting Information *

ABSTRACT: The surface structure of hydroxyapatite (HAP) is crucial for its bioactivity. Using a molecular dynamics simulated annealing method, we studied the structure and its variation with annealing temperature of the HAP (100) surface. In contrast to the commonly used HAP surface model, which is sliced from HAP crystal and then relaxed at 0 K with firstprinciples or force-field calculations, a new surface structure with gradual changes from ordered inside to disordered on the surface was revealed. The disordering is dependent on the annealing temperature, Tmax. When Tmax increases up to the melting point, which was usually adopted in experiments, the disordering increases, as reflected by its radial distribution functions, structural factors, and atomic coordination numbers. The disordering of annealed structures does not show significant changes when Tmax is above the melting point. The thickness of disordered layers is about 10 Å. The surface energy of the annealed structures at high temperature is significantly less than that of the crystal structure relaxed at room temperature. A three-layer model of interior, middle, and surface was then proposed to describe the surface structure of HAP. The interior layer retains the atomic configurations in crystal. The middle layer has its atoms moved and its groups rotated about their original locations. In the surface layer, the atomic arrangements are totally different from those in crystal. In particular for the hydroxyl groups, they move outward and cover the Ca2+ ions, leaving holes occupied by the phosphate groups. Our study suggested a new model with disordered surface structures for studying the interaction of HAP-based biomaterials with other molecules.

1. INTRODUCTION Biomaterials are used in medical devices that interact with biological systems. A suitable biomaterial, particularly its surface, must possess specific biological activity, for example, biocompatibility, functionality, and chemical stability, for usage in vivo.1,2 Biomineralization is such a process that makes the material surface recognizable to the various cells and molecular components in the matrix.3,4 The importance of the surface structures of biomaterials has long been addressed.5,6 Synthetic hydroxyapatite (HAP) is one of the most often used bioactive materials because of its good performances in biocompatibility, nontoxicity, and osteoinductivity.7 Many studies have demonstrated that the growth or inhibition of bone and teeth is related to the biomolecule−HAP interaction at their interfaces.8,9 Moreover, osteoinductivity was observed at the HAP surfaces with disordered micro- and nanopores.10,11 However, explorations to the HAP surface structure are hindered by experimental limits. Although many advanced techniques, such as nuclear magnetic resonance (NMR), photoelectronic energy spectra (PES), and atomic force microscopy (AFM), etc., have been employed in their characterization,12−14 knowledge about the disordered surface structures of bioactive HAP is still scarce. Jäger et al.15 developed an NMR strategy to explore the surface structure of © XXXX American Chemical Society

HAP nanocrystals and found both the composition and the structure in the surface layer were different from those in crystal. Different Ca/P ratios in HAP surface and crystal were observed by Lu et al.16 with X-ray photoelectron spectroscopy (XPS) and time-of-flight secondary ion mass spectrometry (TOF-SIMS). Battistoni et al.17 studied the surface chemical composition and morphology of the HAP coatings by means of XPS, scanning Auger microscopy (SAM), and scanning electron microscopy (SEM), revealing that the cleanness, homogeneity, and roughness of the coatings varied with the environment. Using a combined NMR spectroscopy and firstprinciples calculations approach, Chappell et al.18 found the measured 31P chemical shifts of surface atoms are significantly different from those in the crystal. França et al.19 characterized the surface structure of nano-HAP and its biphasic mixture with β-tricalcium phosphate by means of XPS, X-ray diffraction (XRD), SEM, Fourier-transform infrared (FTIR), TOF-SIMS, and laser granulometry, and found the Ca/P and Ca/O ratios in the outer layer differ significantly from those of bulk. Received: December 22, 2015 Revised: March 25, 2016

A

DOI: 10.1021/acs.langmuir.5b04667 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir On the other hand, numerous computational studies20−22 have been conducted to investigate the HAP structure and its interaction with various biomolecules. Zhu et al.23 have performed periodic density functional theory (DFT) calculations to investigate the energetics of the hydroxyapatite surfaces by using a slab model. At the first-principles level, Lou et al.,24 Corno et al.,25 Almora et al.,26 Filgueiras et al.,27 and Rimola et al.28 have studied the interaction of HAP surface with water, lysine, citric, glycine, and proline, revealing that chemical bonding, electrostatic forces, hydrogen bonding, and van der Waals forces contribute mainly to the interaction between the adsorbates and the HAP surfaces. Meanwhile, the structures and properties of the adsorbates are affected by the HA surface. Molecular dynamics simulations at the force-field level have been conducted to explore the adsorption behaviors of polypeptides and proteins on the HAP surface. Almora et al.,29 Dong et al.,30 and Shen et al.31 have studied the interaction of collagen and bone morphogenetic proteins (BMP-2 and BMP-7), which has great influence on the growth and inhibition of bone, with the HAP surfaces. The selectivity of interacting sites and the preference of protein orientations were addressed. Electrostatic attraction was found dominant in the interaction energy. The disordering of HAP surface structure has also been noted in a couple of computational simulations.32−34 Lee et al.33 studied the structure relaxation of HAP surface with interatomic potentials and found three kinds of distortion in the relaxation process, polarization of layers close to the surface, distortion of the channels in which the hydroxide ions sit, and distortion of the tetrahedral phosphate ions. Using classical molecular dynamics simulations at 310 K, Leeuw34 studied the structure relaxation of hydrated HAP surfaces, which were created from the bulk crystal. The surfaces showed a regular structure at the interface due to the water coordination effect, while some hydroxyl groups near the interface diffused into the solvent. However, all of these studies treated the HAP structure as a perfect crystal using the common slab model in which the HAP surface was created by simply slicing the crystal and intersecting a vacuum layer between adjacent sliced pieces. The created surface was then relaxed using an energy minimization approach at the first-principles or force-field level. The atoms or groups on the surface then have small movements or rotations around their original locations to reach a stationary structure during the relaxation, while the chemical composition and atomic arrangement do not show distinct changes from those in bulk. As mentioned above, such a simplified surface model for HAP surface is in contradiction with the experimental evidence.15−19 It is therefore necessary to reconstruct the HAP surface structure in an approach closer to the experimental approach. Experimentally, bioactive HAP is usually produced from HAP powders with a sintering and cooling process. The HAP powders are heated to a relatively high temperature (lower than or around its melting point) and sintered to have certain shape, strength, and flexibility. The processing at high temperature may make atoms and groups escape from their original locations and move into a more expansive zone, which is not reflected in a usual structure relaxation in simulations. In this work, we mimicked the preparation of HAP ceramics using a molecular dynamics simulated annealing approach. The simulated annealing algorithm was first proposed by Kirkpatrick et al.,35 and has been widely used in the prediction of

molecules, clusters, and solids with unknown structures. The HAP surface structures were constructed using the slab model and then relaxed in a designed annealing approach. The variation of surface structures with annealing temperature and their differences from the HAP crystal structure were highlighted. We found the obtained annealed surfaces have different chemical compositions, atomic arrangements, and stability from the unannealed ones, which will definitely affect the surface reactivity toward other molecules. Our study proposes a reliable approach for the construction of HAP surface, which is useful and critical for a proper simulation of interactions between HAP and biomolecules.

2. METHODS Hydroxyapatite has a hexagonal crystal structure of P63/m space group.36 Its measured lattice parameters are a = b = 9.424 Å, c = 6.879 Å, α = β = 90°, and γ = 120°.37 A slab model was used to construct the initial surface of HAP. The (100) surface was sliced and separated with a vacuum layer of 30 Å, which is wide enough to avoid interaction from adjacent structures under three-dimensional periodic boundary conditions. Hauptmann’s38 potential energy function, which was developed specifically for HAP, was used in our simulation. In this model, the force field consists of an intra-atomic and an interatomic term:

U = Uintra + Uinter

(1)

where

Uij ,inter(rij) =

⎡R + R − r ⎤ CC 1 qiqj i j ij ⎥− i j + w(ρi + ρj ) exp⎢ ⎢⎣ ρi + ρj ⎥⎦ 4πε0 rij rij6 (2)

and

Uintra =

∑ bonds

+

∑ angles

kr (rij − r0)2 + 2

∑ angles

kθ (θij − θ0)2 2

kUB (rik − r0)2 2

(3) 39,40

potential In fact, the widely used Born−Mayer−Huggins (BMH) was defined in Uinter. ρi, Ri, and Ci are BMH parameters, and w = 1.1552 × 10−19 J Å−1.41 ε0, r, and q are vacuum dielectric constant, particle distance, and site charge, respectively. The partial atomic charges and nonbonded force field parameters, which were used for the simulation, are given in Table S1, and the intra-atomic force field parameters are given in Table S2.42 All of these parameters have been commonly used for the HAP-containing systems in many studies.31,43,44 Because surface atoms usually have different statuses from inner atoms, the force field parameters suitable for inner atoms are not necessarily suitable for surface atoms. Therefore, one needs to look into the status change of atoms/groups when their locations change from inner to surface. Defining different atomic types for the surface and inner atoms in the force field is a possible way to meet the requirements. However, the alternation of atomic types for atoms/ groups whose locations change between on surface and inside during simulation processes often leads to convergence difficulty. Using the same parameters for both surface and inner atoms is an approximation adopted in many studies.45,46 We then validated the usage of BMH parameters for HAP surface structures. On the basis of the unannealed and annealed structures, we designed two surface models to compare the changes in atoms/groups status by using density functional theory (DFT) calculations. The computational details and results are given in the Supporting Information. In the unannealed model, Ca ions are exposed directly to the vacuum, and their averaged net charge change is quite different from that of inner Ca ions. In the annealed model, however, all surface B


Article

Langmuir

Figure 1. Unannealed (a) and annealed (b) structures of bulk HAP (Ca = yellow, O = red, P = purple, H = white). atoms including Ca have similar atomic net charge and bond lengths with inner atoms. The OH− migration and Ca2+ retraction make the surface covered by O atoms of phosphate group and H atoms of hydroxyl group. These atoms do not have dangling bonds and have similar statuses with in bulk. Therefore, it is a reasonable approximation to use the same force field parameters for surface and inner atoms of annealed HAP structures. Various supercells of 6 × 3 × 3 were constructed to investigate the surface thickness and temperature effect. The bottom layer (2 × 3 × 3) of the supercells was fixed to mimic the real inner crystalline HAP, while the other upper layers were relaxed in the simulations. All of the MD simulations were carried out in a canonical (NVT) ensemble with the LAMMPS program.47 A Nose−Hoover thermostat was used in each case to couple the system to a heat reservoir,48 and the Verlet leapfrog algorithm was used to integrate the equations of motion.49 The movement of hydroxyl is crucial for the HAP surface structure. A short time step of 0.1 fs was used for a proper description for the hydroxyl groups, which was found necessary for simulations at high temperatures. The process of simulated annealing is illustrated in Figure S1.50 First, the system was relaxed for 50 ps at ambient temperature of 300 K. It then was heated to a target temperature at a speed of 10 K/ps. Several target temperatures (lower, around, and higher than the HAP melting point, 1923.16 K51) were selected in the simulations. After equilibration for 200 ps at the target temperature, the system was cooled with a stepped cooling pattern. In this pattern, the system was cooled at a speed of 100 K/20 ps, and then equilibrated for 20 ps for every 100 K until it reached 300 K. Finally, the system was equilibrated at 300 K for 100 ps to obtain its annealed surface structure, which was used for analysis of its composition, structure, atom and group distribution, and activity. To describe the surface structure of HAP visually, we have examined two-body structural correlations through the pair distribution function (PDF), static structure factor, and coordination number.52,53 PDF describes the atomic distribution around and can be used as a measure of atomic distribution ordering, which is computed with

g (r ) =

1 Nρ0

N

∑ i=1

S(q) =

Useful information can be obtained by integration around the first peak in the PDF. The number distribution of atom β around atom α, or the average coordination number nαβ as a function of r, is defined as

nαβ (R ) = 4πρβ

1 Nαρβ

Nα

∑ i=1

∫0

R

gαβ (r )r 2 dr

(8)

where R is the radius of the first shell of atom β around atom α. This function represents the assemble-averaged number of atom β surrounding atom α at a distance R. The pair distribution function (PDF) and the static structure factor S(q) for each structural model for annealed HAP and crystal were computed with ISAACS.54

3. RESULTS AND DISCUSSION 3.1. Structure of Bulk HAP. To validate the BMH parameters in a simulated annealing approach, we first simulated the bulk HAP structure and compared it to the measurements.18,38 The above-defined procedures of heating and cooling for the surface simulations were applied to the bulk HAP structure without any constraints. Figure 1a shows the annealed structure of bulk HAP, which is in general consistent with the measured structure (Figure 1b). The locations of Ca2+, OH−, and PO43− have tiny changes, while the orientations of hydroxyl change considerably. We found that about one-half of the hydroxyl groups rotate about 90°. The rearrangement of hydroxyl groups in HAP has long been noticed.24,26,55 The hydroxyl groups in HAP usually have different orientations, induced by their surrounding O atoms of phosphate. Figure 2 compares the total pair distribution functions (PDF) of annealed and perfect HAP structures. The PDF of annealed HAP is basically the same as that of HAP crystal. The peaks/troughs locate at almost the same distances and have similar heights/depths. The reorientation of hydroxyl groups does not lead to a remarkable change in PDF because hydroxyl groups are much fewer in number than the

(4)

where N is the number of the atoms in the simulation cell, ρ0 is the number density, and ni(r) is the number of the particles found in the shell from r to r + Δr from a central atom i. In this work, we also used partial pair distribution function to measure the distribution of atom β around atom α, which is computed with

gαβ (r ) =

(7)

α ,β

ni(r ) 4πr 2Δr

∑ (cαcβ)1/2 Sαβ(q)

niβ (r ) 4πr 2Δr

(5)

To describe the correlation between the atom α and atom β, the partial static structure factor Sαβ(q) is defined as Sαβ(q) = δαβ + 4πρ0 (cαcβ)1/2

∫0

∞

sin qr (gαβ (r ) − 1)r 2 dr qr

(6)

where δαβ is the Kronecker delta, and cα = Nα/N. The total static structure factor S(q) is given by

Figure 2. Total pair distribution functions of unannealed (red line) and annealed (black line) structures of bulk HAP. C


Article

Langmuir

Figure 3. Crystal surface (a) and the annealed surface structures of HAP at 1600, 2000, 2400, and 2800 K (b, c, d, and e).

other ions or groups. In addition, the melting point of bulk HAP is about 2000 K (Figure S2), in agreement with the measurements.51 Therefore, the simulated annealing approach and the BMH parameters used in this work produced reliable results for the HAP systems. 3.2. Effect of Annealing Temperature. Now we focused on the surface structure of annealed HAP. First, the effect of annealing temperature was explored. In a usual simulated annealing approach,52,56 the annealing temperature was usually set to temperatures well above the melting point to ensure that the samples were completely melted. Considering that the temperature in the preparation of HAP biomaterials was usually lower than the melting point, a number of annealing temperatures of 1600, 2000, 2400, and 2800 K were set in our simulations to study the temperature effect on the surface structure. The annealed surface structures are shown in Figure 3. The unannealed crystal HAP surface (Figure 3a) was also given for comparison. When the annealing temperature (1600 K, Figure 3b) is lower than the HAP melting point, some atoms on the surface displace remarkably, while atoms inside only have small displacements. Reorientation of some hydroxyls was observed under this temperature, but their locations are almost unchanged. At 2000 K (Figure 3c), the system reaches its melting point. The atomic displacements increase remarkably, in particular to the atoms near to the surface. The reorientation of hydroxyl and phosphate groups escalates, and some of the hydroxyls move toward the surface and leave vacancies for other ions or groups, which was observed in several experimental studies.12,57 When the annealing temperature reaches 2400 and 2800 K (Figure 3d and e), the obtained structures are similar to each other and to that of 2000 K. Although the atoms under higher annealing temperature move faster and more widely, the influence of annealing temperature becomes less significant when it further increases. Similar surface structures were obtained for temperature above the HAP melting point when the systems were cooled. It is apparent that the surface structure is dependent on the annealing temperature. The surface structure becomes more disordered at higher temperatures, which is depicted in their PDFs in Figure 4. Figure 4a compares the total pair distribution function of the surface structure after annealing at different Tmax values. In the r < 3.00 Å region, five curves are similar,

Figure 4. Total pair distribution functions (a) and the partial pair distribution functions of P−Ca (b) of the surface structure of HAP crystal (black line), and annealed at 1600 K (red), 2000 K (blue), 2400 K (dark cyan), and 2800 K (magenta).

indicating that their interatomic arrangements at short ranges remain the same as in crystal. The tiny peak at r = 0.98 Å corresponds to the hydroxyl groups. The second and the third peaks at r = 1.47 and 2.40 Å are for the P−O and O−O bonds in phosphate groups. Therefore, the structures of hydroxyl and phosphate groups do not have remarkable changes after annealing. In the r > 3.00 Å region, the fluctuation of the curves decays and approaches g(r) ≈ 1, indicating that the surface structures become disordered at large r after annealing. The fluctuations at 1600 K are different from those in crystal and to some extent different from those at 2000 K, while the curves at 2000, 2400, and 2800 K are similar. It is evident that the total g(r) changes when Tmax increases until reaching the melting point, but the changes become small if Tmax increases further. In consequence, the surface structure of hydroxyapatite D


Article

Langmuir

defects caused by unsaturated Ca ions are eliminated by the saturated O and H atoms. Second, the total energy is increased by removal of hydroxyl from inner structure, but it is decreased by the ionic bonding between surface Ca2+ and OH−. As a result, the surface becomes more stable with lower surface energy. Third, the surface activity changes after annealing. The exposed and unsaturated Ca ions have high reactivity with adsorbates. Their coverage by hydroxyl and phosphate groups prevents them from interacting with adsorbates directly. The saturate hydroxyl and phosphate groups have low chemical activity. As a result, the annealed surface has lower reactivity than the unannealed one. 3.4. Three-Layer Surface Model. It has been observed that for small-sized HAP particles the surface thickness at disordered state is about 1 nm.15,60 Because of the lack of neighboring atoms, atoms on the surface are usually unsaturated. The interatomic bonding should be different from that in crystal. Therefore, one or several atomic layers near the surface possess amorphous structure with disordered atomic distribution. The predicted thickness of disordered layers is about 10 Å. Now we take the structure annealed at 2400 K as an example to analyze the structure variation from inner to surface. Similar conclusions can be drawn for the structures annealed at 2000 and 2800 K. The ordered-to-disordered structure variation was observed not only at different temperatures, but also at different depths. From inner to surface, the structure disordering increases. So we divided the structure into three layers, interior, middle, and surface, as shown in Figure 6. The interior layer is fixed during

depends on the annealing temperature but has small changes when the temperature is well above the HAP melting point. The disordering in the surface structure is well described in the P−Ca partial pair distribution function, as show in Figure 4b. The curves for the annealed structures are quite different from that in crystal. The two peaks at about 3.50 Å represent the two kinds of P−Ca distance, but both peaks broaden after annealing. Moreover, all of the peaks at large r change their locations and heights and approach g(r)P−Ca ≈ 1. 3.3. Surface Energy. The variations in surface structure are related to the changes in surface energy Esurf, which is defined as Eslab − E bulk (9) 2A in the slab model. Eslab and Ebulk are potential energies of the slab and bulk HAP, respectively. To calculate Eslab and Ebulk, additional calculations were carried out for the systems with all atoms being relaxed. The supercells used in the calculations are the same as those in the annealing approaches. Eslab was then obtained by the annealing and energy minimization processes. A is the surface area of the slab face. Esurf is thus the averaged surface energy of the two faces. Figure 5 displays the surface energy at different annealing temperatures. The surface energy is 4.73 J/m2 at room Esurf =

Figure 5. Surface energy of HAP annealed at different temperatures.

temperature (300 K), and decreases rapidly until 1200 K, which corresponds to sintering temperature often used in experiments.58,59 Esurf decreases slowly from 1200 to 2000 K, which is around the melting point of HAP. When the temperature increases further, Esurf has very small changes. The variations in surface energy are consistent with the variations in surface structure discussed above. The annealed surface structure becomes stable when Tmax is above the melting point, and is almost unchanged even if using a higher temperature. Furthermore, we compared the surface energy of unannealed and annealed structures at DFT level. The computational models were given in Figure S3. The surface energy of the annealed surface is 0.04 J/m2 lower than that of the unannealed one. The surface reconstruction at high temperature leads to a number of changes in HAP surface structures. The inner hydroxyl groups move outward and cover some of the surface Ca ions, while the other Ca ions move inward and are partially buried by O atoms of phosphate groups. As a result, the surface is occupied by O and H atoms. The surface reconstruction certainly leads to some changes in physical and chemical properties. First, in contrast to surface Ca ions with dangling bonds, the surface O and H atoms are saturate. The surface

Figure 6. Ordered-to-disordered structures of HAP annealed at 2400 K were partitioned into three layers: interior, middle, and surface.

annealing, which is used to simulate the HAP internal structure. The middle layer is about 30 Å in thickness with its top about 10 Å from surface. The atoms or groups in the middle layer basically remain unchanged in location except for small moves or rotations for calcium ions, phosphate, and hydroxyl groups. The surface layer is about 10 Å in thickness with disordered structures different from those in the middle and interior layers. The structure ordering or disordering can be measured with the PDF, static structure factor, and atomic coordination number.52,53 Next, we will compare these three quantities for E


Article

Langmuir

Figure 7. Total PDF (a) and partial PDFs of O−Ca (b), P−Ca (c), and P−O (d) of the three layers of HAP annealed at 2400 K, respectively.

Figure 8. Total structure factor (a) and partial structure factor of O−Ca (b),P−Ca (c), and P−O (d) of the three layers of HAP annealed at 2400 K, respectively.

to g(r) ≈ 1. The relative positions of O−Ca and P−Ca in the surface layer are rearranged and become disordered after annealing. The partial PDFs of P−O (Figure 7d) have very small changes through r = 10 Å, indicating that the phosphate groups keep their structures almost unchanged after annealing, although the groups in the surface layer move and rotate. The structural differences in the three layers can be further highlighted with the static structure factor S(q). Because S(q) reflects the medium-range order of a structure, especially for the first sharp diffraction peak, the structure changes in the middle and surface layers are well described with S(q). Although pair distribution function and structure factor contain similar information, S(q) provides an intuitive description that relates directly with scattering experiments. Figure 8 shows the total structure factor and partial structure factor of O−Ca, P−

the three layers of the structure annealed at 2400 K. Figure 7a− d displays the total PDF and partial PDFs of O−Ca, P−Ca, and P−O, respectively. In Figure 7a, the first-, second-, and thirdnearest-neighbor locations and the heights of three layers are nearly the same, corroborating that their short-range order structures are similar. When r > 3.00 Å, the PDFs of the interior and middle layers oscillate in the same way, reflecting the longrange order structures in these two layers. The PDF of surface layer deviates from those of the former two layers at about r = 3.00 Å. Its fluctuations become decay with r and approach to g(r) ≈ 1. For the partial PDFs of O−Ca (Figure 7b) and P−Ca (Figure 7c), their first-nearest-neighbor locations and heights of the three layers are almost same. Starting from the second peak, the PDFs for the interior and middle layers are similar until r = 10 Å, but the PDFs of the surface layers are different and decay F


Article

Langmuir

Figure 9. Coordination number distribution of P, Ca, and O atoms in the three layers of HAP annealed at 2400 K.

Ca, and P−O of the three layers of the structure annealed at 2400 K, respectively. For the Stotal(q) in Figure 8a, the curves are similar for the interior and the middle layers with their peaks at approximately the same locations and heights, which are much different from those of the surface layer. The S(q) of O−Ca correlations (Figure 8b) for the three layers are different in the heights of peaks. The heights decrease slowly for the layers from interior to middle and to surface, corresponding to the remarkable O−Ca correlation changes in three layers. Because O−Ca correlations are affected by the O atoms both in phosphate and in hydroxyl groups, the long displacement of hydroxyl in the surface layer and the rotation of phosphate groups in the middle layer lead to considerable changes in SO−Ca(q). The changes in SP−Ca(q) (Figure 8c) are moderate in comparison with SO−Ca(q), resulting from that P centers at the phosphate groups and has smaller displacements than its surrounding O atoms in rotation. The curves are almost coincident for the interior and middle layers, but different from that for the surface layer in height and location. The changes in SP−O(q) in Figure 8d come from not only the P−O correlation in phosphate groups, but also from that between P and O in hydroxyl. That is why the SP−O(q) is similar to SO−Ca(q). From the pair distribution functions and structure factors, we can quantitatively analyze the coordination numbers of the P, O, and Ca atoms. The average coordination numbers were obtained by integrating the partial pair distribution function. Figure 9 shows the coordination number distribution of P, Ca, and O atoms in the three layers of structure annealed at 2400 K. The pair of atoms A−B in the figure represents the number of B around A. The coordination numbers of P by O atoms in the phosphate groups are identical in the three layers, as are the coordination numbers of H by O atom in the hydroxyl groups, as presented in Table 1. As shown in Figure 9a, 100% of the P atoms in interior and middle layers have nine Ca atoms around, which is consistent with that in HAP hexagonal crystal. In the surface layer, however, most P atoms have eight Ca atoms around, and the others have six, seven, nine, or 10 Ca atoms around. The average coordination number of P−Ca is 8.24 in the surface, which is less than those in the interior and middle

Table 1. Average Coordination Numbers of the Three Layers of HAP Annealed at 2400 Ka A−B pair

surface

middle

interior

P−Ca Ca−P O−Ca Ca−O P−O H−O

8.24 4.12 2.51 5.51 4.00 1.00

9.00 5.40 2.97 7.72 4.00 1.00

9.00 5.40 3.00 7.80 4.00 1.00

a

The pair of atoms A−B represents the average number of B around A. The P−O pair presents the oxygen atoms in phosphate groups around P, and H−O the oxygen atoms in hydroxyl groups around H.

layers (9.00). The coordination of Ca by P atoms (Figure 9b) is identical in the interior and middle layers. 60% of the Ca atoms have five P atoms around, while 40% of the Ca atoms have six P atoms around. In the surface layer, the Ca atoms have three, four, five, and six P atoms around. The average coordination number is 5.40, 5.40, and 4.12 in the three layers from interior to surface, respectively. 100% of the O atoms in the interior layer are coordinated by three Ca atoms (Figure 9c). In the middle layer, 92.5%, 5.3%, and 2.2% of the O atoms have three, two, and four coordinated Ca atoms, respectively, resulting mainly from the reorientation of hydroxyl groups. In the surface, almost 50% and 50% of the O atoms are coordinated by two and three Ca atoms. The average coordination number of O by Ca in the three layers is 3.00, 2.97, and 2.51, respectively. In Figure 9d, 60% of the Ca atoms coordinate with seven O atoms and 40% with nine O atoms in the interior layer. Because of the reorientation of hydroxyl groups, the coordination number of the Ca atoms becomes seven, eight, or nine in the middle layer. Nevertheless, the Ca−O coordination number in surface is four, five, six, or seven, which have nearly equal fractions. The average coordination number is 7.80, 7.72, and 5.51 in the three layers, respectively, as shown in Table 1. The coordination number distribution of P and Ca atoms reveals that the relative locations of P and Ca atoms change slightly in the middle layer, but change remarkably in the G


Article

Langmuir surface layer. A similar finding was revealed by the coordination number distribution of Ca and O atoms. The coordination numbers of atoms are closely related to their reactivity toward other molecules. The different coordination numbers of surface atoms will definitely affect their reactivity. In our surface model, the migration of hydroxyl from inside to the surface changes the coordination number of Ca atom on the surface, lowering the electrophilicity of surface Ca2+. In most previous studies,28,31,43 the unannealed HAP surface models were conventionally used to simulate their interaction with biomolecules. The main difference between the present three-layer model and conventional models is the application of simulated annealing molecular dynamics. With this approach, we observed the surface reconstruction, like OH− migration, of HAP, which was not addressed in previous conventional models. The OH− migration from inside to surface needs to overcome an energy barrier, which can only be realized at high temperature. Because most HAP-based biomaterials were sintered at high temperature, our model mimics a closer situation with experiments than the conventional ones. The conventional HAP surface models were mostly used to study the adsorption of biomolecules, focusing on how the conformers of biomolecules change on the surface. From this aspect, the studies based on the conventional models provide useful information on the bioactivity of HAP-based materials. However, when one focuses on the materials themselves, our model exhibits some more rationality. It is of particular importance to study the surface structure of biomaterials when tissue inductive materials (TIM) are concerned. Regarded as next-generation biomaterials, TIM has particular surface and near-surface structures that may induce the regeneration of damaged tissues or organs.10,11 The HAP surface needs to be modified to have inductivity. From this aspect, our model is useful for studying the surface and near-surface structures of HAP-based biomaterials.

Atoms and groups in the interior layer have the same distribution within bulk crystal. In the middle layer, atoms and groups move only in their vicinities. However, the atoms and groups are rearranged in the surface layer. The hydroxyl groups move outward and bond with the Ca2+ on the surface. The phosphate groups then occupy the holes left by the hydroxyls. The rearrangements of atoms and groups make the surface layer disordering and distinctly different from those in bulk. Moreover, the coverage of hydroxyls on the Ca2+ ions will definitely affect the reactivity of the HAP surface. Our calculations revealed that the surface structures annealed at high temperatures are more stable than those relaxed at 0 K, suggesting a new reasonable surface model for studying the interaction of HAP-based biomaterials with other molecules.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.5b04667. Additional tables, figures, and references (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail:[email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank the National High Technology Research and Development Program of China (no. 2015AA034202) and the National Natural Science Foundation of China (no. 21373140) for financial support. Part of the calculations were carried out at the Analytical & Testing Center in Sichuan University and the National Supercomputing Center in Shenzhen.

4. CONCLUSION The surface structure of HAP was studied using a molecular dynamics simulated annealing method. The surface structures were constructed by slicing the bulk HAP crystal and then annealed at various temperatures of 1200, 1600, 2000, 2400, and 2800 K. The annealed structures were then analyzed in terms of their atom and group arrangements, energetics, pair distribution functions, static structure factors, and coordination numbers. Our study gives new information about the HAP surface structure that has not been revealed in previous studies. First, the atomic arrangements in the surface layer of HAP are dramatically different from those in bulk crystal. The HAP surface structure was usually obtained by slicing the crystal and then optimized at 0 K. Although the atoms were relaxed to have a stationary surface structure, their locations had very small changes. This approach is contrast to the experiments in which the HAP samples are treated at temperatures near to its melting point. In our approach, the structure relaxation was carried out at various temperatures, and disordered surface structures with lower energy were obtained. Second, the disordering of HAP surface structure is related to annealing temperature. The disordering increases with temperature below the melting point. However, the surface structures annealed at temperatures above the melting point are similar to each other. Third, the surface structure can be partitioned into three layers, interior, middle, and surface, according to their degree of disorder.

■

REFERENCES

(1) Thevenot, P.; Hu, W.; Tang, L. Surface Chemistry Influence Implant Biocompatibility. Curr. Top. Med. Chem. 2008, 8, 270−280. (2) Ducheyne, P.; Qiu, Q. Bioactive Ceramics: The Effect of Surface Reactivity on Bone Formation and Bone Cell Function. Biomaterials 1999, 20, 2287−2303. (3) Gorski, J. P. Acidic Phosphoproteins from Bone Matrix: A Structural Rationalization of Their Role in Biomineralization. Calcif. Tissue Int. 1992, 50, 391−396. (4) Pecheva, E.; Pramatarova, L.; Altankov, G. Hydroxyapatite Grown on a Native Extracellular Matrix: Initial Interactions with Human Fibroblasts. Langmuir 2007, 23, 9386−9392. (5) Wang, K.; Zhou, C.; Hong, Y.; Zhang, X. A Review of Protein Adsorption on Bioceramics. Interface Focus 2012, 2, 259−277. (6) Qin, Z.; Gautieri, A.; Nair, A. K.; Inbar, H.; Buehler, M. J. Thickness of Hydroxyapatite Nanocrystal Controls Mechanical Properties of the Collagen−Hydroxyapatite Interface. Langmuir 2012, 28, 1982−1992. (7) Chen, J.; Chu, B.; Hsiao, B. S. Mineralization of Hydroxyapatite in Electrospun Nanofibrous Poly (L-Lactic Acid) Scaffolds. J. Biomed. Mater. Res., Part A 2006, 79, 307−317. (8) Schachschal, S.; Pich, A.; Adler, H.-J. Aqueous Microgels for the Growth of Hydroxyapatite Nanocrystals. Langmuir 2008, 24, 5129− 5134. (9) Cerroni, L.; Filocamo, R.; Fabbri, M.; Piconi, C.; Caropreso, S.; Condò , S. G. Growth of Osteoblast-like Cells on Porous

H


Article

Langmuir Hydroxyapatite Ceramics: An in Vitro Study. Biomol. Eng. 2002, 19, 119−124. (10) Bose, S.; Roy, M.; Bandyopadhyay, A. Recent Advances in Bone Tissue Engineering Scaffolds. Trends Biotechnol. 2012, 30, 546−554. (11) Wang, H.; Zhi, W.; Lu, X.; Li, X.; Duan, K.; Duan, R.; Mu, Y.; Weng, J. Comparative Studies on Ectopic Bone Formation in Porous Hydroxyapatite Scaffolds with Complementary Pore Structures. Acta Biomater. 2013, 9, 8413−8421. (12) Ben Osman, M.; Diallo-Garcia, S.; Herledan, V.; Brouri, D.; Yoshioka, T.; Kubo, J.; Millot, Y.; Costentin, G. Discrimination of Surface and Bulk Structure of Crystalline Hydroxyapatite Nanoparticles by NMR. J. Phys. Chem. C 2015, 119, 23008−23020. (13) Kobayashi, T.; Nakamura, S.; Yamashita, K. Enhanced Osteobonding by Negative Surface Charges of Electrically Polarized Hydroxyapatite. J. Biomed. Mater. Res. 2001, 57, 477−484. (14) Kirkham, J.; Brookes, S. J.; Shore, R. C.; Wood, S. R.; Smith, D. A.; Zhang, J.; Chen, H.; Robinson, C. Physico-Chemical Properties of Crystal Surfaces in Matrix−mineral Interactions during Mammalian Biomineralisation. Curr. Opin. Colloid Interface Sci. 2002, 7, 124−132. (15) Jäger, C.; Welzel, T.; Meyer-Zaika, W.; Epple, M. A Solid-State NMR Investigation of the Structure of Nanocrystalline Hydroxyapatite. Magn. Reson. Chem. 2006, 44, 573−580. (16) Lu, H. B.; Campbell, C. T.; Graham, D. J.; Ratner, B. D. Surface Characterization of Hydroxyapatite and Related Calcium Phosphates by XPS and TOF-SIMS. Anal. Chem. 2000, 72, 2886−2894. (17) Battistoni, C.; Casaletto, M. P.; Ingo, G. M.; Kaciulis, S.; Mattogno, G.; Pandolfi, L. Surface Characterization of Biocompatible Hydroxyapatite Coatings. Surf. Interface Anal. 2000, 29, 773−781. (18) Chappell, H.; Duer, M.; Groom, N.; Pickard, C.; Bristowe, P. Probing the Surface Structure of Hydroxyapatite Using NMR Spectroscopy and First Principles Calculations. Phys. Chem. Chem. Phys. 2008, 10, 600. (19) França, R.; Samani, T. D.; Bayade, G.; Yahia, L.; Sacher, E. Nanoscale Surface Characterization of Biphasic Calcium Phosphate, with Comparisons to Calcium Hydroxyapatite and β-Tricalcium Phosphate Bioceramics. J. Colloid Interface Sci. 2014, 420, 182−188. (20) Wang, Q.; Wang, M.; Wang, K.; Liu, Y.; Zhang, H.; Lu, X.; Zhang, X. Computer Simulation of Biomolecule−biomaterial Interactions at Surfaces and Interfaces. Biomed. Mater. 2015, 10, 032001. (21) Zhou, H.; Wu, T.; Dong, X.; Wang, Q.; Shen, J. Adsorption Mechanism of BMP-7 on Hydroxyapatite (001) Surfaces. Biochem. Biophys. Res. Commun. 2007, 361, 91−96. (22) del Valle, L. J.; Bertran, O.; Chaves, G.; Revilla-López, G.; Rivas, M.; Casas, M. T.; Casanovas, J.; Turon, P.; Puiggalí, J.; Alemán, C. DNA Adsorbed on Hydroxyapatite Surfaces. J. Mater. Chem. B 2014, 2, 6953−6966. (23) Zhu, W.; Wu, P. Surface Energetics of Hydroxyapatite: A DFT Study. Chem. Phys. Lett. 2004, 396, 38−42. (24) Lou, Z.; Zeng, Q.; Chu, X.; Yang, F.; He, D.; Yang, M.; Xiang, M.; Zhang, X.; Fan, H. First-Principles Study of the Adsorption of Lysine on Hydroxyapatite (100) Surface. Appl. Surf. Sci. 2012, 258, 4911−4916. (25) Corno, M.; Busco, C.; Bolis, V.; Tosoni, S.; Ugliengo, P. Water Adsorption on the Stoichiometric (001) and (010) Surfaces of Hydroxyapatite: A Periodic B3LYP Study. Langmuir 2009, 25, 2188− 2198. (26) Almora-Barrios, N.; Austen, K. F.; de Leeuw, N. H. Density Functional Theory Study of the Binding of Glycine, Proline, and Hydroxyproline to the Hydroxyapatite (0001) and (011̅0) Surfaces. Langmuir 2009, 25, 5018−5025. (27) Filgueiras, M. R. T.; Mkhonto, D.; de Leeuw, N. H. Computer Simulations of the Adsorption of Citric Acid at Hydroxyapatite Surfaces. J. Cryst. Growth 2006, 294, 60−68. (28) Rimola, A.; Corno, M.; Zicovich-Wilson, C. M.; Ugliengo, P. Ab Initio Modeling of Protein/Biomaterial Interactions: Glycine Adsorption at Hydroxyapatite Surfaces. J. Am. Chem. Soc. 2008, 130, 16181− 16183.

(29) Almora-Barrios, N.; De Leeuw, N. H. Molecular Dynamics Simulation of the Early Stages of Nucleation of Hydroxyapatite at a Collagen Template. Cryst. Growth Des. 2012, 12, 756−763. (30) Dong, X.; Wang, Q.; Wu, T.; Pan, H. Understanding Adsorption-Desorption Dynamics of BMP-2 on Hydroxyapatite (001) Surface. Biophys. J. 2007, 93, 750−759. (31) Shen, J.-W.; Wu, T.; Wang, Q.; Pan, H.-H. Molecular Simulation of Protein Adsorption and Desorption on Hydroxyapatite Surfaces. Biomaterials 2008, 29, 513−532. (32) Slepko, A.; Demkov, A. A. First Principles Study of Hydroxyapatite Surface. J. Chem. Phys. 2013, 139, 044714. (33) Lee, W. T.; Dove, M. T.; Salje, E. K. H. Surface Relaxations in Hydroxyapatite. J. Phys.: Condens. Matter 2000, 12, 9829. (34) de Leeuw, N. H. Computer Simulations of Structures and Properties of the Biomaterial Hydroxyapatite. J. Mater. Chem. 2010, 20, 5376. (35) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671−680. (36) Kay, M. I.; Young, R. A.; Posner, A. S. Crystal Structure of Hydroxyapatite. Nature 1964, 204, 1050−1052. (37) Renaudin, G.; Laquerrière, P.; Filinchuk, Y.; Jallot, E.; Nedelec, J. M. Structural Characterization of Sol−gel Derived Sr-Substituted Calcium Phosphates with Anti-Osteoporotic and Anti-Inflammatory Properties. J. Mater. Chem. 2008, 18, 3593. (38) Hauptmann, S.; Dufner, H.; Brickmann, J.; Kast, S. M.; Berry, R. S. Potential Energy Function for Apatites. Phys. Chem. Chem. Phys. 2003, 5, 635−639. (39) Fumi, F. G.; Tosi, M. P. Ionic Sizes and Born Repulsive Parameters in the NaCl-Type Alkali Halides-I. J. Phys. Chem. Solids 1964, 25, 31−43. (40) Tosi, M. P.; Fumi, F. G. Ionic Sizes and Born Repulsive Parameters in the NaCl-Type Alkali Halides-II. J. Phys. Chem. Solids 1964, 25, 45−52. (41) Yuen, P. S.; Murfitt, R. M.; Collin, R. L. Interionic Forces and Ionic Polarization in Alkaline Earth Halide Crystals. J. Chem. Phys. 1974, 61, 2383−2393. (42) Bhowmik, R.; Katti, K. S.; Katti, D. Molecular Dynamics Simulation of Hydroxyapatite−polyacrylic Acid Interfaces. Polymer 2007, 48, 664−674. (43) Pan, H.; Tao, J.; Xu, X.; Tang, R. Adsorption Processes of Gly and Glu Amino Acids on Hydroxyapatite Surfaces at the Atomic Level. Langmuir 2007, 23, 8972−8981. (44) Jiang, W.; Pan, H.; Cai, Y.; Tao, J.; Liu, P.; Xu, X.; Tang, R. Atomic Force Microscopy Reveals Hydroxyapatite−Citrate Interfacial Structure at the Atomic Level. Langmuir 2008, 24, 12446−12451. (45) Brogan, A. P. S.; Sessions, R. B.; Perriman, A. W.; Mann, S. Molecular Dynamics Simulations Reveal a Dielectric-Responsive Coronal Structure in Protein−Polymer Surfactant Hybrid Nanoconstructs. J. Am. Chem. Soc. 2014, 136, 16824−16831. (46) Jiang, S.; Jelfs, K. E.; Holden, D.; Hasell, T.; Chong, S. Y.; Haranczyk, M.; Trewin, A.; Cooper, A. I. Molecular Dynamics Simulations of Gas Selectivity in Amorphous Porous Molecular Solids. J. Am. Chem. Soc. 2013, 135, 17818−17830. (47) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (48) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (49) Verlet, L. Computer “Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (50) Lewis, L. J.; De Vita, A.; Car, R. Structure and Electronic Properties of Amorphous Indium Phosphide from First Principles. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 57, 1594. (51) Lu, Y.-P.; Li, M.-S.; Li, S.-T.; Wang, Z.-G.; Zhu, R.-F. PlasmaSprayed Hydroxyapatite+titania Composite Bond Coat for Hydroxyapatite Coating on Titanium Substrate. Biomaterials 2004, 25, 4393− 4403. (52) Duan, G.; Xu, D.; Zhang, Q.; Zhang, G.; Cagin, T.; Johnson, W.; Goddard, W. Molecular Dynamics Study of the Binary Cu46Zr54 I


Article

Langmuir Metallic Glass Motivated by Experiments: Glass Formation and Atomic-Level Structure. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 224208. (53) Vashishta, P.; Kalia, R. K.; Rino, J. P.; Ebbsjö, I. Interaction Potential for SiO2: A Molecular-Dynamics Study of Structural Correlations. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 12197−12209. (54) Le Roux, S.; Petkov, V. ISAACS-Interactive Structure Analysis of Amorphous and Crystalline Systems. J. Appl. Crystallogr. 2010, 43, 181−185. (55) Hochrein, O.; Kniep, R.; Zahn, D. Atomistic Simulation Study of the Order/Disorder (Monoclinic to Hexagonal) Phase Transition of Hydroxyapatite. Chem. Mater. 2005, 17, 1978−1981. (56) Vollmayr, K.; Kob, W.; Binder, K. Cooling-Rate Effects in Amorphous Silica: A Computer-Simulation Study. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 15808. (57) Sato, K.; Kogure, T.; Iwai, H.; Tanaka, J. Atomic-Scale {101¯ 0} Interfacial Structure in Hydroxyapatite Determined by HighResolution Transmission Electron Microscopy. J. Am. Ceram. Soc. 2002, 85, 3054−3058. ́ sarczyk, A.; Paszkiewicz, (58) Rapacz-Kmita, A.; Paluszkiewicz, C.; Sló Z. FTIR and XRD Investigations on the Thermal Stability of Hydroxyapatite during Hot Pressing and Pressureless Sintering Processes. J. Mol. Struct. 2005, 744−747, 653−656. (59) Ramesh, S.; Tan, C. Y.; Tolouei, R.; Amiriyan, M.; Purbolaksono, J.; Sopyan, I.; Teng, W. D. Sintering Behavior of Hydroxyapatite Prepared from Different Routes. Mater. Eng. 2012, 34, 148−154. (60) Bertinetti, L.; Tampieri, A.; Landi, E.; Ducati, C.; Midgley, P. A.; Coluccia, S.; Martra, G. Surface Structure, Hydration, and Cationic Sites of Nanohydroxyapatite: UHR-TEM, IR, and Microgravimetric Studies. J. Phys. Chem. C 2007, 111, 4027−4035.

J


Surface Structure of Hydroxyapatite from Simulated Annealing

Recommend Documents