Molecular Dynamics Exploration of Ordered-to-Disordered Surface

Subscriber access provided by UNIV OF DURHAM

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Molecular Dynamics Exploration of Ordered-to-Disordered Surface Structures of Biomimetic Hydroxyapatite Nanoparticles Qing Xie, Zhi-yu Xue, Haojie Gu, Cuifang Hu, Mingli Yang, Xin Wang, and Dingguo Xu J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b00178 • Publication Date (Web): 08 Mar 2018 Downloaded from http://pubs.acs.org on March 9, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Exploration of Ordered-to-Disordered Surface Structures of Biomimetic Hydroxyapatite Nanoparticles

Qing Xie1, Zhiyu Xue1, Haojie Gu1, Cuifang Hu1, Mingli Yang2,3, Xin Wang1* and Dingguo Xu1,2* 1

College of Chemistry, MOE Key Laboratory of Green Chemistry and Technology, Sichuan

University, Chengdu, Sichuan 610064, P. R. China 2

Genome Research Center of Biomaterial, Sichuan University, Chengdu, Sichuan 610064, P. R.

China 3

Institute of Atomic and Molecular Physics, MOE Key Laboratory of High Energy Density Physics

and Technology, Sichuan University, Chengdu, Sichuan, 610065, P. R. China

_________________ *

To

whom

correspondence

should

be

addressed:

[email protected]

[email protected] (D.X)

1

ACS Paragon Plus Environment

(X.W)

and

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

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 its environment. We, in this work, systematically explored the surface structures of HAP based on two different force fields, BornMayer-Huggins (BMH) parameters and INTERFACE force fields (IFF). 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 are 8.9 Å (BMH) and 10.3 Å (IFF), which are close to experimental measurement of 10-20 Å. Further studies on the surface models of (010) and (100) with different ion-rich terminations show that the disordered layers thickness of P-rich surface are systematically larger than Ca-rich surfaces. In the surface structure analyses, total pair distribution functions, static structure factor, 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 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.

2


Page 2 of 50

Page 3 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction It has been widely accepted that synthetic calcium phosphates containing ceramics can be used as bone grafts since 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 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 tissue inductive materials (TIM) is the one of the major development directions for the next-generation biomaterials.13-16 Therefore particular surface structure bearing tissue inductivity should be critical in future material design. The characterizations of surface structure of biomaterials are thus to be one of crucial issues in bone tissue regeneration and engineering. Currently, HAP particles are mainly obtained by various synthesizing method,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 composition of the HAP surface can be investigated and characterized by versatile instrumental methods, like nuclear magnetic resonance (NMR),18-20 ultrahigh-resolution transmission electron microscopy (UHR-TEM),21, 22 atomic force microscopy (AFM),23 X-ray photoelectron spectroscopy(XPS),24-26 time-of-flight secondary ion mass spectrometry(TOF-SIMS),24,

25

scanning Auger microscopy (SAM),27 scanning electron

3



microscopy (SEM),24, 27 X-ray diffraction (XRD),28, 29 and Fourier-transform infrared (FTIR).30, 31 These observations demonstrate that the great difference between the sintered HAP surface and bulk HAP crystal. Particularly, Jäger et al.

18

found that crystalline HAP core covered by a

disordered surface layer with about thickness of 1 nanometer, and the Ca/P ratio of the surface area is smaller than HAP core area. Bertinetti et al. 21 also investigated the surface features of nanosized HAP by UHR-TEM, and provided evidence that HAP nanoparticles are constituted by a crystalline core and a 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 inorganic material is still hindered by 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 like 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 MD simulations method by Bhowmik et al.37 By MD code, AlmoraBarrios et al.39 found that Hyp-Pro-Gly peptide in collagen prefer to interact with HAP (110) rather than 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

4


Page 4 of 50

