Molecular Dynamics Exploration of Ordered-to-Disordered Surface

Article Cite This: J. Phys. Chem. C 2018, 122, 6691−6703

pubs.acs.org/JPCC

Molecular Dynamics Exploration of Ordered-to-Disordered Surface Structures of Biomimetic Hydroxyapatite Nanoparticles Qing Xie,† Zhiyu Xue,† Haojie Gu,† Cuifang Hu,† Mingli Yang,‡,§ Xin Wang,†,* and Dingguo Xu†,‡,* College of Chemistry, MOE Key Laboratory of Green Chemistry and Technology, and ‡Genome Research Center of Biomaterial, Sichuan University, Chengdu, Sichuan 610064, P. R. China § Institute of Atomic and Molecular Physics, MOE Key Laboratory of High Energy Density Physics and Technology, Sichuan University, Chengdu, Sichuan 610065, P. R. China Downloaded via UNIV OF TEXAS AT EL PASO on October 21, 2018 at 06:34:51 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: Hydroxyapatite (HAP, Ca10(PO4)6(OH)2) constitutes one of the main components in human bones and teeth. Investigations of the surface structures of HAP could help us to understand the interactions between the implant HAP-based biomaterials and their environment. We, in this work, systematically explored the surface structures of HAP based on two different force fields, Born−Mayer−Huggins (BMH) parameters and interface force fields (IFFs). Their different performances in constructing HAP ordered-to-disordered surface structures were compared and discussed. The annealed HAP (001), (010), (100), and (110) surface structures were obtained, and their thickness values of the disordered layers have been computed. The thicknesses have the same order of (001) < (010) < (100) < (110) even at different annealing temperatures by the two force fields. At the annealing temperature 1800 K, the theoretical prediction for the average disordered layer thickness values is 8.9 Å (BMH) and 10.3 Å (IFF), which are close to the experimental measurement of 10−20 Å. Further studies on the surface models of (010) and (100) with different ion-rich terminations show that the disordered layer thickness of the P-rich surface are systematically larger than that of the Ca-rich surfaces. In the surface structure analyses, total pair distribution functions, static structure factors, and the coordination numbers for Ca or P were performed. These analyses show that our present HAP surface models have obvious ordered-to-disordered topology. In the present molecular dynamics exploration on the surface structures of HAP, several new surface models with different terminations for HAP-based biomaterials have been suggested, which may be appropriate for further investigation of interactions between biomolecules and HAP material.

1. INTRODUCTION It has been widely accepted that synthetic calcium phosphates containing ceramics can be used as bone grafts because they have shown excellent biocompatible, osteoinductive, nontoxic, and biodegradable properties. 1 Hydroxyapatite (HAP, Ca10(PO4)6(OH)2)2−4 constitutes the major component of both bone mineral and tooth enamel.5,6 Basically, the preparation of HAP7 is usually sintered at high temperature (lower than the melting point) and cooled down to room temperature, which can be directly and safely implanted into the living body.8−10 During the bone tissue repairing, the importance of understanding interactions of protein-implanted ceramics has been addressed, esp., in the early stage of bone mineralization. Indeed, the growth or inhibition of bone or teeth is related to the biomolecule−HAP interaction at their interface.11,12 Furthermore, it has been realized that the tissueinductive materials (TIMs) is one of the major development directions for the next-generation biomaterials.13−16 Therefore, a particular surface structure bearing tissue inductivity should be critical in future material design. The characterizations of the © 2018 American Chemical Society

surface structure of biomaterials are thus to be one of the crucial issues in bone tissue regeneration and engineering. Currently, HAP particles are mainly obtained by various synthesizing methods,7 and HAP crystals with bioactive and higher crystallinity are obtained through sintering HAP particles at a certain temperature.8 In addition, another method of obtaining HAP crystals is to sinter natural HAP from animal bones at high temperature to remove organic matter and thus be safely applied to the human body.17 Experimentally, the structure and the composition of the HAP surface can be investigated and characterized by versatile instrumental methods, such as nuclear magnetic resonance,18−20 ultrahighresolution transmission electron microscopy (UHR-TEM),21,22 atomic force microscopy,23 X-ray photoelectron spectroscopy,24−26 time-of-flight secondary ion mass spectrometry,24,25 scanning Auger microscopy,27 scanning electron microscoReceived: January 6, 2018 Revised: March 6, 2018 Published: March 8, 2018 6691

DOI: 10.1021/acs.jpcc.8b00178 J. Phys. Chem. C 2018, 122, 6691−6703

Article

The Journal of Physical Chemistry C

Figure 1. Crystal structure and annealed structures at 1400 K based on different force fields for HAP Ca-rich surfaces, and (a−d) represent HAP (001), (010), (100), and (110), respectively. Part I: unannealed HAP structures. Part II: annealed HAP surface structures in BMH. Part III: annealed HAP surface structures in IFF. For clarity, we have labeled vacuum, surface layer, middle layer, and interior layer at the crystal bulk HAP model. Calcium atoms are in green color, oxygen in red, hydrogen in white, and phosphorous atom in purple.

py,24,27 X-ray diffraction (XRD),28,29 and Fourier-transform infrared.30,31 These observations demonstrate the great difference between the sintered HAP surface and the bulk HAP crystal. Particularly, Jäger et al.18 found that the crystalline HAP core is covered by a disordered surface layer with about thickness of 1 nm, and the Ca/P ratio of the surface area is smaller than the HAP core area. Bertinetti et al.21 also investigated the surface features of nanosized HAP by UHRTEM and provided evidence that HAP nanoparticles are constituted by a crystalline core and an amorphous surface layer with 1−2 nm thickness. People22,32 also found that the HAP surface structures from highly to poorly ordered status largely depend on preparation conditions, such as temperature, presence of anionic/cationic substituents, and various biomolecular exposures. Nevertheless, the characterization of biomolecule binding affinity or specificity on an inorganic material is still hindered by the current experimental limit. Theoretical computations, on the other hand, have been extensively applied in the studies of interactions between biomolecules and nonbiological materials, esp., materials such as HAP. First principle calculations and molecular dynamic (MD) simulations have been reported. Majority efforts have been applied to investigate interactions between the different HAP crystal surfaces with water,33−35 organic molecules,36,37 amino acids,38−40 collagen,41 bone morphogenetic proteins (BMPs),42−45 or other proteins.46−49 According to lots of experimental results,22,50−52 HAP (001), (010), (100), and (110) usually adopted as the model

