Calcium Phosphate Prenucleation Complexes in Water by Means of ab Initio Molecular Dynamics Simulations

: The nucleation of calcium phosphates, the main inorganic component of bone and tooth tissues, is thought to proceed by aggregation of prenucleation clusters recently identi ﬁ ed as calcium triphosphate complexes. We have performed ab initio molecular dynamics simulations to elucidate their structures and stabilities in aqueous solution. We ﬁ nd the calcium to be seven-coordinated by two water molecules, two bidentate phosphates, and one monodentate phosphate. Free energy results obtained using umbrella sampling simulations show that the complex with a Ca/P ratio of 1:3 is the most energetically favored and more thermodynamically stable than the free ions.


■ INTRODUCTION
Calcium phosphate (CaP) is the major inorganic component of biological hard tissues, (e.g., bone and teeth), where it is mainly present in the form of hydroxyapatite (nominally, Ca 10 (OH) 2 (PO 4 ) 6 ). However, CaP is also responsible for pathological crystallization, leading to several common diseases including atherosclerosis, dental caries, osteoporosis, and kidney stones. 1 Hydroxyapatite 2 originates from an amorphous calcium phosphate (ACP) precursor in a process of dissolution and reprecipitation. 3 In the 1970s, Posner and Betts proposed that structural units of the form Ca 9 (PO 4 ) 6 would aggregate randomly, with the intercluster space filled with water, into larger spherical particles of ACP. 4 More recently, a clustergrowth model for the formation of ACP was proposed, 5 and aggregating clusters were found to be present in body fluids even before nucleation. 6 XANES and XRD experiments of the early stages of CaP crystallization have suggested that an idealized cluster with the formula Ca 9 (PO 4 ) 6 (H 2 O) 30  , where η 1 and η 2 stand for monodentate and bidentate binding of the phosphate groups). 8 Habraken et al. 9 showed that the prenucleation complexes (PNCs), aggregating in solution to form polymeric structures, have the formula [Ca(η 2 −HPO 4 ) 3 ] 4− . Interestingly, a Posner's cluster can be seen as two deprotonated PNCs in which all negative charges are compensated by complexing calcium ions. To explain the formation of ACP at the supersaturations used in their experiment, the PNCs were proposed to have an excess free energy over that of free ions, which dramatically reduces the nucleation barrier. 10 Despite the CaP PNCs representing a key factor for understanding ACP formation, because of the difficulties in accessing experimental data at such small scale, there is no general consensus about their structures. Here, we therefore present a theoretical investigation based on ab initio molecular dynamics (MD) simulations, which suggests that the resulting PNC geometry differs from the existing models in the literature. In addition, we provide an estimate for the free energy cost required to decompose a PNC into isolated free ions, from which we deduce that calcium triphosphate PNCs are more stable entities than other Ca-to-P ratios.