Page 5 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with experimental proposal that HAP with amorphous surface layer has bioactivity53-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 of different interior and boundary regions as ordered/disordered microstructure is highly desired. We56 have proposed one strategy to construct such kind of model using simulated annealing 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 surfaces models for HAP. Four different surfaces, (001), (110), (010) and (100) were considered here. Our study could thus provide a possible model with 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 force field should be one of the key issues, since classical MD simulation based on force fields approach would be the only choice to tackle the large size system like HAP. Recently, a recent developed INTERFACE force fields (IFF)58 was reported, which is specifically designed to understand interactions between HAP surface and proteins. As suggested by previous work, Born-Mayer-Huggins (BMH) parameters are particularly suitable for investigation of issues related with bulk HAP model59. Therefore, it would be interesting to extend our simulation strategy using other force field, like 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 the understanding the interactions between HAP material and its biological environment.

5



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 MD simulated annealing 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 hydroxyapatite crystal structure. Basically, the hexagonal HAP crystal is in P63/m space group with ten Ca2+ ions, six PO43- and two OH- in the unit cell.65 Its detailed lattice parameters can be found in Table S1 of 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 group22. Therefore, it would be also interesting to further examine their properties upon simulated annealing 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, i.e., interior, middle, and the surface, which are plotted in Figure 1 and 2 for four different crystal surface directions. The interior part was fixed to mimic the real inner core, while the middle and surface parts were freely relaxed during the simulation to check their response under high temperature annealing simulation. In particular, the initial length for the surface layer was set to be about 20 Å, due to experimental 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 the different termination of HAP surface

6


Page 6 of 50

Page 7 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


could lead to different results. In fact, our strategy to build the slab model could result in either Carich or P-rich surface model. Hence, to test the effects of different termination pattern of HAP surface, we have then included models with both termination types in the simulation, but only for (010) and (100) surfaces. 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 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 annealed HAP surface structure. The trajectory relaxed at 300 K was saved for final analysis. In this work, atoms in HAP model were described using two force fields, i.e., BMH and IFF. Their corresponding potential energy forms are summarized in 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 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.

7



Once we have constructed an annealed HAP surface model, several structural analysis methods were applied to understand structural properties in this work, e.g., 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:

1

〈 〉 , 1 =

4 ∆

where is the number of the atoms in slab model, is the number density, and is the number of the 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:

1 〈 〉. 2 = 4 ∆

The static structure factor is one of 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. Static structure factor can be obtained by Fourier transformation of PDF. The partial static structure factor between the atom α and atom β is defined as: /

= &4 !" " #

%

∞

sin ! − 1# d, 3

in which is the Kronecker delta function, and cα = α /. The total static structure factor is then determined by: /

= !" " # ,

. 4

8


Page 8 of 50

Page 9 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In the characterization of micro-structure, the coordination number is used to describe the average atoms number in the first shell of the central atom, 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: .

- = 4 % d, 5