substrates of bimolecular adsorptions. For instance, Wu et al.44,45 investigated the interactions between BMP-7 and the HAP (001) surface. Interfacial actions between the HAP (001), (010), (100) planes and polyacrylic acid via the MD simulation method by Bhowmik et al.37 By the MD code, Almora-Barrios et al.39 found that the Hyp−Pro−Gly peptide in collagen prefers to interact with HAP (110) rather than with the HAP (001) plane. However, it should be emphasized that the HAP model used in these studies are mostly based on the model directly cut from the crystal structure, which is clearly inconsistent with the experimental proposal that HAP with the amorphous surface layer has bioactivity.53−55 Moreover, the importance of the biomaterial surface properties has been realized to be one of the key issues in the development of TIM. Therefore, a suitable HAP model comprising different interiors and boundary regions as the ordered/disordered microstructure is highly desired. We56 have proposed one strategy to construct such kind of model using a simulated annealing (SA) algorithm, in which the HAP (100) surface model is described using the force field parameters reported by Hauptmann et al.57 We then extend our strategy to more surface models for HAP. Four different surfaces, (001), (110), (010), and (100) were considered here. Our study could thus provide a possible model with an amorphous surface layer for studying the interactions between HAP-based biomaterials and biomolecules. Basically, to correctly describe the HAP model and corresponding properties in the biological process, accurate 6692


Article


Figure 2. Same as Figure 1, but the annealing temperature is set as 1600 K.

OH− in the unit cell.65 Its detailed lattice parameters can be found in Table S1 of the Supporting Information. The slab models of the HAP (001), (110), (010), and (100) surfaces were sliced from the bulk HAP crystal. It has to been pointed out that although (100) plane is isostructural with (010), distinct amorphous surface properties have been identified by experimental group.22 Therefore, it would be also interesting to further examine their properties upon the SA MD approach. In our simulation, each model was added a 40 Å vacuum layer. This is to simulate the surface effectively and avoid interaction between atomic layers under the threedimensional periodic boundary condition. Followed by our previous work,56 the obtained models were putatively divided into three layers, that is, interior, middle, and the surface, which are plotted in Figures 1 and 2 for four different crystal surface directions. The interior part was fixed to mimic the real inner core, whereas the middle and surface parts were freely relaxed during the simulation to check their response under hightemperature annealing simulation. In particular, the initial length for the surface layer was set to be about 20 Å because of the experimentally observed 10−20 Å amorphous layer. The original thickness values for the total and three layers of all HAP surface models are summarized in Table S2. Another important issue is that different terminations of the HAP surface could lead to different results. In fact, our strategy to build the slab model could result in either Ca-rich or P-rich surface model. Hence, to test the effects of different termination patterns of the HAP surface, we have then included models with both termination types in the simulation, but only for (010) and (100) surfaces.

force field should be one of the key issues because classical MD simulation based on the force field approach would be the only choice to tackle the large size system like HAP. Recently developed interface force fields (IFFs)58 were reported, which is specifically designed to understand interactions between the HAP surface and proteins. As suggested by the previous work, Born−Mayer−Huggins (BMH) parameters are particularly suitable for investigation of issues related with the bulk HAP model.59 Therefore, it would be interesting to extend our simulation strategy using other force field, such as IFF. On the other hand, it is also well-known that the HAP surface has different terminations based on the Ca/P ratio. Indeed, features of these surfaces were systematically examined by many groups.33,60−63 More importantly, there is no direct experimental evidence to support which kind of surface should be present when the material is implanted into human body. In addition, a correct HAP surface model system is particularly useful in understanding the interactions between the HAP material and its biological environment. In this paper, we will try to extend our strategy to four different surface directions and compare performances of BMH and IFF parameters in constructing the surface structures of HAP via the MD SA method.64 In particular, to get more general pictures, (010) and (100) with different Ca/P terminations will be examined here.

2. COMPUTATIONAL DETAILS It has been well-documented about the geometry properties for the HAP crystal structure. Basically, the hexagonal HAP crystal is in the P63/m space group with 10 Ca2+ ions, 6 PO43−, and 2 6693


Article


Figure 3. Calculated IR spectra based on annealed models via BMH (a) and IFF (b) parameters at 1400 K. For simplicity, the IR spectra using the HAP (100) surface model were displayed here, and experimental spectra76,77 at 1100 °C are also included for comparison.

All MD simulations were carried out using the LAMMPS suite of program66 in a canonical (NVT) ensemble with the Nose−Hoover thermostat,67 and the Verlet leapfrog algorithm was applied to integrate the equations of motion.68 All HAP surface models were relaxed for 50 ps at an initial temperature of 300 K, and then the target temperature was reached at the heating rate of 10 K/ps. The integration time interval is set as 0.1 fs. Practically, the sintering temperature is usually around 1400−1600 K. However, in our simulation, several target temperatures (lower, around, and above the melting point of HAP, 1923.16 K69) are selected to examine the temperature dependence for different force fields. The much high annealing temperatures around or over the melting point can accelerate the surface structure reconstruction. Basically, after 200 ps equilibration at the target temperature, the system was cooled step by step with 100 K per 20 ps until it reached 300 K. In this cooling process, a 50 ps MD relaxation will be carried out every 100 K. Finally, the system was further equilibrated at 300 K for 50 ps to obtain the annealed HAP surface structure. The trajectory relaxed at 300 K was saved for final analysis. In this work, atoms in the HAP model were described using two force fields, that is, BMH and IFF. Their corresponding potential energy forms are summarized in the Supporting Information with full set of parameters. The computational strategy for the thickness value can be found in Figure S1. To examine the viability of our models after annealing, we then calculate the infrared (IR) spectrum based on Fourier transform of the velocity correlation function.37,70,71 The IR spectra were computed for both HAP crystal (a 2 × 2 × 2 supercell model) and annealed structures. Once we have constructed an annealed HAP surface model, several structural analysis methods were applied to understand structural properties in this work, for example, the pair distribution function (PDF), static structure factor, and coordination number.72,73 The PDF can reflect the probability of other particles appearing from an arbitrarily designated “center” particle to the location of the radius r, which can measure the atomic distribution ordering and describe the atomic distribution. The PDF is defined as g (r ) =

1 Nρ0

N

∑ i=1

gαβ (r ) =

∑ i=1

niβ (r ) 4πr 2Δr

(2)

The static structure factor is one of the parameters to connect the experimental data and the simulations. For diffraction analysis, it represents the scattering ability of materials to radiation and reflects the average information of structures. The static structure factor can be obtained by Fourier transformation of PDF. The partial static structure factor Sαβ(q) between the atom α and atom β is defined as Sαβ(q) = δαβ &4πρ0 (cαcβ)1/2

∫0

∞

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

in which δαβ is the Kronecker delta function and cα = Nα/N. The total static structure factor S(q) is then determined by S(q) =

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

α ,β

In the characterization of microstructure, the coordination number is used to describe the average atoms number in the first shell of the central atom and reflects the central atom with other atoms binding ability and coordination relationship. In general, the coordination number is obtained by integrating the first peaks of PDF. The number distribution of atom β around atom α is defined as nαβ (R ) = 4πρβ

∫0

R

gαβ (r )r 2 dr

(5)

where R is the radius of the first shell of atom β around atom α. ISAACS program74 was used to compute the PDF and the static structure factor S(q) for the HAP crystal and each annealed surface structure.