■ THEORETICAL METHODS
The simulations have been carried out using the QUICKSTEP program of CP2K, 11−16 which solves the Kohn−Sham equations of density functional theory (DFT) by combining a Gaussian basis set for the wave functions with an auxiliary plane wave basis set for the density. We have used Goedecker−Teter−Hutter (GTH)-type pseudopotentials 17−19 together with the PBE functional. 20 It is known that gradient-based functionals give an overstructured water system with a lower diffusivity compared to those of experiments. 21,22 Although hybrid functionals can improve the structural and dynamic properties of water, 23 their application in this already highly compute-intensive work would make the simulations beyond reasonable high performance computing capabilities. We note that the adoption of higher simulation temperatures or the inclusion of empirical van der Waals corrections to PBE can lead to a softening of the structure and to higher mobilities. 21,24 However, because it is not clear how they would affect the free energy barriers investigated, we have not pursued these routes in this work.
The [Ca(HPO 4 ) 3 ] 4− PNC was optimized in the gas phase before the resulting structure was solvated with 93 water molecules using the Visual Molecular Dynamics (VMD) solvation tool 25 in a cell of ∼15 × 15 × 15 Å 3 . Two calcium ions were added to neutralize the −4 charge of the complex. The system was first equilibrated in an NPT ensemble for 57.5 ps at a constant temperature of 300 K and 1 bar of pressure to regulate the density of water and then for a further 60 ps in an NVT ensemble (T = 300 K). The time step was set to 0.5 fs.
The stability of the calcium triphosphate complex in water has been evaluated through the analysis of free energy profiles obtained with the umbrella sampling (US) technique. 26,27 In this method, a bias potential is introduced to ensure efficient sampling along reaction coordinate ξ. A series of MD simulations is performed. In each, the system is restrained at a specific value of the reaction coordinate ("target") with the help of the bias potential. The value of the reaction coordinate oscillates around the target, producing a Gaussian-like distribution called the US window. The US windows, centered at increasing values of the target, are then combined with the weighted histogram analysis method (WHAM). 28−31 The potential used in this work was harmonic with strength k where ξ i ref is the target of each window i. The system was equilibrated at each target value for 7.5 ps, and only the last 5 ps were considered for the WHAM analysis. In some cases, more than 2.5 ps were required for equilibrating the system at the new target. A value of 100 kcal mol −1 Å −2 was chosen for k in most of the US windows. Around the minima, lower values of k were sufficient to keep the reaction coordinate oscillating around the target, whereas in the steepest regions of free energy, k was increased up to 150 kcal mol −1 Å −2 .
■ RESULTS AND DISCUSSION Prenucleation Complex in Water. Figure 1(I) shows the structure obtained after the gas phase optimization. Three bidentate phosphate groups are present, giving a calcium coordination number of six. A similar structure with three bidentate phosphate groups, and ligating one water molecule, has been found by both ab initio calculations with a number of explicit water molecules 9 and classical MD simulations in water. 32 However, during our NVT dynamics starting from the PNC depicted in Figure 1(I), two water molecules entered the first coordination sphere and during the simulation never exchanged with other water molecules. For almost 90% of the simulation time, the PNC remained in a Ca (η 2 − HPO 4 2− ) 2 (η 1 −HPO 4 2− )(H 2 O) 2 distorted pentagonal bipyramidal geometry (Figure 1(II)). A similar arrangement of the calcium ligands has been found during classical MD for several types of calcium−phosphate ion pairs. 32 Here, two η 2 − HPO 4 2− groups and one water molecule occupy the equatorial positions, whereas the η 1 −HPO 4 2− oxygen and a second water molecule occupy the axial sites. As the two equatorial groups could undergo a (η 2 −HPO 4 respectively. This result is in line with the findings of Almora-Barrios and de Leeuw for the nucleation of hydroxyapatite at a collagen template by means of classical MD. 32 33 which could be attributed to the lower coordination of calcium (six) or to the less pronounced steric effect of the (bi)carbonate ligand. The Ca−P g(r) shows two maxima signaling bidentate and monodentate phosphate ligands, whose Ca−P distances are 3. 15

Crystal Growth & Design
Article are approximately 0.1 Å larger than their upper limits. 9 In general, the computational data tend to overestimate the distances for Ca−O, P−O, and P−P with respect to the experiments on ACP, 4 but the latter were obtained in the solid phase without solvent effects. Our results are summarized in Table S1 of the Supporting Information, where the coordination numbers, obtained by integrating g(r) between the relevant minima, are consistent with the distorted pentagonal bipyramidal geometry described above.
Stability of Secondary Configurations. To better evaluate the stability of the six-coordinated PNC ( Figure  1(III)) present during approximately 10% of the MD trajectory, we have calculated the free energy profile for the Ca(η 2 − HPO 4 2− ) 2 (η 1 −HPO 4 2− ) to Ca(η 2 −HPO 4 2− )(η 1 −HPO 4 2− ) 2 transition. We show it in the upper panel of Figure 4, where the reaction coordinate is given by the distance between the calcium and a phosphorus from an equatorial group. The minimum of the six-coordinated structure (Figure 1(III)) is only 1.9 kcal/mol higher in energy than that of the sevencoordinated structure (Figure 1(II)). An additional 5 ps simulation with the reaction coordinate restrained at the target, corresponding to the six-coordinated minimum, revealed no significant coordination changes, confirming that a Ca(η 2 − HPO During the entire 60 ps NVT run, two water molecules were tightly ligated to the PNC, whereas we observed a number of attempts by the three phosphate groups to bind simultaneously to calcium in a bidentate mode (Figure 2). To shed light on the possibility of formation of a Ca(η 2 −HPO 4 2− ) 3 PNC, we have studied the free energy profile corresponding to the monodentate-bidentate transition of the central phosphate group. We present it in the lower panel of Figure 4, where the reaction coordinate is the distance between the central calcium and the phosphorus of the monodentate phosphate. We observe that the structure found after equilibration ( Figure  1(II)) is slightly more stable than the structure at shorter Ca−P distances (Figure 1(IV)) in which the central phosphate group is bidentate. However, during the US, we observed the η 2 − cluster. This behavior can be explained by considering the bridging between the phosphate groups: in both minima, the central phosphate is linked to an equatorial phosphate via a direct hydrogen bond and a bridging water and to the other equatorial group via a second bridging water (types (g) and (h), respectively, in Figure S1). For retaining the bridging structures, the lateral bidentate phosphate in the structure in Figure 1(II) had to become monodentate in the structure in Figure 1(IV).
Other than water bridging between phosphate groups, the constant presence of one monodentate ligand can be explained in terms of charge transfer: when the 2+ cation is surrounded by three 2− anions, its partial positive charge is decreased, and therefore, there is a loss of charge−charge interaction with the ligands. 34 The weakening of metal−phosphate interactions increases the importance of other factors, such as the effect of hydrogen bonds with water, which favors the monodentate binding mode, as it maximizes the interaction with the first solvation sphere.
The Ca(η 2 −HPO 4 2− ) 2 (η 1 −HPO 4 2− )(H 2 O) 2 geometry emerging from our MD simulations differs from that suggested by Habraken et al. 9 in the substitution of a (η 2 −HPO 4  4 2− or a water molecule with calcium in a coordination of six. However, we note that a seven-coordinated calcium remains consistent  We set the reaction coordinate to be the distance between the central calcium ion and the phosphorus atom of the η 1 − HPO 4 2− detaching group. Starting from the equilibrated structure (Figure 1(II)), a series of US simulations were performed with increasing target distances up to the complete removal of the phosphate from the Ca 2+ first coordination sphere. The free energy profile depicted in the upper panel of Figure 5 confirms that the PNC is very stable in water, as the energetic cost for moving the η 1 −HPO 4 2− to the second coordination sphere of calcium is approximately 22 kcal/mol. The calcium biphosphate complex in Figure 6(V) appears to be in a shallow minimum located at a distance of 6.4 Å. Because of a third water molecule replacing the η 1 −HPO 4 2− , the calcium coordination is still seven, and the complex is in a Ca(η 2 − HPO 4 2− ) 2 (H 2 O) 3 arrangement. The detached HPO 4 2− group interacts with the calcium biphosphate via two bridging water molecules (type (g) in Figure S1) and a direct hydrogen-bond, as shown in Figure 6 (V). The significant increase in free energy for the dissociation process described above suggests that the PNC is stable with respect to the isolated free ions. To verify this hypothesis, we continued the fragmentation process, removing a second phosphate group from the calcium according to the reaction ) 2 is 2.6 kcal/mol, a value slightly higher than that of the transition in the central phosphate (see Figure 4, lower panel). Furthermore, the position of the minimum for the dissociated structure (5.1 Å) is lower than that in the upper panel of Figure 5 (6.4 Å), where two bulky ligands form hydrogen bonds and water bridges with the leaving group. These interactions oppose the separation of the leaving ion and are reflected in the shape of the free energy curve. Only after breaking the bridges with one of the equatorial phosphates does the curve smooth down to a shallow minimum. The structure with the calcium ion coordinating one bidentate phosphate and five water molecules ( Figure 6(VI)) has already been found during classical molecular dynamics simulations of calcium phosphate complexes. 32 The authors reported a similar arrangement of the   Finally, starting from the minimum in Figure 6(VI), we have detached the last phosphate group, obtaining the isolated ions in solution shown in Figure 6 The corresponding free energy profile is illustrated in the lower panel of Figure 5. Similar to the profile in the middle panel of Figure 5, its main features are a minimum at ∼3.2 Å (bidentate binding mode), a flattening at ∼3.6 Å (monodentate binding mode), and a third depression corresponding to the egress of the phosphate group from the first coordination sphere of the calcium. The free energy barrier is ∼10 kcal/mol, whereas the free energy difference between the monodentate and bidentate binding modes is 2.7 kcal/mol. We note that all the positions and energy differences of the minima involved in reactions 3 and 4 are very similar. The isolated calcium ion shown in Figure  6(VIIa) is six-coordinated, and the structure is a distorted octahedron. Di Tommaso and de Leeuw reported an equivalent structure during Car−Parrinello molecular dynamics simulations of calcium carbonates in water. 33 Despite the similarity, we found a Ca−O wat average distance of 2.45 Å, 0.1 Å longer than theirs, which presumably is due to the different choice of pseudopotentials.
Taking into account the three free energy profiles in Figure 5, we can deduce that each process of a phosphate leaving the calcium coordination shell has an important energetic cost (≥9 kcal/mol) and that the complex is more stable than the free ions. It is worth noting that the energetic expense of the detachment of the first phosphate group was approximately two times more costly than the subsequent detachments, suggesting that calcium triphosphate complexes are more stable entities than other Ca/P ratios, in agreement with experimental data. 9

■ CONCLUSIONS
We have studied the structures and stabilities of [Ca-(HPO 4 ) 3 ] 4− prenucleation complexes in water by means of ab initio molecular dynamics simulations and umbrella sampling techniques. We have found that, in the most stable configuration, one phosphate group ligated to the calcium is monodentate, whereas two are bidentate and placed at opposite sides of the calcium. In addition, two water molecules also bind the calcium ion, which prefers to stay in a (η 2 −HPO 4 2− ) 2 (η 1 − HPO 4 2− )(H 2 O) 2 arrangement. The calcium triphosphate PNC is more stable than the isolated ions, as evidenced by the free energy profiles simulating the dissociation process. Moreover, our data suggest that a Ca/P ratio of 1:3 is thermodynamically favored over other ratios, supporting the experimental findings in the literature. 9 Because nucleation processes take place on a molecular scale, a future theoretical investigation of the aggregation of the clusters would be useful to further support experimental findings, and the interaction of the PNCs with positive counterions (Na + or Ca 2+ ) may improve our understanding of the aggregation of these negatively charged entities. We intend to examine these issues in future work.