where R is the radius of the first shell of atom β around atom α. ISAACS program74 was used to compute the pair distribution function and the static structure factor S(q) for 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, e.g., crystal cell unit parameters,57 infrared spectra,58 melting point mean square displacement (MSD), elastic modulus,75 cleavage energies and immersion energies58 etc. As we pointed above, the surface of the sintered HAP appeared the disorder layer of 1-2 nm (10-20 Å), but the core part is nearly in its crystal form.18, 21, 22 Simulated annealing 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 HAP (100) plane based on BMH force field. Here, we will further examine 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 infrared (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 9



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, while nearly 3600 cm-1 in IFF case can be observed58. More importantly, relatively cleaner peak structure on this IR spectra obtained using 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 calculation and experimental measurements can be still established. Current IR results could ensure us that both BMH and IFF parameters might be suitable for further construction of layered HAP model in the annealing simulation. 3.2 Ordered-to-disordered HAP models We have in 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 2400K 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 synthesis of the HAP is about 1400~1600 K. Hence, we first set the target temperatures at 1400 K and 1600 K here in our simulated annealing (SA) computations. Similarly, only those surfaces with Ca-rich termination were considered in this section. Of course, such temperature 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.

10


Page 10 of 50

Page 11 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The initial and annealed structures of HAP (001), (010), (100) and (110) in two types force fields with Ca rich termination are displayed in Figure 1 and 2, and corresponding thickness values of the disordered layer for all plane directions are listed in Table 1. Indeed, although the annealing temperature (1400 K or 1600 K) is much lower than the melting point of HAP, obvious ordered-todisordered 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 1400 K 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 HAP slab model in BMH parameters is 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. But the thicknesses for those disordered layers are systematically larger than BMH based models. The difference becomes bigger with increasing of annealing starting 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 nonstoicheiometric terminations for HAP surfaces. In fact, theoretical simulation has suggested that stoichiometric surface like (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 nonstoichiometric terminations, according to different Ca/P ratio, we can have so-called Ca-rich and P-rich surfaces. Recently, Chiatti et al60 have investigated

11



the CO adsorption mechanism on HAP nanoparticle surface structure using 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 termination should be included here. 3.3 Termination and Temperature Dependence In this work, the effects of different terminations in the 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 K, 1600 K, 1800 K 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 surface 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 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 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 P-rich surfaces from Table 2. For example, in 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, respectively. Same situation is also for the models described using IFF parameters, but systematically more chaotic representation

12


Page 12 of 50

Page 13 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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, while 10.7 Å and 11.4 Å for (010)+P and (100)+P surfaces when the annealing temperature is 1400 K. Nearly the same trends can be found for the case obtained at 1600 K, although a little changes 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 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, while 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 BMH form than in 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 temperature (1400 K 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 last work built a so-called three-layered HAP models based on BMH force field parameters, but only Ca-rich HAP (100) surface at 2400 K was employed in our work56. A chaotic surface layer less than 10 Å was obtained. In this work, two annealing temperatures of 1800 K 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 of annealing temperatures. The anisotropic characteristics of HAP are also clearly shown by the different temperature dependence. Typically, for the models in IFF parameters, for (010)+Ca surface, there is only 2.8 Å increment from 1400 K to 2400 K, while as large as 7.1 Å increment for

13



(100)+Ca surface. More interestingly, for temperatures around the melting point (1923.16 K), e.g., 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 BMH form is calculated to be 8.9 Å, while 10.3 Å for the models in IFF parameters. Of course, different performance 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 for 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, since the IFF parameters are particularly designed for simulation of interface of 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 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 HAP model described using IFF might not be as stable as the case in BMH, esp., for the part close to surface. Therefore, in comparison with 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 like 2400 K, dramatic different performance 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

14


Page 14 of 50

Page 15 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


become completely chaotic and instable. However, models described in BMH form maintain very well except (100)+P. In other words, the thickness values of disordered layers for P-rich surface are systematically larger than Ca-rich surfaces, which might suggest that 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 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 thickness of the disordered layer, some common features can be observed for both parameters. For the HAP annealed structures, the locations of the ions (Ca2+, OHand PO43-) in the middle layers have tiny changes, while 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 has been reported in some literatures38, 40, 83. To illustrate motions of Ca2+ or PO43-, we then plotted the surface topology in Figure 6 and 7 from the top view. The phosphate groups seem to stay their original positions but undergo some rotations. Most of calcium ions, on the contrary, left their original positions and some cavities could come up. This could definitely lead to 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 unsaturated coordination status. However, during the annealing simulation, since OH- is much lighter than Ca2+ and PO43-, drastic movements can be observed. For the model subjected to annealing simulation, we could find that

15



some hydroxyl groups move from inside to the surface, and finally cover the surface calcium ions as shown in Figure 8 and 9. Such behavior could definitely decrease the surface energy and finally stabilize the surface zone, which is confirmed in our last work. In order to further illustrate the conformational changes of the annealed structures from interior to surface layer, we then calculated the total pair distribution functions (PDFs) for both HAP crystal and all annealed structures using 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 the 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 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 bears its crystal properties with atomic arrangement is pretty regular. In the medium-range area (r = 3~7 Å), obvious difference from interior part can thus be observed, which means significant atomic movements could occur during the simulation. Generally, the g(r) curves of 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 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.

16


Page 16 of 50

Page 17 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 X-ray diffraction experiments. Here, S(q) is used to tell us the structural difference among three layers of the HAP annealed structures expressing in different force fields. We calculated 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 were plotted in Figure 11 and 12 as well. All other surface directions have nearly the same trend, and are thus included in the supporting material (Figure 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 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 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 Figure 11 and 12 too. Clearly, curves only at short distance part (r < 3 Å) show nearly the same performance, which suggests atomic arrangement are relatively regular at 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 middle layer. Generally, much different static structure factors calculated here can clearly support the formation of a three-layers HAP model after

17



annealing simulation, and thus reflect the disordering difference between the middle and the surface layers. Once the ordered-to-disordered model via 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, since they can be easily obtained via integration of partial pair distribution function. Basically, the coordination number around the center atom could shed lights on the structural changes under simulation. Again, HAP (100) surface model at 1600 K was taken as the template system to address the changes of 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, while other atom pairs in middle and surface layers change from 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, O atoms calculated using IFF parameters, which were then illustrated in Figure 13. The average coordination numbers for atom pairs were summarized in Table 3. For comparison, coordination numbers obtained at 2400 K using BMH parameters are also included. Obviously, dramatic changes for numbers of calcium ion around 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, since this layer is fixed and pretty much the same as a HAP crystal structure. But in middle layers, due to increment of chaos, two types of coordination for P-Ca can be seen, i.e., 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 very stable bulk structure described in BMH parameters under simulated

18


Page 18 of 50

Page 19 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


annealing 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, for IFF case, about 17.4%, 7.9%, 7.9%, 66.8% of the phosphate groups are coordinated with six, seven, eight, and nine Ca atoms, which results in an average coordination number of 8.24 for P-Ca pair. Such kind of change of coordination number should be reasonable, since some drastic movements of ions at the ceramic surface must occur under high temperature and thus cause natural disordering for those atoms at 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 surface zone. Nearly identical coordination number between interior and middle layer 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 interface force field 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 HAP amorphous surface model, which could be more appropriate for further investigation of interactions between biomolecules and HAP material. 4. Conclusion 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

19



with biomolecules. Such treatment apparently ignore the fact of an amorphous surface layer (Ca/P ratio close to 1.5) identified experimentally. We have in last work56 proposed one strategy to reconstruct the surface structure using a simulated annealing MD simulation. The Hauptmann’s potential energy functional form was applied and a distinguished ordered-to-disordered structure along HAP (100) direction can be established. However, the Hauptmann’s parameters57 are suggested to be suitable for the bulk structure, and its applicability to HAP interfacial characteristics is thus questionable. Recently, Lin et al58 proposed another type of parameters that is specifically designed to simulate process occurring at material interfacial area. In this work, we then systematically explore the surface structures of hydroxyapatite (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 temperature 1400 K and 1600 K. Interestingly, HAP models based on IFF parameters have the same thickness layer sequence. But the thicknesses by IFF are systematically larger than those for BMH based models. Moreover, to understand effects of temperatures and terminations with nonstoicheiometric 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 of annealing temperatures. The anisotropic characteristics of HAP are also observed by the different temperature dependence. In addition our calculated results shown that the thickness values of disordered layers for P-rich surface are systematically larger than Ca-rich surfaces, which indicates that Ca-rich surface is more stable. Much different performances between two types of force fields have been found. The

20


Page 20 of 50

Page 21 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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 structures calculations. Interestingly, at 1800K, 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 theoretical predictions for the average disordered layer thickness values are 8.9 Å (BMH) and 10.3 Å (IFF), which are close to experimental measurement of 10-20 Å. Further surface structure analyses show that the ions (Ca2+, OH- and PO43-) in the middle layers have tiny changes, while about one-half of hydroxyl groups deflect about 90° in BMH force field and much less than 90° using IFF parameters. The total pair distribution functions were computed to further illustrate the conformational changes of the annealed structures from interior to surface layer. Moreover, our calculations reveal that the coordination numbers at the surface layer are much different from bulk structure, which is due to drastic motions of ions. According to these analyses, we suggest that IFF might be more suitable for constructing HAP amorphous surface model than BMH. In present molecular dynamics 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. Acknowledgment 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.

21



Supporting Information 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. This material is available free of charge via the internet at http://pubs.acs.org.

22


Page 22 of 50

Page 23 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) Dorozhkin, S. V.; Epple, M., Biological and Medical Significance of Calcium Phosphates. Angew. Chem.-Int. Edit. 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.; Gonzá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-Muñoz, E. M.; Velázquezcastillo, R.; Cabrera-Torres, J. L., Morphological Analysis of Hydroxyapatite Particels Obtained by Different Methods. Mater. Sci. Forum 2010, 681-686. (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. Design 2012, 34, 148-154. (9) Champion, E., Sintering of Calcium Phosphate Bioceramics. Acta Biomater. 2013, 9, 58555875. (10) Karnalanathan, 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.

23



(12) Cerroni, L.; Filocamo, R.; Fabbri, M.; Piconi, C.; Caropreso, S.; Condò, S. G., Growth of Osteoblast-like Cells on Porous Hydroxyapatite Ceramics: an in Vitro Study. Biomol. Eng. 2002, 19, 119-124. (13) Mourino, 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 Del. 2013, 10, 1353-1365. (14) Pérez, R. A.; Won, J.-E.; Knowles, J. C.; Kim, H.-W., Naturally and Synthetic Smart Composite Biomaterials for Tissue Regeneration. Adv. Drug Deliver 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) Ruhé, 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. Mat. Sci. Eng. C-Bio. S 2008, 28, 1381-1387. (18) Jä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.

24


Page 24 of 50

Page 25 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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. F.; Robinson, C., Physico-chemical Properties of Crystal Surfaces in Matrix-mineral Interactions During Mammalian Biomineralisation. Curr. Opin. Colloid In. 2002, 7, 124-132. (24) França, R.; Samani, T. D.; Bayade, G.; Yahia, L. H.; Sacher, E., Nanoscale Surface Characterization of Biphasic Calcium Phosphate, with Comparisons to Calcium Hydroxyapatite and β-Tricalcium Phosphate Bioceramics. J. Colloid Inter. 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.; Sjöberg, S., A Solubility and Surface Complexation Study of a Non-stoichiometric Hydroxyapatite. Geochim. Cosmochim. Acta 2009, 73, 257-267.

25



(27) Battistoni, C.; Casaletto, M. P.; Ingo, G. M.; Kaciulis, S.; Mattogno, G.; Pandolfi, L., Surface Characterization of Biocompatible Hydroxyapatite Boatings. Surf. Interface Anal. 2000, 29, 773781. (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, 25382543. (29) Tanaka, H.; Yasukawa, A.; Kandori, K.; Ishikawa, T., Surface Modification of Calcium Hydroxyapatite with Hexyl and Decyl Phosphates. Colloid Surf. A-Physicochem. Eng. Asp. 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 Interf. 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.; Delle Piane, M.; 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.; Stottt, M. J., First-principles Study of Hydroxyapatite Surfaces and Water Adsorption. Phys. Rev. B 2008, 78, 075427. (35) Zahn, D.; Hochrein, O., Computational Study of Interfaces between Hydroxyapatite and Water. PhyS. Chem. Chem. Phys. 2003, 5, 4004-4007.

26


Page 26 of 50

Page 27 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(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 Hydroxyapatitepolyacrylic 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 (011̅0) 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., Firstprinciples 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.

27



(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, 12841290. (47) Liao, C.; Xie, Y.; Zhou, J., Computer Simulations of Fibronectin Adsorption on Hydroxyapatite Surfaces. Rsc Adv. 2014, 4, 15759-15769. (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. L.; Leng, Y.; Chen, J. Y.; Zhang, Q. Y., A Comparative Study of Calcium Phosphate Formation on Bioceramics in Vitro and in Vivo. Biomaterials 2005, 26, 6477-6486. (51) Xin, R. L.; 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. Y.; Qu, S. X., TEM Study of Calcium Phosphate Precipitation on HA/TCP Ceramics. Biomaterials 2003, 24, 2125-2131. (53) Lukić, M. J.; Škapin, S. D.; Marković, S.; Uskoković, D., Processing Route to Fully Dense Nanostructured HAp Bioceramics: From Powder Synthesis to Sintering. J. Am. Ceram. Soc 2012, 95, 3394-3402. (54) Veljović, D.; Jokić, B.; Jankovic-Castvan, I.; Smiciklas, I.; Petrović, R.; Janaćković, D., Sintering Behaviour of Nanosized HAP Powder. Key Eng. Mater. 2007, 330-332, 259-262.

28


Page 28 of 50

Page 29 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(55) Veljović, D.; Jokić, B.; Petrović, R.; Palcevskis, E.; Dindune, A.; Mihailescu, I. N.; Janaćković, 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.

29



(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} P-Rich/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) Nosé, S., A Unified Formulation of the Constant Temperature Molecular-Dynamics Methods. J. Chem. Phy. 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., Plasma-Sprayed Hydroxyapatite Plus Titania Composite Bond Coat for 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 Spectrain Molecular-Dynamics Simulations of Anharmonic Solids and Surfaces. Surf. Sci. 1992, 262, 189-207.

30


Page 30 of 50

Page 31 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(72) Duan, G.; Xu, D. H.; Zhang, Q.; Zhang, G. Y.; 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 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 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. M.; Zhang, X. D.; Chen, J. Y.; Zeng, S. X.; Degroot, 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.

31



(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, 19781981.

32


Page 32 of 50

Page 33 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Crystal structure and annealed structures at 1400K based on different force fields for HAP Ca-rich surfaces, and (a), (b), (c), (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. 33



34


Page 34 of 50

Page 35 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

35



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

36


Page 36 of 50

Page 37 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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.

37



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

38


Page 38 of 50

Page 39 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

39



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

40


Page 40 of 50

Page 41 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

41



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

42


Page 42 of 50

Page 43 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 10. Total pair distribution functions of the crystal and annealed structures for the HAP Carich models at 1600K in BMH and IFF types force fields, (a) = HAP (010), (b) = HAP (100).

43



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], P-O [S(q)P-O] of the three layers of the HAP (100) Ca-rich structure annealed at 1600K in IFF force field, (a) = S(q)total, (b) = S(q)O-Ca, (c) = S(q)P-Ca, (d) = S(q)P-O.

44


Page 44 of 50

Page 45 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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

45



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

46


Page 46 of 50

Page 47 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1. Calculated thickness values (Å) of disordered surface layers for HAP surface models.

Temp/K 1400 1600

(001) 4.4 5

BMH (010) (100) 5.7 9.9 5.8 10.3

(110) 10.1 11.3

(001) 5.9 6.7

47


IFF (010) 7 7.5

(100) 10.7 11.8

(110) 12 15.2


Table 2. Calculated thickness values (Å) of disordered surface layers for HAP surface models. Temp/K 1400 1600 1800 2400

BMH (010)+Ca (100)+Ca (010)+P (100)+P 5.1 5.7 9.9 10.1 5.4 5.8 10.3 10.3 5.9 6 11.7 11.9 6.6 7.1 12.5 24.7

IFF (010)+Ca (100)+Ca (010)+P (100)+P 5.7 7 10.7 11.4 6.1 7.5 11.8 12 6.7 8.2 12.6 13.5 8.5 14.1 44.2 88

48


Page 48 of 50

Page 49 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3 The average coordination number of three layers for HAP (100) annealed in BMH at 2400K[Ref.56] and IFF at 1600K parameters.a

A-B pair P-Ca Ca-P O-Ca Ca-O P-O H-O a The pair of

BMH[Ref.56] Interior

Middle

9.00 5.40 3.00 7.80 4.00 1.00 atoms A−B

9.00 5.40 2.97 7.72 4.00 1.00 represents

IFF

Surface

Interior

8.24 9.00 4.12 5.40 2.51 3.00 5.51 7.80 4.00 4.00 1.00 1.00 the average number of B

Middle 8.95 5.35 2.80 7.68 4.00 1.00 around A.

Surface 8.24 4.94 2.63 6.98 4.00 1.00 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.

49



TOC Graphic

50


Page 50 of 50

Molecular Dynamics Exploration of Ordered-to-Disordered Surface

Recommend Documents