3. RESULTS AND DISCUSSION 3.1. IR Spectra. The validation of two parameters has been systematically done by several groups for many fundamental properties of HAP, for example, crystal cell unit parameters,57 IR spectra,58 melting point mean square displacement, elastic modulus,75 cleavage energies, immersion energies,58 and so forth. As we pointed above, the surface of the sintered HAP appeared in the disorder layer of 1−2 nm (10−20 Å), but the core part is nearly in its crystal form.18,21,22 A SA approach has been suggested to be a suitable method to mimic the mineralization process. Indeed, in our previous work, we have constructed a three-layer structure for the HAP (100) plane based on the BMH force field. Here, we will further examine

ni(r ) 4πr 2Δr

Nα

1 Nαρβ

(1)

where N is the number of atoms in the slab model, ρ0 is the number density, and ni(r) is the number of particles found in the shell from r to r + Δr from a central atom i. In our work, the partial PDF is used to measure the distribution of atom β around atom α, which is defined as 6694


Article

The Journal of Physical Chemistry C two different force field parameters in constructing HAP models with disordered surface layer and ordered core. Before we further discuss the models via BMH and IFF potentials, we calculated the IR spectra of the annealed HAP model at 1400 K using HAP (100) as the template system. As displayed in Figure 3, the overall shapes of IR spectra of HAP by BMH and IFF are in good agreement with experimental spectra obtained at 1100 °C.76,77 In particular, the peak around 3600 cm−1 should be assigned to be the stretching vibration of OH− experimentally. For the spectra calculated based on BMH parameters, such vibrational peaks occur around 3500 cm−1, whereas nearly 3600 cm−1 in the IFF case can be observed.58 More importantly, relatively cleaner peak structure on this IR spectra obtained using the IFF parameter might suggest that IFF parameters have better performance than BMH parameters in describing intramolecular vibrations in HAP. On the other hand, a relatively wide peak for P−O stretching around 1000− 1100 cm−1 can be seen from the figures. In particular, more than 200 cm−1 difference can be found for BMH potential and around 50 cm−1 for IFF parameters for both crystal and annealed structures. Similar results can be found for the O−P− O bending mode occurring around 500−600 cm−1. Nevertheless, overall IR spectra shapes are still similar between calculated and experimental measurements. Current IR results could ensure us that both BMH and IFF parameters might be suitable for further construction of the layered HAP model in the annealing simulation. 3.2. Ordered-to-Disordered HAP Models. We have in the last work built a so-called three-layered HAP models based on BMH force field parameters, but only for Ca-rich HAP (100).56 The detailed pair distributions and static structure factors at 2400 K have also been systematically analyzed. The middle layer represents a gradual zone for the change from a crystal status to a chaotic state. However, as we have mentioned above, the usual sintering temperature for the synthesis of HAP is about 1400−1600 K. Hence, we first set the target temperatures at 1400 and 1600 K here in our SA computations. Similarly, only those surfaces with Ca-rich termination were considered in this section. Of course, such temperatures cannot cause the HAP crystal surface zone to be completely melted, but they should be enough to cause some chaos for ions of HAP at the surface layer. The initial and annealed structures of HAP (001), (010), (100), and (110) in two types of force fields with Ca-rich termination are displayed in Figures 1 and 2, respectively, and corresponding thickness values of the disordered layer for all plane directions are listed in Table 1. Indeed, although the

1400 or 1600 K. Experimentally, it is impossible to identify which surface direction is dominant; thus, we can simply give an average thickness values for the disordered layer. The average thickness values for the disordered layer over four planes for the HAP slab model in BMH parameters are calculated to be 7.5 Å (1400 K) and 8.1 Å (1600 K). Interestingly, HAP models based on IFF parameters have the same thickness layer sequence. However, the thicknesses for those disordered layers are systematically larger than those for BMH-based models. The difference becomes bigger with increasing starting annealing temperature. The average thickness values for the disordered layer over four planes in IFF parameters are calculated to be 8.9 Å (1400 K) and 10.3 Å (1600 K). However, we should also emphasize that previous studies have suggested that there are stoichiometric and nonstoichiometric terminations for HAP surfaces. In fact, theoretical simulation has suggested that stoichiometric surface such as (010) is not the most stable termination compared to its nonstoichiometric counterpart.33 More importantly, (010) stoichiometric surface was suggested to react with water spontaneously.34,78 Such results might indicate that we should not simply use a stoichiometric surface model in the investigation of interactions between biomolecules and HAP materials. On the other hand, among the nonstoichiometric terminations, according to different Ca/P ratios, we can have so-called Ca-rich and P-rich surfaces. Recently, Chiatti et al.60 have investigated the CO adsorption mechanism on the HAP nanoparticle surface structure using the first principle method. In addition, their different performances in the interactions with proteins have been addressed by several groups using various simulation approaches.79−81 Obviously, to get a complete understanding of the surface structure, different surface terminations should be included here. 3.3. Termination and Temperature Dependence. In this work, the effects of different terminations in constructing the ordered-to-disordered HAP surface model were also examined, esp., the chaotic degrees of surface layers. HAP (100) and (010) surfaces were selected as the templates. Moreover, to tackle the temperature effects to the amorphous surface layer, several starting annealing temperatures (1400, 1600, 1800, and 2400 K) are included here. All MD setup protocols are essentially the same as above. We have to point out that the preparation of different termination HAP surfaces has several choices. Here, we simply used the way of “cutting” from the crystal parallel to the surface direction to obtain corresponding slab models. Therefore, we can easily have either more Ca ions or phosphate groups on the surface exposed to vacuum. For the HAP (100) Ca-rich and P-rich surfaces, the Ca/P ratios were calculated to be 1.71 and 1.60, respectively. In addition, the computed Ca/P ratios are 1.83 and 1.60 for HAP (010) Ca-rich and P-rich surfaces, respectively. For clarity, we simply define (010) Ca-rich surface as (010) + Ca, and (010) + P for the P-rich surface, as well as (100) surface. The surface models with different ion-rich patterns are then given in Figure 4. In general, the thickness values for Ca-rich surfaces are systematically thinner than those for P-rich surfaces from Table 2. For example, in the BMH case, when the annealing temperature is 1400 K, about 5.1 Å for (010) + Ca and 5.7 Å for (100) + Ca surface. However, much larger values of 9.9 Å for (010) + P and 10.1 Å for (100) + P surface models can be found. Same situation is also for the models described using

Table 1. Calculated Thickness Values (Å) of Disordered Surface Layers for HAP Surface Models BMH

IFF

Temp/K

(001)

(010)

(100)

(110)

(001)

(010)

(100)

(110)

1400 1600

4.4 5.0

5.7 5.8

9.9 10.3

10.1 11.3

5.9 6.7

7.0 7.5

10.7 11.8

12.0 15.2

annealing temperature (1400 or 1600 K) is much lower than the melting point of HAP, obvious ordered-to-disordered topology can be found for every surface model but the thickness degree of the disordered layer shows different performances. In BMH cases, the thicknesses of the disordered layers have the order of (001) < (010) < (100) < (110) according to Table 1, no matter the annealing temperature is 6695


Article


