Surface Energetics of the Hydroxyapatite Nanocrystal–Water Interface

Oct 14, 2014 - *E-mail: [email protected]. .... Method Evaluations for Adsorption Free Energy Calculations at the Solid/Water Interface through Metadyn...
0 downloads 0 Views 8MB Size
Download PDF
Article pubs.acs.org/Langmuir

Surface Energetics of the Hydroxyapatite Nanocrystal−Water Interface: A Molecular Dynamics Study Weilong Zhao,† Zhijun Xu,† Yang Yang,‡ and Nita Sahai*,† †

Department of Polymer Science, University of Akron, 170 University Ave., Akron, Ohio 44325-3909, United States Department of Chemistry and Biochemistry, Rowan University, 201 Mullica Hill Rd., Glassboro, New Jersey 08028, United States

‡

S Supporting Information *

ABSTRACT: Face-specific interfacial energies and structures of water at ionic crystal surfaces play a dominant role in a wide range of biological, environmental, technological, and industrial processes. Nanosized, plate-shaped crystals of calcium phosphate (CaP) with nonideal stoichiometry of hydroxyapatite (HAP, ideal stoichiometry Ca10(PO4)6(OH)2) comprise the inorganic component of bone and dentin. The crystal shape and size contribute significantly to these tissues’ biomechanical properties. Plate-shaped HAP can be grown in the presence of biomolecules, whereas inorganically grown HAP crystals have a needlelike shape. Crystal morphology reflects the relative surface areas of the faces and, for an ideal inorganically grown crystal, should be governed by the surface energies of the faces with water. Interfacial energies and dynamics also affect biomolecule adsorption. Obtaining face-specific surface energies remains experimentally challenging because of the difficulty in growing large HAP single crystals. Here we employed molecular dynamics (MD) simulations to determine nanocrystalline HAP−water interfacial energies. The (100) face was found to be the most favorable energetically, and (110) and (004) were less hydrophilic. The trend in increasing interfacial energy was accompanied by a decrease in the average coordination number of water oxygen to surface calcium ions. The atomic-level geometry of the faces influenced interfacial energy by limiting lateral diffusion of water and by interrupting the hydrogen bond network. Such unfavorable interactions were limited on (100) compared to the other faces. These results provide a thermodynamic basis for the empirically observed trends in relative surface areas of HAP faces. The penetration of charged biomolecules through the interfacial water to form direct interactions with HAP faces, thus potentially influencing morphology, can also be rationalized.

1. INTRODUCTION The physicochemical properties of the mineral−water interface play a critical role in a diversity of biological, geochemical, technological, and industrial processes.1,2 Interfacial energies and dynamics of water, inorganic ions, and organic molecules at mineral surfaces influence phenomena such as the nucleation of the mineral crystal directly from solution or of a precursor metastable phase which then transforms to the stable mineral, crystal growth, biomineralization, self-assembly of natural and synthetic organic−mineral composite materials, adsorption and surface precipitation of toxic metals from the environment, geomicrobial remediation of environmental contaminants, and so forth. Understanding these fundamental principles also provides a basis for the intelligent of design of bioinspired materials for orthopedic and dental applications. In the present study, we focus on determining the interfacial energies and dynamics of water at different faces of hydroxyapatite (HAP). In calcified connective tissue, such as bone and dentin, nonstoichiometric HAP nanocrystals of platelike morphology occur in specific structural registry with the extracellular matrix protein, collagen.3,4 Thus, these tissues form a natural hierarchical, composite material endowed with unique mechan© XXXX American Chemical Society

ical and biochemical properties compared to their individual components alone.4−9 For example, bone has higher hardness than proteins and higher fracture resistance than inorganic calcium phosphate (CaP) phases. HAP crystals found in human bone exhibit an aspect ratio of 5−7, with the largest face being (300).4 It has been proposed that HAP morphology modulation is a result of the face-specific adsorption of acidic noncollagenous proteins (ANCPs) and small molecules, such as bone sialoprotein and citrate, respectively.10 The in vivo adsorption is believed to be mediated through acidic carboxylate groups (e.g., glutamate) of the ANCPs and citrate, and this mechanism has also been shown in vitro.11−17 Recent efforts by computer simulation have enabled the observation of such interactions with atomic-level resolution, by using stoichiometric HAP as model surface.18−20 In principle, the stronger binding of acidic biomolecules to the (100) face than the other faces could inhibit crystal growth in the [100] direction, resulting in a larger surface area of (100) compared Received: August 7, 2014 Revised: October 9, 2014

A

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