increment from 1400 to 2400 K, whereas as large as 7.1 Å increment for the (100) + Ca surface. More interestingly, for temperatures around the melting point (1923.16 K), for example, 1800 K, HAP models in both force fields can still maintain quite stable, but with obvious increasing chaos in the surface zone no matter in BMH or in IFF parameters. From Table 2, average thickness value for the models described in the BMH form is calculated to be 8.9 Å, whereas 10.3 Å for the models in IFF parameters. Of course, different performances of two force fields can be identified, but the thickness values are systematically larger for the models in IFF parameters than those described in BMH parameters in all cases. To understand this, we should take a close look at the detailed parameters of two force fields. Basically, it has been suggested that BMH force field parameters are particularly suitable to simulate the bulk HAP model.57,59 Hence, current results should be understandable that the surface layers described using BMH parameters are not entirely disordered or quite thin. However, because the IFF parameters are particularly designed for the simulation of interface of the HAP model, it could be expected that the more flexibility can be observed under annealing simulations. More importantly, the charge on Ca2+ ion is set to be around 1.5 e, which is smaller than 2.0 e in the BMH form. Such difference could obviously cause weaker electrostatic interactions involved Ca2+. In fact, people have realized that the charge on surface calcium ion is different from that in the inner core.82 Moreover, 9−6 L-J potential used in IFF could also indicate much stronger repulsion of van der Waals force. The stability of the HAP model described using IFF might not be as stable as the case in BMH, esp., for the part close to the surface. Therefore, in comparison with the BMH model, a little bit more chaos in the surface layer based on IFF parameters should be reasonable. Despite results presented above, we have to point out that once the annealing temperature is much over the melting point such as 2400 K, dramatic different performances of two force fields can be seen according to Table 2. Clearly, for both (010) + P and (100) + P planes in IFF parameters, the overall HAP models become completely chaotic and instable. However, models described in the BMH form maintain very well except (100) + P. In other words, the thickness values of disordered layers for the P-rich surface are systematically larger than those for Ca-rich surfaces, which might suggest that the Ca-rich surface is more stable. Our results are not inconsistent with previous theoretical studies that the Ca-rich surface is much more stable than the P-rich surface.33 Considering the experimental estimation18 of 10−20 Å for the amorphous layer of HAP, it would be more suitable to construct this kind of HAP model around temperature closer to melting point. 3.4. Surface Structure Analysis. Except the different thicknesses of the disordered layer, some common features can be observed for both parameters. For the HAP annealed structures, the locations of the ions (Ca2+, OH−, and PO43−) in

Figure 4. HAP (010) and (100) surface models with different ion-rich patterns from side view. Color codes are the same as Figure 1, except PO3− is labeled in a triangle shape.

IFF parameters, but systematically more chaotic representation for both Ca-rich and P-rich terminations can be located. In particular, about 5.7 and 7.0 Å are calculated for (010) + Ca and (100) + Ca models, respectively, whereas 10.7 and 11.4 Å for (010) + P and (100) + P surfaces, respectively, when the annealing temperature is 1400 K. Nearly the same trends can be found in the case obtained at 1600 K, although a little change of thickness values of disordered layer can be seen. It should be pointed out that there is no experimental evidence to say which kind of plane or termination when the material is implanted into human body. In other words, the experimental measurements of thickness values of the surface disordered layer might be an average of all possible surface directions. Interestingly, the average values are calculated to be 7.7 Å (1400 K) and 8.0 Å (1600 K) for the models in BMH parameters, whereas 8.7 Å (1400 K) and 9.4 Å (1600 K) for the models in IFF parameters. Overall, it seems that the surface area is more stable for the models described in the BMH form than in the IFF form. On the other hand, it has been found by experiments that there are about 10−20 Å disordered surface layer for nanocrystalline HAP.18,21 Clearly, the simulated disordered surface layer at low annealing temperatures (1400 and 1600 K) cannot fully reproduce experimental measurements. Therefore, we further extend our annealing temperature to much higher values around or much over the HAP melting point. We have in the last work built a so-called three-layered HAP models based on BMH force field parameters, but only a Ca-rich HAP (100) surface at 2400 K was employed in our work.56 A chaotic surface layer less than 10 Å was obtained. In this work, two annealing temperatures of 1800 and 2400 K are further used for the evaluation of temperature dependence. According to Table 2, we can easily find that the disordered layers become thicker with increasing annealing temperatures. The anisotropic characteristics of HAP are also clearly shown by different temperature dependences. Typically, for the models in IFF parameters, for the (010) + Ca surface, there is only 2.8 Å

Table 2. Calculated Thickness Values (Å) of Disordered Surface Layers for HAP Surface Models BMH

IFF

Temp/K

(010) + Ca

(100) + Ca

(010) + P

(100) + P

(010) + Ca

(100) + Ca

(010) + P

(100) + P

1400 1600 1800 2400

5.1 5.4 5.9 6.6

5.7 5.8 6.0 7.1

9.9 10.3 11.7 12.5

10.1 10.3 11.9 24.7

5.7 6.1 6.7 8.5

7.0 7.5 8.2 14.1

10.7 11.8 12.6 44.2

11.4 12.0 13.5 88.0

6696


Article

The Journal of Physical Chemistry C the middle layers have tiny changes, whereas about one-half of hydroxyl groups deflect about 90° in BMH force field and much less than 90° using IFF parameters, as displayed in Figure 5. Such kind of deflection phenomenon of hydroxyls in HAP

Figure 5. (a) Unannealed structure and (b) annealed structures of HAP (100) in the BMH force filed and (c) in IFF parameters. Color codes for atoms are the same as Figure 1.

has been reported in some literatures.38,40,83 To illustrate motions of Ca2+ or PO43−, we then plotted the surface topology in Figures 6 and 7 from the top view. The phosphate groups

Figure 7. Same as Figure 6, except the structures are for HAP (010) + P (left) and (100) + P (right) from top view.

definitely decrease the surface energy and finally stabilize the surface zone, which is confirmed in our last work.

Figure 6. Structures of HAP (010) + Ca (left) and (100) + Ca (right) from top view, (a) unannealed, (b) annealed at 1600 K in BMH, and (c) annealed at 1600 K in IFF. Color codes are the same as Figure 1, except PO3− is labeled in a triangle shape.

Figure 8. Annealed structures of HAP (010) + Ca (left) and (100) + Ca (right) from side view, (a) annealed at 1600 K in BMH and (b) annealed at 1600 K in IFF. Color codes are the same as Figure 1, except PO3− is labeled in a triangle shape.

seem to stay in their original positions but undergo some rotations. Most of the calcium ions, on the contrary, left their original positions, and some cavities could come up. This could definitely lead to the irregular coordination state around either calcium ions or phosphate groups. Meanwhile, such rearrangements of ions should result the charge density distribution of the annealed surface different from the pure crystal plane and cause some roughness on the surface. Besides, the calcium ions at the surface and core parts are not in the same chemical environment. Surface calcium ions of crystal HAP are in the unsaturated coordination status. However, during the annealing simulation because OH− is much lighter than Ca2+ and PO43−, drastic movements can be observed. For the model subjected to annealing simulation, we could find that some hydroxyl groups move from inside to the surface and finally cover the surface calcium ions, as shown in Figures 8 and 9. Such behavior could

Figure 9. Same as Figure 8 with the annealed structures of HAP (010) + P (left) and (100) + P (right) from side view. 6697


Article


Figure 10. Total PDFs of the crystal and annealed structures for the HAP Ca-rich models at 1600 K in BMH and IFF types force fields, (a) HAP (010) and (b) HAP (100).

Figure 11. Total structure factor [S(q)total] and the partial static structure factor of O−Ca [S(q)O−Ca], P−Ca [S(q)P−Ca], and P−O [S(q)P−O] of the three layers of the HAP (100) Ca-rich structure annealed at 1600 K in the IFF force field, (a) S(q)total, (b) S(q)O−Ca, (c) S(q)P−Ca, and (d) S(q)P−O.

Figure 12. Same as Figure 11, but the force field is BMH. 6698


Article


Figure 13. Coordination number distraction of P, Ca, and O atom in the three layers of HAP (100) Ca-rich structures annealed at 1600 K in the IFF force field, (a) P−Ca, (b) Ca−P, (c) O−Ca, and (d) Ca−O.

plotted in Figures 11 and 12 as well. All other surface directions have nearly the same trend and are thus included in the supporting material (Figures S4−S6). Interestingly, overall shapes of the S(q) curves and peak locations of the three layers are almost the same for both BMH and IFF models. Although the bond length and bond angle between neighboring atoms in annealed models deviate a little bit from the crystal, the interaction of neighboring atoms still shows some regularity in the short distance range. Indeed, nearly identical peak positions and heights for the total static structure factor of S(q) around short and middle distance ranges (r < 7 Å) can be observed for all layers, which suggests an ordered atomic arrangements even in the surface layer. For the long distance zone (r > 7 Å), a little bit irregular occurrence of peaks can be observed for middle and surface layers compared with the crystal structure. To get more clear view of the structure, we further calculated and plotted partial static structure factors of S(q)O−Ca, S(q)P−Ca, and S(q)P−O in Figures 11 and 12 too. Clearly, curves only at short distance part (r < 3 Å) show nearly the same performance, which suggests that the atomic arrangement are relatively regular at the short distance zone. Significant differences between three layers can be located for the part of r > 3 Å, esp., for r > 7 Å. In particular, the curves of S(q)O−Ca and S(q)P−O for middle and surface layers in the r > 8 Å region become smooth and finally close to 0, but the surface layer is smoother than the middle layer. Generally, much different static structure factors calculated here can clearly support the formation of a threelayer HAP model after annealing simulation and thus reflect the disordering difference between the middle and surface layers. Once the ordered-to-disordered model via an annealing approach has been constructed for HAP, it would be worthwhile to characterize detailed surface structural patterns. Here, the coordination numbers for Ca or P deserve special attentions because they can be easily obtained via integration of partial PDF. Basically, the coordination number around the center atom could shed lights on the structural changes under simulation. Again, the HAP (100) surface model at 1600 K was taken as the template system to address the changes of

To further illustrate the conformational changes of the annealed structures from interior to surface layer, we then calculated the total PDFs for both HAP crystal and all annealed structures using an annealing temperature of 1600 K, as shown in Figure 10. The detailed pair distributions and static structure factors at 2400 K have also been systematically analyzed in our last work.56 The middle layer represents a gradual zone for changing from a crystal to a chaotic state. Notably, the positions of distinguished peaks of crystal and annealed structures in both BMH and IFF are nearly identical in the region of r < 3 Å. It means that the atomic arrangements of the annealed systems in short ranges are essentially regular and close to the crystal structure. Moreover, the first tiny peaks can be assigned to the hydroxyl groups, the second (about 1.5 Å) and the third peaks (about 2.4 Å) to the P−O and O−O bonds in phosphate groups, respectively. Clearly, the structures in short distance zone still bear their crystal properties with atomic arrangement is pretty regular. In the medium-range area (r = 3−7 Å), obvious difference from the interior part can thus be observed, which means significant atomic movements could occur during the simulation. Generally, the g(r) curves of the middle part are subjected to some fluctuations and with some small peaks. For the long-distance region (r > 7 Å), the values of pair distribution finally approach to 1, which indicates that the annealed surface structures in the larger r region might undergo some changes away from the crystal structure. In a word, total PDFs can partially confirm that the chaotic structures should exist in annealed structures but cannot tell which part is in chaos status. The static structure factor S(q) is another physical quantity to identify the structural information of the system, which can be compared directly with the results of XRD experiments. Here, S(q) is used to tell us the structural difference among three layers of the HAP annealed structures expressed in different force fields. We calculated the total static structure factor of S(q)total for different structure layers. Similarly, the annealed HAP (100) model at 1600 K was analyzed for both BMH and IFF cases, and corresponding performances are 6699


Article


coordination number between the interior and middle layers suggests BMH force field might not be well-suitable for constructing this kind of order-to-disorder surface model. On the contrary, models described using IFF parameters show clear thicker disorder layer than BMH, and gradual changes of coordination numbers from inside to outside structure. In fact, this IFF is specifically designed to account for the interaction occurring at the interface area of HAP. Combined all above information, we then tentatively suggest that IFF parameters might be more suitable for constructing the HAP amorphous surface model, which could be more appropriate for further investigation of interactions between biomolecules and HAP material.

coordination numbers. Not surprisingly, the coordination number of P−O is calculated to be 4, and 1 for H−O pair in all three layers, whereas other atom pairs in middle and surface layers change from the interior to surface layer. We have observed such phenomena in our previous work,56 when the annealing temperature was set to be as high as 2400 K. In this work, we only reported the coordination numbers of the P, Ca, and O atoms calculated using IFF parameters, which were then illustrated in Figure 13. The average coordination numbers for atom pairs are summarized in Table 3. For comparison, Table 3. Average Coordination Number of Three Layers for HAP (100) Annealed in BMH at 2400 K (Ref 56) and IFF at 1600 K Parametersa BMH (ref 56)