surface water coverage, but the calculation was limited to only a few ( (110) > (001).4 In some cases, the (001) face is suppressed entirely by other faces. This empirical observation would suggest that the (100) is energetically the most favorable face in contact with water. The interfacial energy of HAP micrometer-sized bulk powders obtained by precipitation from aqueous solution has been measured using the thin-layer wicking technique,31 but face-specific information is missing. Experimentally determined face-specific interfacial energies are lacking, because of the difficulty in growing large HAP single crystal with ideal faces to measure contact angles. A similar problem limits the determination of the structure, dynamics, and properties of water at HAP surfaces by experimental techniques such as synchrotron X-ray reflectivity. The latter method has been applied, however, to the (100) face of fluorapatite (FAP, Ca5(PO4)3F), which is much easier to grow as large ideal single crystals.32 The situation is even more challenging for nanometer-sized HAP crystals. In the present study, we approach this fundamental issue by applying molecular dynamics (MD) simulation to calculate face-specific HAP− water interfacial energies of ideal HAP nanocrystals. Previous attempts have been made using quantum mechanical density functional theory (QM-DFT) or molecular mechanical MD simulation to attain surface energies and small molecule binding energies on ideal, stoichiometric faces of HAP, but these studies have limitations in various aspects. DFT has been employed to calculate the surface energies of HAP faces in vacuum while surface hydration was not included.33,34 Astala and Stott35 have introduced DFT scheme to study the surface energy change of different faces of HAP as a function of

2. METHODS 2.1. Preparation of Crystal Surface. HAP bulk crystals were prepared using X-ray crystallography structure reported by Wilson et al.48 OH− groups are aligned in the [001] direction within channels constituted by calcium ions (Ca2+). For hexagonal HAP, the P63/m B

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Figure 1. Structures of HAP surfaces used in the present study. Black dashed lines shows the position of cleaving bulk HAP crystal, and the numbers in parentheses indicate the surface designation. Atoms of HAP crystal are represented by ball-and-stick model. Color designation: cyan, calcium; brown, phosphorus; red, oxygen; white, hydrogen. space group renders OH− groups in adjacent parallel channels to possess the same orientation, whereas within the same channel adjacent OH− groups take alternating orientation through a mirror plane of symmetry.49 Slabs with a particular crystallographic zone orientation were obtained from HAP bulk crystals by cleaving the crystal at the plane perpendicular to the respective crystallographic direction. As with all previously published simulation studies for HAP, ideal surfaces were used in the MD simulations without any chemical modifications, such as protonation and hydroxylation of phosphate (PO43−) and Ca2+, or consideration of topological features, such as steps, terraces, and kinks. The reason is that this type of detailed information is experimentally lacking or limited for HAP, and hence, initial structures cannot be constructed.50 Four surface terminations were selected for the current work (Figure 1). Indices in square brackets refer to the zone axis. For example, [100] refers to the direction parallel to the a-axis. The Miller index of each face denotes the relative distance of that face along the unit cell axis where the termination is chosen. Thus, (100) means that the surface is cleaved at one unit cell distance along the a-axis of the crystal, and (200) indicates a termination at half the distance of the unit cell in the same direction. 2.2. Molecular Dynamics Simulation. 2.2.1. Force Field. The potential energy parameters for bulk HAP were developed by Hauptmann et al. using the rigid bond and angle model.45 The intermolecular interaction in the original work was described using the Born−Mayer−Huggins potential. In our previous study, the BMH potential was modified to fit the 12-6 Lennard−Jones potential, and the refitted force field reproduces HAP bulk properties with high accuracy.18 The TIP3P model was employed to treat water molecules.52 For interactions at the HAP−water interface, the Lorentz−Berthelot mixing rule was used to obtain cross-terms for interactions between nonbonded atoms. Point charges for summation of electrostatic potential were taken from Hauptmann et al.’s force field and were employed by recent work of RosettaSurface docking application.38 Benchmarking of MD simulation using this force field against experimental data and DFT calculations have proven the accuracy of the current scheme in describing the HAP−water and HAP−water−biomolecule interface.18,32,50 2.2.2. Simulation Details. Simulations were performed using GROMACS 4.5.5 package.51 HAP slabs were cleaved from bulk crystal. The slab dimensions in x- and y-directions were chosen to be 6 × 6 times of the unit cell lengths along the a- and b-axes of hexagonal HAP. The coordinates were transformed to orient the direction perpendicular to the target surface along the z-direction of the simulation box. Depending on the different faces studied, the slab thickness in the z-direction ranged from 1.62 to 3.44 nm. The HAP

crystal with infinite periodic images in x-, y-, and z-directions was relaxed by carrying out 100 ps NVT ensemble simulations at a temperature of 300 K. The z-dimension of the simulation box was then increased by 100 Å, creating a vacuum gap between two identical surfaces in order to eliminate interactions with the adjacent simulation cell. This simulation box was used subsequently for the solid−vacuum and the solid−water interfacial energy calculations. To minimize the surface energies of the slabs, a 10 ns NVT simulation was performed for each surface. The solid−vacuum interfacial energy per unit surface area, γSV, was then calculated as

γSV =

Eslab − E bulk,HAP (1)

2A

where Eslab is the potential energy of the slab in contact with the vacuum gap, Ebulk,HAP is the potential energy of the bulk HAP crystal normalized to the same number of molecules as in the structure used for calculating Eslab, and A is the surface area of the crystal face. All energy components were averaged over the last 5 ns, and the error was estimated with the block average method. An independent 10 ns NPT ensemble simulation was performed for hexagonal HAP bulk crystal in order to calculate Ebulk,HAP. To calculate solid−water interfacial energy, the vacuum gap on top of each slab was filled with water molecules by using the GROMACS build-in utility. A water layer of 5 nm thickness ensures that sufficient molecules are included to represent accurately the bulk water phase beyond the interfacial region. The solid−water interface energy per unit surface area, γSL, was calculated according to the following equation:

γSL =

ESL − (E bulk,HAP + E bulk,water) 2A

(2)

where ESL is the potential energy of the solid−water system and Ebulk,water is the potential energy of bulk water normalized to the same number of molecules as in the solid−water system. The value for ESL was averaged for the last 5 ns of the simulation. Ebulk,water was calculated as the average of a 10 ns trajectory for a pure water system with simulation box size of 40 Å × 40 Å × 40 Å. The potential energy obtained for individual water molecule (−41.1 kJ·mol−1) is consistent with the experimental value (−41.5 kJ·mol−1).52 The solid−water system was equilibrated initially for 100 ps using the NVT ensemble and then simulated for another 500 ps with the NPT ensemble. For NPT equilibration, semi-isotropic pressure coupling18 was applied to ensure proper relaxation of the pressure tensor in the z-direction of the box. Finally, a 15 ns NPT MD simulation was performed and the last 5 ns trajectories were used for data production of interfacial water C

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Figure 2. Atomistic density profile of water oxygen (Ow) and hydrogen (Hw) along the z-axis of the simulation box for HAP faces: (a) (100), (b) (110), (c) (001), and (d) (004). The position of the HAP−water interface along the z-axis is determined by the thickness of different HAP slabs. The values for Ow are scaled by a factor of 2 for clarity. Gray dashed lines indicate the density corresponding to that of bulk water (1 g·cm−3).

Figure 3. Water adsorption structures for the first two hydration layers on four HAP surfaces: (a) (100), (b) (110), (c) (001), and (d) (004). Water molecules in the first and second layers are shown by yellow and purple licorice representation, respectively. Only the first hydration layer is displayed on the (110) face. Water molecules in the bulk aqueous solution are not shown for clarity. Color designation for atoms in HAP is the same as in Figure 1. properties. Nosé−Hoover and Parrinello−Rahman methods were selected, respectively, for temperature and pressure coupling schemes.53−55 The particle mesh Ewald (PME) method is used to treat long-range electrostatic interaction.56 A 14 Å cutoff was applied for dividing real and reciprocal space. A switch radius of 12 Å was used for the Lennard−Jones potential. Initial velocity of particles in the system was assigned randomly on the basis of Maxwell distribution. 2.3. Analysis of Simulation Results. The results of MD simulation were analyzed based on the last 5 ns MD trajectories, using standard methods for interfacial water density normal to the surface, planar water density parallel to the surface, and radial distribution function (RDF). In addition, the orientation of a single water molecule at the surface was determined by calculating order parameter S10 = cos θ, where θ is the angle between the direction normal to the surface and the water dipole moment (μ⃗ ).57 A value of 1.0 indicates that the plane of the water molecule is oriented perpendicular to the surface, and a value of 0 indicates that the water dipole is at 90° with respect to the surface normal (i.e., parallel to the surface). Another relevant parameter is cos β, where β is the angle

between the normal to the surface and the vector normal to the plane of the water molecule (n⃗). A plot of cos β versus cos θ for water molecules within the first and second hydration layers provides a probability distribution of water orientation. The point at which cos2 θ + cos2 β = 1 represents that both n⃗ and μ⃗ are in the same plane with the normal to the surface.

3. RESULTS 3.1. HAP−Water Interface Characterization. Results of HAP−water interface characterization are presented below in terms of density, snapshots of the structure at the interfacial region and RDF (Figures 2−4). Interfacial Water Density Normal to Surface. The number density profiles of water oxygen (Ow) and water hydrogen (Hw) atoms at the (100), (110), (001), and (004) faces are shown in Figure 2. High density peaks (>100 atoms.nm−3) and fluctuations in density are seen for all the surfaces except D

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

(110). Two surface hydration layers can be recognized for all the surfaces, but the (110) face is again the exception (Figures 2 and 3). Taking the z-coordinate of the topmost Ca2+ ions as the reference, the first and second density peaks of Ow on the (100) face appear at 1.2 and 2.4 Å, respectively; the corresponding peaks of Hw appear at 0.6 and 2.4 Å (Table 1). Furthermore, a close examination of the density profile Table 1. Hydration Layer Thickness and Position of Density Profile Peaks Shown in Figure 2a density peak position (Å) hydration layer thickness (Å)

surface (100) (110) (001) (004)

layer layer layer layer layer layer layer layer

1 2 1 2 1 2 1 2

peak density (atoms· nm−3)

Ow

Hw

Ow

Hw

3.6 2.4 6.4

1.2 2.4

0.6 2.4

140 85 80

87 90 75

3.8 2.3 3.9 2.4

0.7 3.6 2.2 4.4

0.5 3.4 1.5

100 84 105 80

88 80 97

Figure 4. (a) Radial distribution function (RDF) of Ow−Ow pair for the first hydration water layers on the (100), (001), and (004) faces. (b) RDF of Ow−Ow pair for bulk water. The gray box highlights the change in shape of the primary peak on the three faces compared to that of bulk water. The black arrow indicates that the height of the secondary peak follows the order of (100) < (001) < (004).

first RDF peaks shifts slightly from 2.9 Å of pure water, to 3.1 Å of the first hydration layers on all the three faces (Figure 4, gray box). This phenomenon suggests that the normal water packing structure is disrupted to some extent and the hydrogen bond network is partially broken. For comparison with another wellstudied biomineral, the disrupting effect of HAP surfaces to the structure of adsorbed water is less pronounced than the effect of the calcite (104) surface.57 Beyond the first hydration layer the RDF for the three faces exhibit distinct structural features, suggesting different types of water adsorption behavior (Figures 3 and 4). In particular, the secondary peak is more pronounced at the (001) and (004) faces compared to the (100) face (Figure 4, black arrow). We propose that this result is the evidence of difference water adsorption behavior, as also seen in subsection 3.3. 3.2. Surface Energetics. For the convenience of understanding the critical role of solid−water interfacial energy (γSL) in determining HAP crystal morphology, we provide here the definition of Gibbs free energy change as a function of the nucleus radius: 4 ΔG = πR3ΔG bm + 4πR2γSL (3) 3

a

The height of the density peaks above each face is determined using the z-coordinate of the topmost Ca2+ ion as the reference surface. The z-coordinate of the Ca2+ ions is reported as the average of the last 5 ns of MD production stage.

reveals a bimodal distribution for Ow in the first hydration layer, which suggests that water molecules adopt both oxygen-down and oxygen-up conformation at the HAP surface (Figure 3). This result is attributed to the electrostatic interactions of surface Ca2+ and PO43− ions with Ow and Hw, respectively. The density profiles also reveal that water adsorbed on the (100) face has the highest density compared to the other three faces (Figure 2). In detail, the density is about 2 times higher than that of pure bulk water, suggesting the strong hydrophilic nature of the (100) face. In contrast the density peaks obtained for Ow in the first hydration layers of the (001) and (004) faces are unimodal and the position of the Ow peak overlaps with that of Hw. Taken together, these results suggest randomly distributed orientation of the water dipole on the (001) and (004) faces. Water exhibits only a single adsorption layer on the (110) face and has a low density peak, which reflects unfavorable interactions between this surface and water. Furthermore, the peak is not as well-defined as those on other three faces. Examining the trajectory of the (110) surface atoms during the MD simulation shows that the crystal structure at the surface is disturbed, presumably as a result of interaction with water (Figure 3). Because of the difficulty in finding a well-defined hydration layer on the (110) face, we exclude this surface from subsequent analysis of water adsorption structure. Interfacial Water RDF. The calculated RDFs for the Ow−Ow pair in the first surface hydration layer of (100), (001), and (004) faces are shown in Figure 4. Because we are dealing with a two-dimensional situation with respect to the geometry of the water layer, the element volume for integrating number density is considered to be a cylindrical shell instead of a spherical shell. The reference density to calculate the RDF is the density of the hydration layer itself rather than that of bulk water. The interfacial RDFs are compared to that of the Ow−Ow pair in a pure bulk water system. It is evident that the position of the

where ΔGm b . is the Gibbs free energy of the bulk crystal per unit size and R is the radius of the nucleus. The calculated values of interfacial energies, γSL and γSV (eq 1 and 2), for several major faces observed in inorganically grown HAP crystals are reported in Table 2. The γSV value obtained for the (001) face is 1.29 ± 0.07 J·m−2, which is in good Table 2. Calculated γSV and γSL Values for Each Surface of HAPa surface

slab thickness (nm)

(100)_Ca (110)_Ca−OH (001)_Ca (004)_Ca−PO4

1.62 2.86 3.44 3.44

γSV (J·m−2) 1.93 2.38 1.29 2.10

± ± ± ±

0.06 0.07 0.07 0.07

γSL (mJ·m−2) 324.8 698.6 −58.9 1046.2

± ± ± ±

4.2 33.1 9.7 52.4

a

On the nomenclature employed, atoms which appear after the face name indicate their presence at the top of the surface.

E

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

agreement with DFT calculations (1.34 J·m−2).33 This result provides additional support for the accuracy of our force field. It is also evident that the calculated interfacial energies are highly dependent on the specific face (Table 2). For example, γSV is 1.93 J·m−2 for (100), 2.38 J·m−2 for (110), and 1.29 J·m−2 for (001). The (100) face has lower γSL value (324.8 mJ·m−2) than the values of the (110) and (004) faces (698.6 and 1046.2 mJ·m−2, respectively). This result is consistent with previous experimental reports, in which high-resolution transmission electron microscopy studies revealed that the major face of inorganically grown HAP is (100), and the (110) and (004) faces have very small surface areas.4 The (110) and (004) faces are similar in having high γSL, even though their γSV values are very different. The high γSL of these two faces provides an unfavorable contribution to the overall free energy of crystal growth (eq 3). As any face enlarges by incorporating free ions from solution into the crystal structure, the enthalpy of the bulk crystal decreases and the entropy increases. However, the high γSL is an unfavorable contribution that leads to an increase of ΔG, thus preventing the further enlargement of these faces. In comparison, γSL of the (100) face is smaller than that of others faces, so increasing its surface area does not lead to a significant increase in ΔG. Accordingly, the (100) face becomes the largest crystal face. A negative γSL value (−58.9 mJ·m−2) was obtained for the (001) face. This appears nonintuitive, but a negative surface energy between CaP phase and water has been reported previously in the literature, and it was interpreted as an indication of thermodynamic instability of the interface.58 A similar result for the (110) face was obtained for θ-Al2O3, which was related to the dissociative adsorption of water.59 It is important to note that, in both cases, the interfacial energy not only reflects the reversible work required for the creation of new surface areas, but also the intermolecular intactions between the solid and the adsorbed water, which could lead to spontaneous surface relaxation. In particular, the interaction between Ca2+ and Ow is responsible for the negative enthalpy contribution (−ΔH), resulting in a negative interfacial energy. Furthermore, empirical observation supports out results in that the (001) face is rarely observed in both inorganic and biogenic HAP crystals, Therefore, we conclude that the negative γSL value obtained in the present study for the (001) face indicates the thermodynamically instability of this face when brought in contact with aqueous environment. 3.3. Correlation between Water Adsorption Behavior and Surface Energy. Ca2+−Ow RDF. For an ionic solid and a polar solvent, such as HAP and water respectively, it may be expected that the atomic-level structure of the various faces is highly related to the adsorption behavior of the solvent and to the calculated interfacial energies. This is because the electrostatic force between Ca2+ and Ow largely influences the overall potential energy, rather than other kinds of surface water interactions such as van der Waals force. Therefore, the possible effect of HAP surface structure on water adsorption behavior and γSL is investigated. We first calculate Ca2+−Ow RDF, as well as the coordination number of Ow to Ca2+ for the (100), (001), and (004) faces (Figure 5). The average coordination number of Ow to the topmost Ca2+ layer is calculated by selecting Ow within a distance of 2.7 Å and by normalizing to the number of Ca2+ present at the surface. The RDF of Ca2+−Ow pair of the (001) face (Figure 5a) shows the highest value for the first peak compared to the (100) and (004) faces. Furthermore, the largest coordination number of

Figure 5. RDFs of the Ca2+−Ow pair (a) and coordination numbers of Ow to Ca2+ in the topmost atomic layer of the (100), (001), and (004) faces. Results shown are averaged for the last 5 ns of the MD production stage.

Ow to Ca2+ is obtained for the (001) face (Figure 5b). Both results indicate that the greatest number of Ca 2+ −O w intermolecular bonds are formed at the (001) surface. The (004) face shows the lowest first peak of Ca2+−Ow RDF and the smallest coordination number of Ow to Ca2+. These results suggest that water adsorption to the (004) face is less favorable than on the (100) face, which is in accordance with the solid− water interfacial energy trend (Figure 6).

Figure 6. Correlation between coordination numbers of Ow to Ca2+ and γSL values for the (100), (001), and (004) faces. Error bar indicates standard error of γSL. Linear fitting equation: Y = 1914− 540X; R2 = 0.95.

The results shown in Figures 5 and 6 suggest that water adsorption behavior is governed by the coordination of Ow to Ca2+ on each surface. Strong adsorption is correlated with high coordination number of water molecules to surface Ca2+, which is expected when the surface is highly hydrophilic and has a low value of γSL. However, when the coordination number of Ow to Ca2+ becomes too high, as in the case of the (001) face, a negative γSL value is obtained. This correlation between water adsorption behavior and surface energy is well illustrated by the linear fitting in Figure 6. Planar Water Density. The interaction of adsorbed water to HAP surface can be visualized by calculating the twodimensional density distribution of the first hydration layers on the (100), (001), and (004) faces in the x−y plane (Figure 7a−c). The corresponding structure of the topmost atom layer of HAP and water molecules of the first hydration layer are shown in Figure 7d−f. Note that, despite the high density of Ca2+ cations and PO43− anions on the (100) face, the variations in planar density of water across the surface are not large and F

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

the lateral correlation with atom lattice positions is not obvious. The main feature of the planar density profile is that the high density regions correspond to the lattice “void” created by three PO43−, which is a result of strong hydrogen bond interaction between Hw and oxygen of PO43−. In contrast, large local variations are observed for the water layer on the (001) face and the low water density region overlaps with the position of Ca2+ on the (001) surface. Thus, the emergence of the low water density region is a result of volume occupation by Ca2+. Areas of high density are seen at the perimeter of each low density location, indicating the position of water molecules coordinated by Ca2+. Comparing the crystal structures of the (100) and (001) faces, we find that similar vertical positions of PO43− and Ca2+ on the (100) face provide a smoother surface compared to the (001) (Figure 3a and c). The lower atomic roughness of the (100) face results in low coordination number of Ow to Ca2+. The (001) face is rougher and the PO43− are not fully coordinated with Ca2+. This leaves the top layer of Ca2+ to be highly coordinated by Ow (Figure 5). We infer that the rough structure of the (001) face is the origin of the negative interfacial energy. As a consequence, the crystal surface undergoes remodeling by water molecules to eliminate local hydrophobic regions, and the surface Ca2+ provides potential reaction sites for adsorbed water. Driven by such reaction potential, the (001) face would grow primarily by nonequilibrium process, such as establishing new steps and kinks on the surface. These atomic-level details of surface lattice position result in the small surface area of (001) facets seen in inorganically grown HAP crystals. The adsorbed water on (004) face exhibits the highest degree of planar density variation among the three surfaces. Small localized regions show density as high as 8 water molecules/ nm3, while the water density approaches zero on large areas of the (004) surface (Figure 7c). This observation indicates that the normal lateral diffusion of water is suppressed at the interface, which contradicts our intuitive expectation of full water coverage on a charged surface. However, this result helps to explain the high γSL value calculated for (004) and suggests unfavorable interactions between this face and water. Interfacial Water Order Parameter. The preferential orientation of water molecules on HAP surfaces affects the interfacial energy. Therefore, we extend our study of HAP− water interface by calculating orientation order parameter of

Figure 7. Planar water density profile of the first hydration layer (a−c) and the corresponding adsorbed water conformation (d−f) for (100), (001), and (004) faces. For water density profile, blue color indicates high density (8 water molecules per nm3) and red color indicates low density (0 water molecule per nm3). The corresponding structure of the topmost atom layer of HAP is shown in van der Waals sphere and water molecules of the first hydration layer are shown in yellow licorice. Color designation: cyan, calcium; brown, phosphorus; red, oxygen; pink, hydrogen; yellow, water. For clarity, only the topmost layer of atoms on each surface are shown in color and the rest are represented by “ghost” atoms.

Figure 8. (a) Schematic illustration of the definition of the angles θ and β related to the water orientation order parameter. Z is the direction normal to the crystal surface. The white box represents the molecular plane of water. θ is the angle between Z and the vector normal to the molecular plane (n⃗). β is the angle between Z and the dipole moment of water (μ⃗ ). (b−e) Water orientation order parameter of the first two hydration layers on (100) and (001) faces of HAP. Each dot represents an individual water molecule. The regions with high dot density indicate high probability. G

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

4. DISCUSSION The growth of a crystal involves numerous kinetic steps, which may be influenced by several factors, such as interfacial water behavior, incorporation of nonstoichiometric ions, and interactions with organic molecules. For an ideal crystal growing at thermodynamic equilibrium, the relative interfacial energies of the faces will govern their relative surface areas and hence the crystal morphology. We employed all-atom equilibrium MD simulations to study the surface energies, water structure, and diffusion behavior of water at different faces of a single nanocrystal of HAP. These results are summarized below for each face in order to facilitate the discussion on morphology of inorganically grown crystals and biomolecule interactions with HAP faces in biomineralization. The (100) face was found be the most stable, having the smallest positive solid−water interfacial energy and the hydrogen bond network was the least disturbed in the first hydration layer compared to the (001) and (004) faces. The coordination number of Ow to surface Ca2+ was smaller than that of the (001) face, but larger than for the (004) face. This trend was correlated with the relative stability of the various faces. The planar density profile of the first hydration layer on the (100) face suggested a homogeneous distribution of adsorbed water over the entire face. The structure order parameter distribution indicated that water molecules primarily take “oxygen up” or “oxygen down” conformations on (100) because of attractive forces with PO43− or Ca2+. The results of interfacial energy and water adsorption structures are consistent with the empirical observation of the relatively large surface area of the (100) face in inorganically grown HAP crystal. A negative interfacial energy was obtained for the (001) face. The coordination number of Ow to surface Ca2+ was the highest for this face, suggesting that strong electrostatic interactions were responsible for the negative energy. This inference was confirmed by the binding structure obtained in which a single water molecule interacts simultaneously with both PO43− and Ca2+. The planar density profile showed that the strong attraction of water molecules limited the lateral diffusion of water. Thus, the (001) face is thermodynamically unstable and subject to surface relaxation in aqueous environment, which explains why this face is rarely observed in HAP crystals grown inorganically. The highest positive interfacial energy and the lowest coordination number of Ow to surface Ca2+ were obtained for the (004) face. The lateral diffusion of interfacial water was highly limited as indicated by the planar density profile. Thus, the (004) face is less hydrophilic than the (100), resulting in a relatively small surface area. The unique crystal shapes of many biominerals compared to their inorganic counterparts has been attributed in the literature to the effect of charged biomolecules, which may interact preferentially with specific crystal faces, thus affecting the growth rate in those directions. Previous computational studies have pointed out that strong binding of charged biomolecules requires penetration into the interfacial water layer. Using MD simulation, it was shown that a Ti-binding peptide possesses high binding affinity as a result of incorporating the charged side chains into the interfacial water layer upon peptide adsorption.60 Similarly, the common feature of the HAPbinding peptides is the presence of charged side chains of amino acid residues, which form intermolecular bonds with Ca2+ and PO43− in the presence of interfacial water.

water molecules in the hydration layers as described in Methods (Figure 8a). We choose two representative surfaces, (100) and (001). Water orientation on the (001) and (004) faces are very similar, so the latter is not presented here. In the regions with a high density of dots, different kinds of binding modes of water can be identified for (100) and (001) faces (Figure 9). For the first layer on the (100) face, water dipole

Figure 9. Different binding modes of water molecule on HAP surface: (a) θ = 0° on (100) face; (b) θ = 180° on (100) face; (c) θ = 0° on (001) face; and (d) θ = 90° on (001) face.

orientation has two major distributions at cos θ = 1 and −1 (Figure 8b). These two values correspond to two configurations in which two Hw form hydrogen bonds with oxygen of two PO34− (Figure 9a), and Ow forms intermolecular bonds with two Ca2+ (Figure 9b), respectively. The observation of these two preferential orientations is consistent with the bimodal atomistic density profile shown in Figure 3a. In contrast, the second layer on the (100) face does not show preferential orientation of water dipole (Figure 8c). On the (001) face, two configurations in which cos θ = 1 and cos θ varies from −0.5 to 0, are identified for the first hydration layer (Figure 8d). The value at cos θ = 1 indicates the formation of hydrogen bond between Hw and the oxygen of a surface PO43− (Figure 9c). A more stable binding of water to HAP surface is observed in Figure 9d. The Ow forms an intermolecular bond with Ca2+ while the Hw forms a hydrogen bond with the oxygen of PO43−, resulting in cos θ in the range of −0.5 to 0. Nonetheless, the water dipole in the second layer is less orientated compared with that in the first layer (Figure 8e), which is similar to the observation for the (100) face. The formed ionic and hydrogen bonds are responsible for the strong binding of water in the first hydration layer to the (001) surface, as well as for the suppressed ability for water diffusion. This result further helps to explain the negative solid−water interfacial energy and the observed local hydrophobicity of the (001) surface (Figure 7). Thus, the binding conformation of water to HAP surfaces, as revealed by the analysis of interfacial water order parameter, provides additional insight into the solid−water interfacial energy. H

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Zachara, J. M. Metal Oxide Surfaces and Their Interactions with Aqueous Solutions and Microbial Organisms. Chem. Rev. 1998, 99, 77−174. (2) Hazen, R. M.; Sverjensky, D. A. In The Origins of Life; Deamer, D. W., Szostak, J. W., Eds.; Cold Springs Harbor Perspectives in Biology; Cold Spring Harbor Laboratory Press: Cold Spring Harbor, NY, 2010. (3) Landis, W. J. Mineral Characterization in Calcifying Tissues: Atomic, Molecular and Macromolecular Perspectives. Connect. Tissue Res. 1996, 34, 239−246. (4) Heywood, B. R.; Sparks, N. H.; Shellis, R. P.; Weiner, S.; Mann, S. Ultrastructure, Morphology and Crystal Growth of Biogenic and Synthetic Apatites. Connect. Tissue Res. 1990, 25, 103−119. (5) Fratzl, P.; Gupta, H. S.; Paschalis, E. P.; Roschger, P. Structure and Mechanical Quality of the Collagen-Mineral Nano-Composite in Bone. J. Mater. Chem. 2004, 14, 2115−2123. (6) Landis, W. J. The Strength of a Calcified Tissue Depends in Part on the Molecular Structure and Organization of Its Constituent Mineral Crystals in Their Organic Matrix. Bone 1995, 16, 533−544. (7) Lee, D. D.; Glimcher, M. J. Three-Dimensional Spatial Relationship between the Collagen Fibrils and the Inorganic Calcium Phosphate Crystals of Pickerel (Americanus americanus) and Herring (Clupea harengus) Bone. J. Mol. Biol. 1991, 217, 487−501. (8) Traub, W.; Arad, T.; Weiner, S. Three-Dimensional Ordered Distribution of Crystals in Turkey Tendon Collagen Fibers. Proc. Natl. Acad. Sci. U. S. A. 1989, 86, 9822−9826. (9) Weiner, S.; Traub, W. Bone Structure: From Angstroms to Microns. FASEB J. 1992, 6, 879−885. (10) George, A.; Veis, A. Phosphorylated Proteins and Control over Apatite Nucleation, Crystal Growth, and Inhibition. Chem. Rev. 2008, 108, 4670−4693. (11) Addadi, L.; Weiner, S. Interactions between Acidic Proteins and Crystals: Stereochemical Requirements in Biomineralization. Proc. Natl. Acad. Sci. U. S. A. 1985, 82, 4110−4114. (12) Boskey, A. L. Osteopontin and Related Phosphorylated Sialoproteins: Effects on Mineralization. Ann. N. Y. Acad. Sci. 1995, 760, 249−256. (13) Hu, Y. Y.; Rawal, A.; Schmidt-Rohr, K. Strongly Bound Citrate Stabilizes the Apatite Nanocrystals in Bone. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 22425−22429. (14) Silverman, L. D.; Saadia, M.; Ishal, J. S.; Tishbi, N.; Leiderman, E.; Kuyunov, I.; Recca, B.; Reitblat, C.; Viswanathan, R. Hydroxyapatite Growth Inhibition by Osteopontin Hexapeptide Sequences. Langmuir 2010, 26, 9899−9904. (15) Wu, Y.-J.; Tsai, T. W. T.; Chan, J. C. C. Asymmetric Crystal Morphology of Apatite Induced by the Chirality of Dicarboxylate Additives. Cryst. Growth Des. 2012, 12, 547−549. (16) Wu, Y.-J.; Tseng, Y.-H.; Chan, J. C. C. Morphology Control of Fluorapatite Crystallites by Citrate Ions. Cryst. Growth Des. 2010, 10, 4240−4242. (17) Xie, B.; Nancollas, G. H. How to Control the Size and Morphology of Apatite Nanocrystals in Bone. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 22369−22370. (18) 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. (19) Friddle, R. W.; Battle, K.; Trubetskoy, V.; Tao, J.; Salter, E. A.; Moradian-Oldak, J.; De Yoreo, J. J.; Wierzbicki, A. Single-Molecule Determination of the Face-Specific Adsorption of Amelogenin’s CTerminus on Hydroxyapatite. Angew. Chem., Int. Ed. Engl. 2011, 50, 7541−7545. (20) 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. (21) Corno, M.; Busco, C.; Civalleri, B.; Ugliengo, P. Periodic Ab Initio Study of Structural and Vibrational Features of Hexagonal Hydroxyapatite Ca10(PO4)6(OH)2. Phys. Chem. Chem. Phys. 2006, 8, 2464−2472.

In the present work, the detailed examination of interfacial energies, water structure and dynamics has revealed significantly distinct surface specificity, thus providing a plausible explanation for how biomolecules selectively bind to HAP nanocrystal surfaces and control the morphology.44 Our previous calculations showed that adsorption of both phosphoserine and glutamate to the (100) face is more favorable than to the (001) face.18 According to the planar water density distribution determined in the present study, the entire (100) face is hydrophilic and water can diffuse relatively freely across the surface. When the charged side chain of a biomolecule approaches such a surface, the free energy barrier is lower for the interfacial water to diffuse away. This enables the side chain to penetrate the interfacial layer and directly interact with the surface ions. In contrast water molecules have limited ability to diffuse on the (001) face, thus raising the energy barrier of desolvating surface and preventing direct contact of the surface with a charged side chain. Combined with our previous study,18 the present results provide a thermodynamic basis for the morphology control of inorganic crystals and the adsorption of charged biomolecules to HAP surfaces.

5. CONCLUSION Using well-benchmarked MD simulation, we have shown that the relative interfacial energies of the (100), (110), (001), and (004) faces correspond to their relative surface areas. Detailed analyses of water structure and diffusion at the faces provide atomic-level mechanisms for the behavior of the adsorbed water. The contribution of interfacial water to the binding affinities of charged small biomolecules at HAP faces can also be rationalized based on the results of the present work. Thus, interfacial energies and water behavior play a critical role in the growth of both inorganic and biogenic HAP.

■

ASSOCIATED CONTENT

S Supporting Information *

The interfacial water structures and small molecule adsorption free energies on HAP surfaces obtained by using current HAP− water MD force fields are benchmarked against DFT calculations and experimental measurements from previous literature. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Telephone: 330-972 5795. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Prof. Qiang Cui, University of WisconsinMadison, and Prof. Mesfin Tsige and Ms. Ziqiu Wang, University of Akron, for helpful discussions. Funding was provided by startup funds from University of Akron and NSF DMR Biomat Grant 0906817 to N.S. Computational resources were provided by the Ohio Supercomputer Center (OSC) Grant PBS0286.

■

REFERENCES

(1) Brown, G. E.; Henrich, V. E.; Casey, W. H.; Clark, D. L.; Eggleston, C.; Felmy, A.; Goodman, D. W.; Grätzel, M.; Maciel, G.; McCarthy, M. I.; Nealson, K. H.; Sverjensky, D. A.; Toney, M. F.; I

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

(22) Alves Cardoso, D.; Jansen, J. A.; Leeuwenburgh, S. C. Synthesis and Application of Nanostructured Calcium Phosphate Ceramics for Bone Regeneration. J. Biomed. Mater. Res., Part B 2012, 100, 2316− 2326. (23) Choi, A. H.; Ben-Nissan, B.; Matinlinna, J. P.; Conway, R. C. Current Perspectives: Calcium Phosphate Nanocoatings and Nanocomposite Coatings in Dentistry. J. Dent. Res. 2013, 92, 853−859. (24) Rele, S.; Song, Y.; Apkarian, R. P.; Qu, Z.; Conticello, V. P.; Chaikof, E. L. D-Periodic Collagen-Mimetic Microfibers. J. Am. Chem. Soc. 2007, 129, 14780−14787. (25) Hartgerink, J. D.; Beniash, E.; Stupp, S. I. Self-Assembly and Mineralization of Peptide-Amphiphile Nanofibers. Science 2001, 294, 1684−1688. (26) Flade, K.; Lau, C.; Mertig, M.; Pompe, W. OsteocalcinControlled Dissolution−Reprecipitation of Calcium Phosphate under Biomimetic Conditions. Chem. Mater. 2001, 13, 3596−3602. (27) Tang, R.; Darragh, M.; Orme, C. A.; Guan, X.; Hoyer, J. R.; Nancollas, G. H. Control of Biomineralization Dynamics by Interfacial Energies. Angew. Chem., Int. Ed. Engl. 2005, 44, 3698−3702. (28) Johnsson, M. S. A.; Nancollas, G. H. The Role of Brushite and Octacalcium Phosphate in Apatite Formation. Crit. Rev. Oral Biol. Med. 1992, 3, 61−82. (29) Mahamid, J.; Sharir, A.; Addadi, L.; Weiner, S. Amorphous Calcium Phosphate Is a Major Component of the Forming Fin Bones of Zebrafish: Indications for an Amorphous Precursor Phase. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 12748−12753. (30) Jiang, W.; Chu, X.; Wang, B.; Pan, H.; Xu, X.; Tang, R. Biomimetically Triggered Inorganic Crystal Transformation by Biomolecules: A New Understanding of Biomineralization. J. Phys. Chem. B 2009, 113, 10838−10844. (31) Liu, Y.; Wu, W.; Sethuraman, G.; Nancollas, G. H. Intergrowth of Calcium Phosphates: An Interfacial Energy Approach. J. Cryst. Growth 1997, 174, 386−392. (32) Park, C.; Fenter, P.; Zhang, Z.; Cheng, L.; Sturchio, N. C. Structure of the Fluorapatite (100)-Water Interface by HighResolution X-Ray Reflectivity. Am. Mineral. 2004, 89, 1647−1654. (33) Chiatti, F.; Corno, M.; Ugliengo, P. Stability of the Dipolar (001) Surface of Hydroxyapatite. J. Phys. Chem. C 2012, 116, 6108− 6114. (34) Zhu, W.; Wu, P. Surface Energetics of Hydroxyapatite: A DFT Study. Chem. Phys. Lett. 2004, 396, 38−42. (35) Astala, R.; Stott, M. J. First-Principles Study of Hydroxyapatite Surfaces and Water Adsorption. Phys. Rev. B 2008, 78, 075427. (36) Filgueiras, M. R. T.; Mkhonto, D.; de Leeuw, N. H. Computer Simulations of the Adsorption of Citric Acid at Hydroxyapatite Surfaces. J. Cryst. Growth 2006, 294, 60−68. (37) Zahn, D.; Hochrein, O. Computational Study of Interfaces between Hydroxyapatite and Water. Phys. Chem. Chem. Phys. 2003, 5, 4004−4007. (38) Yang, Y.; Mkhonto, D.; Cui, Q.; Sahai, N. Theoretical Study of Bone Sialoprotein in Bone Biomineralization. Cells Tissues Organs 2011, 194, 182−187. (39) O’Young, J.; Liao, Y.; Xiao, Y.; Jalkanen, J.; Lajoie, G.; Karttunen, M.; Goldberg, H. A.; Hunter, G. K. Matrix Gla Protein Inhibits Ectopic Calcification by a Direct Interaction with Hydroxyapatite Crystals. J. Am. Chem. Soc. 2011, 133, 18406−18412. (40) 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. (41) Lee, W. T.; Dove, M. T.; Salje, E. K. H. Surface Relaxations in Hydroxyapatite. J. Phys.: Condens. Matter 2000, 12, 9829. (42) Mkhonto, D.; de Leeuw, N. H. A Computer Modelling Study of the Effect of Water on the Surface Structure and Morphology of Fluorapatite: Introducing a Ca10(PO4)6F2 Potential Model. J. Mater. Chem. 2002, 12, 2633−2642. (43) Pan, H.; Tao, J.; Xu, X.; Tang, R. Adsorption Processes of Gly and Glu Amino Acids on Hydroxyapatite Surfaces at the Atomic Level. Langmuir 2007, 23, 8972−8981.

(44) Dorvee, J. R.; Veis, A. Water in the Formation of Biogenic Minerals: Peeling Away the Hydration Layers. J. Struct. Biol. 2013, 183, 278−303. (45) 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. (46) Heinz, H. Computational Screening of Biomolecular Adsorption and Self-Assembly on Nanoscale Surfaces. J. Comput. Chem. 2009, 31, 1564−1568. (47) Heinz, H.; Lin, T. J.; Mishra, R. K.; Emami, F. S. Thermodynamically Consistent Force Fields for the Assembly of Inorganic, Organic, and Biological Nanostructures: The Interface Force Field. Langmuir 2013, 29, 1754−1765. (48) Wilson, R. M.; Elliott, J. C.; Dowker, S. E. P. Rietveld Refinement of the Crystallographic Structure of Human Dental Enamel Apatites. Am. Mineral. 1999, 84, 1406−1414. (49) Elliott, J. C.; Mackie, P. E.; Young, R. A. Monoclinic Hydroxyapatite. Science 1973, 180, 1055−1057. (50) Yang, Y.; Xu, Z.; Cui, Q.; Sahai, N. In Biomineralization Sourcebook: Characterization of Biominerals and Biomimetic Materials; DiMasi, E., Gower, L. B., Eds.; CRC Press: Boca Raton, FL, 2014. (51) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (52) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (53) Nose, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (54) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (55) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (56) Darden, T.; Perera, L.; Li, L.; Pedersen, L. New Tricks for Modelers from the Crystallography Toolkit: The Particle Mesh Ewald Algorithm and Its Use in Nucleic Acid Simulations. Structure 1999, 7, R55−60. (57) Marry, V.; Rotenberg, B.; Turq, P. Structure and Dynamics of Water at a Clay Surface from Molecular Dynamics Simulation. Phys. Chem. Chem. Phys. 2008, 10, 4802−4813. (58) Liu, Y.; Nancollas, G. H. Crystallization and Colloidal Stability of Calcium Phosphate Phases. J. Phys. Chem. B 1997, 101, 3464−3468. (59) Lodziana, Z.; Topsoe, N.-Y.; Norskov, J. K. A Negative Surface Energy for Alumina. Nat. Mater. 2004, 3, 289−293. (60) Schneider, J.; Colombi Ciacchi, L. Specific Material Recognition by Small Peptides Mediated by the Interfacial Solvent Structure. J. Am. Chem. Soc. 2011, 134, 2407−2413.

J

dx.doi.org/10.1021/la503158p | Langmuir XXXX, XXX, XXX−XXX