4. CONCLUSIONS Traditionally, people often applied slab models directly cutting from the HAP crystals to mimic the implanted materials, which are further used to investigate corresponding interactions with biomolecules. Such treatment apparently ignore the fact of an amorphous surface layer (a Ca/P ratio close to 1.5) identified experimentally. We have in the last work56 proposed one strategy to reconstruct the surface structure using a SA MD simulation. The Hauptmann’s potential energy functional form was applied, and a distinguished ordered-to-disordered structure along the HAP (100) direction can be established. However, the Hauptmann’s parameters57 are suggested to be suitable for the bulk structure, and their applicability to HAP interfacial characteristics is thus questionable. Recently, Lin and Heinz58 proposed another type of parameters that is specifically designed to simulate the process occurring at the material interfacial area. In this work, we then systematically explore the surface structures of HAP based on these two force fields. The annealed surface structures of HAP (001), (010), (100), and (110) were obtained, and their thickness values of the disordered layers have been calculated. In BMH cases, the thicknesses have the same order of (001) < (010) < (100) < (110) at the annealing temperatures 1400 and 1600 K. Interestingly, HAP models based on IFF parameters have the same thickness layer sequence. However, the thicknesses by IFF are systematically larger than those for BMH-based models. Moreover, to understand effects of temperatures and terminations with the nonstoichiometric Ca/P ratio in constructing surface models, (010) and (100) with different Ca and P-rich terminations have further been investigated at four different starting annealing temperatures here. The results show that the disordered layers become thicker with increasing annealing temperatures. The anisotropic characteristics of HAP are also observed by the different temperature dependences. In addition, our calculated results shown that the thickness values of disordered layers for the P-rich surface are systematically larger than those for Ca-rich surfaces, which indicates that the Ca-rich surface is more stable. Much different performances between two types of force fields have been found. The thickness values are systematically larger for the models in IFF parameters than those described in BMH parameters in all cases. The differences of detailed parameters of two force fields were compared and discussed to understand their different performances for HAP surface structure calculations. Interestingly, at 1800 K, a temperature close to the melting point of HAP, surface models in both force fields can still maintain quite stable, but with obvious increasing chaos in the surface zone no matter in BMH or in IFF parameters. At the temperature, the

IFF

A−B pair

interior

middle

surface

interior

middle

surface

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

9.00 5.40 3.00 7.80 4.00 1.00

9.00 5.40 2.97 7.72 4.00 1.00

8.24 4.12 2.51 5.51 4.00 1.00

9.00 5.40 3.00 7.80 4.00 1.00

8.95 5.35 2.80 7.68 4.00 1.00

8.24 4.94 2.63 6.98 4.00 1.00

a

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

coordination numbers obtained at 2400 K using BMH parameters are also included. Obviously, dramatic changes for numbers of calcium ions around the phosphate group can be observed from inside to outside. For example, 100% of the phosphate groups in interior layers are surrounded by nine Ca atoms because this layer is fixed and pretty much the same as a HAP crystal structure. However, in middle layers, because of the increment of chaos, two types of coordination for P−Ca can be seen, that is, nearly 95% of the P atoms are coordinated with nine Ca atoms surrounded and 5% with eight Ca atoms. This could lead to a slightly smaller averaged coordination number of 8.95. However, according to Table 3, nearly no change of coordination number for the annealed model in BMH parameters can be observed, which clearly means that a very stable bulk structure described in BMH parameters under SA simulation. In addition, according to Figure 13 in this work and previous work in ref 56 for models obtained in IFF and BMH parameters, much more coordination possibilities for P− Ca pair can be found in the surface layer. For example, in the IFF case, about 17.4, 7.9, 7.9, and 66.8% of the phosphate groups are coordinated with six, seven, eight, and nine Ca atoms, respectively, which results in an average coordination number of 8.24 for P−Ca pair. Such kind of change of coordination number should be reasonable because some drastic movements of ions at the ceramic surface must occur under high temperature and thus cause natural disordering for those atoms at the surface layer. Same phenomena can be observed for other atom pairs, as shown in Figure 13. Overall, even at higher temperature, systems undergo tiny deviations from the crystal structure for the model described using BMH parameters. This is consistent with the disorder thickness calculations presented above. Physically, it would be natural to believe that when the HAP ceramics is under sintering, particles arrangement away from its original positions might change gradually from the core to the surface zone. Nearly identical 6700


Article


(8) 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. Des. 2012, 34, 148−154. (9) Champion, E. Sintering of Calcium Phosphate Bioceramics. Acta Biomater. 2013, 9, 5855−5875. (10) Kamalanathan, P.; Ramesh, S.; Bang, L. T.; Niakan, A.; Tan, C. Y.; Purbolaksono, J.; Chandran, H.; Teng, W. D. Synthesis and Sintering of Hydroxyapatite Derived from Eggshells as a Calcium Precursor. Ceram. Int. 2014, 40, 16349−16359. (11) Schachschal, S.; Pich, A.; Adler, H.-J. Aqueous Microgels for the Growth of Hydroxyapatite Nanocrystals. Langmuir 2008, 24, 5129− 5134. (12) Cerroni, L.; Filocamo, R.; Fabbri, M.; Piconi, C.; Caropreso, S.; Condò , S. G. Growth of Osteoblast-like Cells on Porous Hydroxyapatite Ceramics: an in Vitro Study. Biomol. Eng. 2002, 19, 119−124. (13) Mouriño, V.; Cattalini, J. P.; Roether, J. A.; Dubey, P.; Roy, I.; Boccaccini, A. R. Composite Polymer-Bioceramic Scaffolds with Drug Delivery Capability for Bone Tissue Engineering. Expert Opin. Drug Delivery 2013, 10, 1353−1365. (14) Pérez, R. A.; Won, J.-E.; Knowles, J. C.; Kim, H.-W. Naturally and Synthetic Smart Composite Biomaterials for Tissue Regeneration. Adv. Drug Deliv. Rev. 2013, 65, 471−496. (15) Nakamura, M.; Hiratai, R.; Hentunen, T.; Salonen, J.; Yamashita, K. Hydroxyapatite with High Carbonate Substitutions Promotes Osteoclast Resorption through Osteocyte-like Cells. ACS Biomater. Sci. Eng. 2016, 2, 259−267. (16) Ruhé, P. Q.; Kroese-Deutman, H. C.; Wolke, J. G. C.; Spauwen, P. H. M.; Jansen, J. A. Bone Inductive Properties of rhBMP-2 Loaded Porous Calcium Phosphate Cement Implants in Cranial Defects in Rabbits. Biomaterials 2004, 25, 2123−2132. (17) Barakat, N. A. M.; Khalil, K. A.; Sheikh, F. A.; Omran, A. M.; Gaihre, B.; Khil, S. M.; Kim, H. Y. Physiochemical Characterizations of Hydroxyapatite Extracted from Bovine Bones by Three Different Methods: Extraction of Biologically Desirable HAP. Mater. Sci. Eng., C 2008, 28, 1381−1387. (18) 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. (19) 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. (20) 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−606. (21) 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. (22) Sakhno, Y.; Bertinetti, L.; Iafisco, M.; Tampieri, A.; Roveri, N.; Martra, G. Surface Hydration and Cationic Sites of Nanohydroxyapatites with Amorphous or Crystalline Surfaces: A Comparative Study. J. Phys. Chem. C 2010, 114, 16640−16648. (23) 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. (24) 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. (25) 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. (26) Bengtsson, Å.; Shchukarev, A.; Persson, P.; Sjöberg, S. A Solubility and Surface Complexation Study of a Non-stoichiometric Hydroxyapatite. Geochim. Cosmochim. Acta 2009, 73, 257−267.

theoretical predictions for the average disordered layer thickness values are 8.9 Å (BMH) and 10.3 Å (IFF), which are close to the experimental measurement of 10−20 Å. Further surface structure analyses show that the ions (Ca2+, OH−, and PO43−) in the middle layers have tiny changes, whereas about one-half of hydroxyl groups deflect about 90° in the BMH force field and much less than 90° using IFF parameters. The total PDFs were computed to further illustrate the conformational changes of the annealed structures from the interior to surface layer. Moreover, our calculations reveal that the coordination numbers at the surface layer are much different from the bulk structure, which is due to drastic motions of ions. According to these analyses, we suggest that IFF might be more suitable for constructing the HAP amorphous surface model than BMH. In present MD exploration on the ordered-to-disordered surface structure of HAP, several new surface models with different terminations for HAP-based biomaterials have been characterized to some details, which may be ready for further investigation of interactions between biomolecules and HAP material.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b00178. Details of the BMH and IFF potential functional forms and corresponding parameters; the PDFs of three-layer HAP annealed models; and the computational strategy of thickness value for the disordered surface layer (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (X.W.). *E-mail: [email protected] (D.X.). ORCID

Mingli Yang: 0000-0001-8590-8840 Dingguo Xu: 0000-0002-9834-8296 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by National Key Research and Development Program (no. 2016YFB0700801). Some of results described in this paper were obtained on the National Supercomputing Center of Guangzhou.

■

REFERENCES

(1) Dorozhkin, S. V.; Epple, M. Biological and Medical Significance of Calcium Phosphates. Angew. Chem., Int. Ed. 2002, 41, 3130−3146. (2) Mann, S. Molecular Recognition in Biomineralization. Nature 1988, 332, 119−124. (3) Mann, S. Molecular Tectonics in Biomineralization and Biomimetic Materials Chemistry. Nature 1993, 365, 499−505. (4) Simmer, J. P.; Fincham, A. G. Molecular Mechanisms of Dental Enamel Formation. Crit. Rev. Oral Biol. Med. 1995, 6, 84−108. (5) Vallet-Regí, M.; González-Calbet, J. M. Calcium Phosphates as Substitution of Bone Tissues. Prog. Solid State Chem. 2004, 32, 1−31. (6) Dorozhkin, S. V. Calcium Orthophosphates. J. Mater. Sci. 2007, 42, 1061−1095. (7) Rivera-Muñoz, E. M.; Velázquez-Castillo, R.; Cabrera-Torres, J. L. Morphological Analysis of Hydroxyapatite Particels Obtained by Different Methods. Mater. Sci. Forum 2010, 638−642, 681−686. 6701


Article


(48) Liao, C.; Zhou, J. Replica-Exchange Molecular Dynamics Simulation of Basic Fibroblast Growth Factor Adsorption on Hydroxyapatite. J. Phys. Chem. B 2014, 118, 5843−5852. (49) Makrodimitris, K.; Masica, D. L.; Kim, E. T.; Gray, J. J. Structure Prediction of Protein-Solid Surface Interactions Reveals a Molecular Recognition Motif of Statherin for Hydroxyapatite. J. Am. Chem. Soc. 2007, 129, 13713−13722. (50) Xin, R.; Leng, Y.; Chen, J.; Zhang, Q. A Comparative Study of Calcium Phosphate Formation on Bioceramics in Vitro and in Vivo. Biomaterials 2005, 26, 6477−6486. (51) Xin, R.; Leng, Y.; Wang, N. In Situ TEM Examinations of Octacalcium Phosphate to Hydroxyapatite Transformation. J. Cryst. Growth 2006, 289, 339−344. (52) Leng, Y.; Chen, J.; Qu, S. TEM Study of Calcium Phosphate Precipitation on HA/TCP Ceramics. Biomaterials 2003, 24, 2125− 2131. (53) Lukić, M. J.; Škapin, S. D.; Marković, S.; Uskoković, D. Processing Route to Fully Dense Nanostructured HAp Bioceramics: From Powder Synthesis to Sintering. J. Am. Ceram. Soc. 2012, 95, 3394−3402. (54) Veljović, D.; Jokić, B.; Jankovic-Castvan, I.; Smičiklas, I.; Petrović, R.; Janaćković, D. Sintering Behaviour of Nanosized HAP Powder. Key Eng. Mater. 2007, 330−332, 259−262. (55) Veljović, D.; Jokić, B.; Petrović, R.; Palcevskis, E.; Dindune, A.; Mihailescu, I. N.; Janaćković, D. Processing of Dense Nanostructured HAP Ceramics by Sintering and Hot Pressing. Ceram. Int. 2009, 35, 1407−1413. (56) Wu, H.; Xu, D.; Yang, M.; Zhang, X. Surface Structure of Hydroxyapatite from Simulated Annealing Molecular Dynamics Simulations. Langmuir 2016, 32, 4643−4652. (57) 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. (58) Lin, T.-J.; Heinz, H. Accurate Force Field Parameters and pH Resolved Surface Models for Hydroxyapatite to Understand Structure, Mechanics, Hydration, and Biological Interfaces. J. Phys. Chem. C 2016, 120, 4975−4992. (59) Xu, Z.; Yang, Y.; Wang, Z.; Mkhonto, D.; Shang, C.; Liu, Z.-P.; Cui, Q.; Sahai, N. Small Molecule-Mediated Control of Hydroxyapatite Growth: Free Energy Calculations Benchmarked to Density Functional Theory. J. Comput. Chem. 2014, 35, 70−81. (60) Chiatti, F.; Corno, M.; Sakhno, Y.; Martra, G.; Ugliengo, P. Revealing Hydroxyapatite Nanoparticle Surface Structure by CO Adsorption: A Combined B3LYP and Infrared Study. J. Phys. Chem. C 2013, 117, 25526−25534. (61) Jimenez-Izal, E.; Chiatti, F.; Corno, M.; Rimola, A.; Ugliengo, P. Glycine Adsorption at Nonstoichiometric (010) Hydroxyapatite Surfaces: A B3LYP Study. J. Phys. Chem. C 2012, 116, 14561−14567. (62) Slepko, A.; Demkov, A. A. First Principles Study of Hydroxyapatite Surface. J. Chem. Phys. 2013, 139, 044714. (63) Sakhno, Y.; Ivanchenko, P.; Iafisco, M.; Tampieri, A.; Martra, G. A Step toward Control of the Surface Structure of Biomimetic Hydroxyapatite Nanoparticles: Effect of Carboxylates on the {010} PRich/Ca-Rich Facets Ratio. J. Phys. Chem. C 2015, 119, 5928−5937. (64) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671−680. (65) Wang, Q.; Wang, M.-h.; Wang, K.-f.; Liu, Y.; Zhang, H.-p.; Lu, X.; Zhang, X.-d. Computer Simulation of Biomolecule-Biomaterial Interactions at Surfaces and Interfaces. Biomed. Mater. 2015, 10, 032001. (66) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1−19. (67) Nosé, S. A Unified Formulation of the Constant Temperature Molecular-Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (68) Verlet, L. Computer “Experiments” on Classical Fluids.I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (69) Lu, Y.-P.; Li, M.-S.; Li, S.-T.; Wang, Z.-G.; Zhu, R.-F. PlasmaSprayed Hydroxyapatite Plus Titania Composite Bond Coat for

(27) 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. (28) Chen, J.-Z.; Shi, Y.-L.; Wang, L.; Yan, F.-Y.; Zhang, F.-Q. Preparation and Properties of Hydroxyapatite-containing Titania Coating by Micro-arc Oxidation. Mater. Lett. 2006, 60, 2538−2543. (29) Tanaka, H.; Yasukawa, A.; Kandori, K.; Ishikawa, T. Surface Modification of Calcium Hydroxyapatite with Hexyl and Decyl Phosphates. Colloids Surf., A 1997, 125, 53−62. (30) Yasukawa, A.; Ouchi, S.; Kandori, K.; Ishikawa, T. Preparation and Characterization of Magnesium-calcium Hydroxyapatites. J. Mater. Chem. 1996, 6, 1401−1405. (31) Tanaka, H.; Futaoka, M.; Hino, R. Surface Modification of Calcium Hydroxyapatite with Pyrophosphoric Acid. J. Colloid Interface Sci. 2004, 269, 358−363. (32) Celotti, G.; Tampieri, A.; Sprio, S.; Landi, E.; Bertinetti, L.; Martra, G.; Ducati, C. Crystallinity in Apatites: How Can a Truly Disordered Fraction Be Distinguished from Nanosize Crystalline Domains? J. Mater. Sci.: Mater. Med. 2006, 17, 1079−1087. (33) Chiatti, F.; Piane, M. D.; Ugliengo, P.; Corno, M. Water at Hydroxyapatite Surfaces: the Effect of Coverage and Surface Termination as Investigated by All-electron B3LYP-D* Simulations. Theor. Chem. Acc. 2016, 135, 54. (34) Astala, R.; Stott, M. J. First-principles Study of Hydroxyapatite Surfaces and Water Adsorption. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 075427. (35) Zahn, D.; Hochrein, O. Computational Study of Interfaces between Hydroxyapatite and Water. Phys. Chem. Chem. Phys. 2003, 5, 4004−4007. (36) de Leeuw, N. H.; Rabone, J. A. L. Molecular Dynamics Simulations of the Interaction of Citric Acid with the Hydroxyapatite (0001) and (011̅ 0 ) Surfaces in an Aqueous Environment. CrystEngComm 2007, 9, 1178−1186. (37) Bhowmik, R.; Katti, K. S.; Katti, D. Molecular Dynamics Simulation of Hydroxyapatite-polyacrylic Acid Interfaces. Polymer 2007, 48, 664−674. (38) 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 (0110̅ ) Surfaces. Langmuir 2009, 25, 5018−5025. (39) Almora-Barrios, N.; de Leeuw, N. H. Modelling the Interaction of a Hyp-Pro-Gly Peptide with Hydroxyapatite Surfaces in Aqueous Environment. CrystEngComm 2010, 12, 960−967. (40) 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. (41) Almora-Barrios, N.; de Leeuw, N. H. A Density Functional Theory Study of the Interaction of Collagen Peptides with Hydroxyapatite Surfaces. Langmuir 2010, 26, 14535−14542. (42) 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. (43) Dong, X.-L.; Qi, W.; Tao, W.; Ma, L.-Y.; Fu, C.-X. The Dynamic Behaviours of Protein BMP-2 on Hydroxyapatite Nanoparticles. Mol. Simulat. 2011, 37, 1097−1104. (44) 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. (45) 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. (46) Chen, X.; Wang, Q.; Shen, J.; Pan, H.; Wu, T. Adsorption of Leucine-rich Amelogenin Protein on Hydroxyapatite (001) Surface through -COO- Claws. J. Phys. Chem. C 2007, 111, 1284−1290. (47) Liao, C.; Xie, Y.; Zhou, J. Computer Simulations of Fibronectin Adsorption on Hydroxyapatite Surfaces. RSC Adv. 2014, 4, 15759− 15769. 6702


Article

The Journal of Physical Chemistry C Hydroxyapatite Coating on Titanium Substrate. Biomaterials 2004, 25, 4393−4403. (70) Ikoma, T.; Yamazaki, A.; Nakamura, S.; Akao, M. Preparation and Structure Refinement of Monoclinic Hydroxyapatite. J. Solid State Chem. 1999, 144, 272−276. (71) Kleinhesselink, D.; Wolfsberg, M. The Evaluation of Power Spectra in Molecular-Dynamics Simulations of Anharmonic Solids and Surfaces. Surf. Sci. 1992, 262, 189−207. (72) Duan, G.; Xu, D.; Zhang, Q.; Zhang, G.; Cagin, T.; Johnson, W. L.; Goddard, W. A. Molecular Dynamics Study of the Binary Cu46Zr54 Metallic Glass Motivated by Experiments: Glass Formation and Atomic-level Structure. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 224208. (73) 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. (74) Le Roux, S.; Petkov, V. ISAACSInteractive Structure Analysis of Amorphous and Crystalline Systems. J. Appl. Crystallogr. 2010, 43, 181−185. (75) Ma, X.; Zhang, L.; Wu, H.; Zhang, X.; Yang, M. Computer Simulations on the Mechanical Behaviors of Biphasic Calcium Phosphates. J. Mol. Model. 2017, 23, 156. (76) Ooi, C. Y.; Hamdi, M.; Ramesh, S. Properties of Hydroxyapatite Produced by Annealing of Bovine Bone. Ceram. Int. 2007, 33, 1171− 1177. (77) Zhou, J.; Zhang, X.; Chen, J.; Zeng, S.; De Groot, K. High Temperature Characteristics of Synthetic Hydroxyapatite. J. Mater. Sci.: Mater. Med. 1993, 4, 83−85. (78) Corno, M.; Rimola, A.; Bolis, V.; Ugliengo, P. Hydroxyapatite as a Key Biomaterial: Quantum-mechanical Simulation of Its Surfaces in Interaction with Biomolecules. Phys. Chem. Chem. Phys. 2010, 12, 6309−6329. (79) Dong, X.-L.; Zhou, H.-L.; Wu, T.; Wang, Q. Behavior Regulation of Adsorbed Proteins via Hydroxyapatite Surface Texture Control. J. Phys. Chem. B 2008, 112, 4751−4759. (80) Libonati, F.; Nair, A. K.; Vergani, L.; Buehler, M. J. Mechanics of Collagen−Hydroxyapatite Model Nanocomposites. Mech. Res. Commun. 2014, 58, 17−23. (81) Huang, B.; Yuan, Y.; Ding, S.; Li, J.; Ren, J.; Feng, B.; Li, T.; Gu, Y.; Liu, C. Nanostructured Hydroxyapatite Surfaces-Mediated Adsorption Alters Recognition of BMP Receptor IA and Bioactivity of Bone Morphogenetic Protein-2. Acta Biomater. 2015, 27, 275−285. (82) de Leeuw, N. H. A Computer Modelling Study of the Uptake and Segregation of Fluoride Ions at the Hydrated Hydroxyapatite (0001) Surface: Introducing Ca10(PO4)6(OH)2 Potential Model. Phys. Chem. Chem. Phys. 2004, 6, 1860−1866. (83) 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.

6703


Molecular Dynamics Exploration of Ordered-to-Disordered Surface

Recommend Documents