Physical origin underlying the prenucleation-cluster-mediated nonclassical nucleation pathways for calcium phosphate

Xiao Yang a, Mingzhu Wang a, Yang Yang b, Beiliang Cui c, Zhijun Xu *a and Xiaoning Yang a
aCollege of Chemical Engineering, State Key Laboratory of Materials-Oriented Chemical Engineering, Nanjing Tech University, Nanjing 210009, China. E-mail:
bDepartment of Chemistry, Lehigh University, 6 East Packer Avenue, Bethlehem, PA 18015, USA
cNetwork Information Center, Nanjing Tech University, Nanjing 210009, China

Received 15th February 2019 , Accepted 2nd April 2019

First published on 3rd April 2019

The involvement of prenucleation clusters (PNCs) in crystallization from a supersaturated solution has been recently admitted within the framework of nonclassical nucleation theory; however, little is known about PNCs, at the quantitative level, for their formation mechanism and stability, the new phase formed by them, as well as their impact on nucleation barriers. Herein, using the sophisticated free energy calculations with a cumulative simulation time of over 5 μs, we identify a thermodynamically favored pathway of the PNC-mediated nucleation for calcium phosphate, starting with the ion pair association in solution. We demonstrate that such an ion association occurs not only between cations and anions, but also for the polyatomic species with charges of the same sign, which, in fact, leads to PNC formation via the consecutive coordination of the phosphate ions to calcium. The free energy decomposition calculations illustrate that the water phase is capable to either hinder or promote ion association for the abovementioned processes, and its specific role is intricately related to the characteristics of the hydration shell around calcium ions. The favorable interactions between the highly charged species play a crucial role in stabilizing the PNC complexes and the aggregates formed by PNCs. Furthermore, our present work reveals that the uptake of an extra calcium ion is the first and mandatory step to trigger PNC aggregation into amorphous calcium phosphate (ACP) by eliminating the related free energy barriers. Our theoretical study successfully provides quantitative explanations to a large set of experimental data in the field, which is currently under intense discussions associated with the nonclassical nucleation mechanism. The combination of computational methods developed in our present work offers a feasible and general solution to quantitatively and systematically study ion associations and crystal nucleation/growth in an aqueous solution at the atomic level, which are normally inaccessible to most of the existing experimental acquisitions.


Crystallization plays a central role in various scientific disciplines such as chemistry, biology, and materials science.1 The formation of a solid phase by crystallization from an aqueous solution always starts with nucleation, arising from the dynamic and stochastic aggregations of ions.2,3 Studies involving the most important biominerals, such as silica,4,5 calcium carbonates,6–11 and calcium phosphates,12–15 along with investigations involving protein crystallization16,17 have revealed the existence of stable solute species, often termed as prenucleation clusters (PNCs) at the initial stage of nucleation.2,11,18 These nanometer-sized clusters can exist in homogeneous solutions and have been observed to serve as building blocks in the formation of a new phase.4–6,12,13,17,19

In particular, calcium phosphate (Ca–P), the major inorganic constituent of bones and teeth with the most thermodynamically stable crystalline structure,20 has attracted tremendous attention in unraveling its mineralization mechanism at the atomic level, with a considerable variety of applications in bone graft substitutes, orthopedics, and dentistry.21–27 Recently, it was demonstrated that the precipitation of Ca–P starts with the aggregation of PNCs, leading to the formation of amorphous calcium phosphate (ACP), which, later, transforms into the crystalline phase.12–15 However, due to the small size of the PNCs, it is extremely challenging to directly visualize them via model characterization techniques. As a result, the structural details of the clusters and their role as building blocks for ACP have been under intense discussion for many years.28,29

In the early 1970s, the “Posner's cluster” with a chemical composition of Ca9(PO4)6 has been postulated as a structural unit of the ACP, as validated by X-ray diffraction studies,30,31 which provided the initial evidence for the cluster-based nucleation of ACP. Later, several experimental studies have confirmed the presence of an isolated unit with an average diameter of 0.7–1.0 nm, which was in good agreement with the size of “Posner's cluster.” This suggested a cluster-growth model involving nanometer-sized PNCs.32–36 In spite of such an agreement, the atomic details of the chemical composition and structures of the proposed Ca–P PNCs have remained elusive.37 The phosphate ions at different protonation states have been discovered in the non-apatitic components of bone and synthetic ACP from several studies.12,13,38–40 It has been further demonstrated that the Ca–P nucleation clusters comprised hydrated calcium and protonated phosphate ions, HPO42−, which played an important role in the formation and structural transitions of ACP, particularly at the very earliest stages of nucleation.39–42 In particular, Habraken et al.13 recently focused on the early-stages apatite nucleation using a range of in situ techniques. They demonstrated that the initial species in the nanometer size is the ion association complex consisting of a single calcium ion and three coordinated phosphate ions (calcium triphosphate), which is proposed to further aggregate into polymeric assemblies via hydrogen bond interactions among multiple clusters. This study represents a significant achievement toward the understanding of the structural details of clusters in the native hydrated state. However, at present, there is still neither any information regarding the thermodynamic stability of such highly charged complexes nor a detailed formation mechanism at the atomic level, such as the pathways leading to ion associations and corresponding kinetic features. Furthermore, the molecular origin underlying the aggregations of these highly charged clusters into subsequent species remains uninvestigated from the perspective of colloidal chemistry and physical chemistry.2

In principle, the ion association in an aqueous solution intricately correlates with the balance between favorable cation–anion interactions and ion–water interactions (e.g., dehydration).42 A quantitative study of the ion association is a nontrivial task, since the association process is speculated to experience a number of metastable states that are strongly related to the subtle rearrangement of the surrounding water molecules.43 For the ion pair association, several states/modes are possible, namely, the double solvent-separated pair (DSSP) (where the anion and cation possess their individual hydration shells), the solvent-shared pair (SSP) (where the two ions share a common hydration shell), and direct contact ion pair (CIP).44 The presence of water molecules can create free energy barriers preventing frequent transitions between these states, and common techniques have difficulties in distinguishing and quantifying these states.45,46 Ions with opposite charges could either associate or repel in an aqueous solution depending on various competing factors, such as Coulomb interaction, dispersion, hydrogen bonding, and the complicated solvation effect.47,48 Currently, ion specificity has posed a challenge on the classical electrolyte theory to explain ion pairing in an aqueous solution.44,49

For the prenucleation cluster studied here, ion association involving multiple-ion interactions presumably experiences a much more intricate rearrangement of the surrounding water molecules; therefore, the complexity increases exponentially from ion pairs to ion clusters and aggregates.43 Even if the advanced experimental techniques have been used to investigate the microscopic structure of the PNCs with remarkable achievements,12,13,32,33,36,50 direct structural information regarding the nature of such species at the angstrom scale is still very scarce. On the other hand, computational methods have served as highly important and supplemental means to theoretically investigate the structure, dynamics, and energetics of condensed-phase systems. For example, molecular dynamics (MD) simulations of a supersaturated solution have successfully demonstrated that the PNCs of calcium carbonate are stable in solution (exist as a strongly hydrated and flexible ionic polymer consisting of chains, branches, and rings).10 Wallace et al.7 further investigated the structure–energetics relationships of the hydrated calcium carbonate clusters and observed a downhill trend in the free energy profile with an increase in cluster size. Very recently, Mancardi et al.51 employed ab initio MD simulations to demonstrate that the calcium triphosphate cluster is thermodynamically stable in their free energy study with umbrella sampling (US) simulations for 5 ps. However, the relatively short simulation time limits the relaxation and rearrangement of water molecules within the ion coordination shell, and the calculated results may significantly depend on the initial cluster configurations, too.45

Herein, we perform comprehensive free energy calculations combined with the available experimental data to investigate, at the atomic level, the formation mechanisms of the ion association complexes of Ca–P and the new phase formed by them in aqueous solutions. The free energy landscapes for the proposed pathways starting from initially free ions in solution to the ACP are carefully computed to identify the most probable PNC-mediated nucleation pathway. The decompositions of the free energies into water-induced contributions and ion–ion interactions allow us to precisely assess the roles of water/ion in Ca–P nucleation. The event triggering the aggregation of highly charged species has been successfully identified by combining large-scale MD simulations at different ion concentrations with specifically designed strategies of free energy calculations. We aim to provide an in-depth understanding of the PNC-mediated nonclassical nucleation pathways for the formation and aggregation of these nanometer-sized species.


General simulation details

The aim of this study is not only to quantitatively characterize the pathway of PNC formation, but also to reveal the molecular mechanism of how these charged complexes aggregate to form the new phase. Therefore, we employ a set of sophisticated free energy calculations and large-scale MD simulations for various conditions. GROMACS 4.5.5 package52 is used to perform all the MD simulations with the periodic boundary condition applied in three directions. The TIP3P model,53 and the classical CHARMM force field,54 are used to describe the water molecules and ions (calcium and phosphate; see below for more details), respectively. These force fields have been successfully used to study the association of calcium and phosphate in the presence of proteins in the condensed phase.55–59 During the preparation of our manuscript, a new force field for Ca–P was successfully developed with universal applications for the solid crystal phase and free ions in the solution.45 When compared to the results reported therein, the force fields used in the present work are capable to effectively reproduce the thermodynamics of ion pairing (see the discussions below).45 It has been demonstrated that the CHARMM force field provides a more realistic description of the mineral ions in an aqueous solution in comparison to that with DS95, a typical force field used to model salt solutions.60 We anticipate, therefore, that the general trends obtained in this study will hold and are independent of force field selection. The van der Waals interactions between the nonbonded atoms are treated by the Lorentz–Berthelot mixing rule. The long-range electrostatic interactions are calculated using the particle mesh Ewald (PME) method.61 A cutoff radius of 12 Å is adopted for the real-space electrostatic interaction and dispersion interaction. The LINCS algorithm is applied to bonds involving hydrogen atoms in all the simulations. The Nosé–Hoover thermostat62 and Parrinello–Rahman method63 are used to control the system temperature and pressure in the NPT ensemble, respectively. An integration time step of 1 fs is used and the trajectories are saved every 1 ps during the production run. The total cumulative simulation time is about 5 μs in the present study (ESI, Table S1). The details regarding each of the simulation systems including the number of water molecules, box size, and specific simulation time are given below.

Binding free energy calculations for PNC formation

Our earlier study has systemically evaluated several well-established computational algorithms used to compute surface–peptide binding free energy.64 The US65 technique, one of the most rigorous free energy calculation strategies, exhibits the best performance for charged interface systems. On the other hand, metadynamics (MTD)—one of the most popular enhanced sampling methods—has attracted considerable attention due to its efficiency and ease of use, particularly for calculating multidimensional free energy. Therefore, these two computational algorithms are employed in our current study to quantify Ca–P binding/association.

Although experimental studies elucidated that the chemical composition of PNC for Ca–P could be calcium triphosphate, the pathway along which the highly charged PNC forms has never been identified. In principle, there are two possible pathways to form an ion association complex: (1) calcium ions could sequentially recruit phosphate ions in three consecutive steps, as shown in Fig. S1a, ESI; (2) two or three phosphate ions simultaneously bind to one calcium ion (Fig. S1b, ESI). Two different strategies, in line with these distinct pathways, for the binding free energy calculations are adopted to characterize the thermodynamics of ion association. Although both hydrogen phosphate and dihydrogen phosphate were detected in Ca–P nucleation, hydrogen phosphate extensively exists in several stages with different phases, such as polymeric assemblies, spheres, and ribbons.13 Therefore, we mainly focus on the ion associations with hydrogen phosphates and subsequent phase formation to unravel the PNC-mediated nucleation mechanism. Therefore, as a special note, “phosphate” in our present work in fact refers to “hydrogen phosphate.”

In Fig. S1a (ESI), where three phosphate ions are proposed to consecutively bind with one calcium ion in a multistep fashion, US together with the weighted histogram analysis method (WHAM),66 is used to calculate the binding free energy. A typical simulation system contains 4136 TIP3P water molecules in a cubic box of 5 × 5 × 5 nm3, which is first equilibrated for 6 ns in an NPT ensemble followed by a 2 ns NVT simulation prior to implementing the US simulations.

The free energy profile as a function of a particular reaction coordinate (RC) is defined as the potential of mean force (PMF). The RC defined here is the separation distance between the calcium ion and center of mass of the incoming phosphate ion for all the three steps. A total number of around 34 window simulations for each binding free energy profile (Table S1, ESI) are performed along the RC with a 0.01 or 0.05 nm interval between the adjacent windows with respective bias potential65 constants of 30[thin space (1/6-em)]000 and 4000 kJ mol−1 nm−2, respectively. The smaller separation distance with a larger umbrella potential applied in the simulations ensures to sufficiently sample the configuration space, including the rearrangement of water molecules around ions. In the present study, the number of water molecules in the solvation shell of ions is not adopted as the second RC in the free energy calculations.45 Previous studies on surface–peptide binding interactions have revealed that the US simulation with a large number of highly overlapping windows is particularly advantageous for configuration sampling in the region close to the surface, where the adsorbates penetrate into the surface-bound water layers with a complicated rearrangement/replacement of water molecules.64,67 Each window simulation is performed for 30 ns in the NVT ensemble, and the error analysis confirms that the obtained results effectively converge, which are suitable for the subsequent free energy decomposition calculations (more details are provided below).

An atomic configuration of the Ca–P complex from an ensemble corresponding to the global free energy minimum is used as the initial structure for subsequent US simulations to further evaluate the free energy landscape for the same Ca2+ to recruit one more HPO42− (Fig. S1, ESI). Additional potential is introduced to restrain the distance between the calcium ion and the already-bound phosphate in the subsequent steps. Such a treatment would minimize the effect due to the structural variations of the existing Ca–P complex; therefore, we could focus on the free energy evaluation for the coordination of the incoming phosphate to calcium.

For the pathway shown in Fig. S1b (ESI), where two or three phosphate ions are simultaneously sequestered to one calcium, we perform well-tempered MTD simulations68 in the NVT ensembles to obtain the multidimensional free energy surfaces (FES) using the PLUMED plugin.69 The setup of the simulation system is similar to that in the US simulations. The RCs are defined as the separation distances between the calcium ion and each of the phosphate ions. Gaussians are deposited every 0.5 ps with a bias factor of 10. The width and initial height of the Gaussian are 0.01 nm and 1 kJ mol−1, respectively, based on our previous experiences and trial tests.64 The MTD simulations are performed for 300 ns, which is sufficiently long to achieve the convergence by monitoring the amplitude of the global free energy minimum. The 1-D free energy profile for the formation of CaHPO4 is also calculated by employing MTD simulation with only one RC to validate its accuracy under the above setup. The obtained result effectively matches that obtained from the US calculations, as shown in Fig. S2 (ESI).

PNC stability

It is particularly worthwhile to investigate whether the ion association complexes (calcium triphosphate) are thermodynamically stable, and therefore, additional MD simulations are further carried out to investigate the possible molecular rearrangements of the PNC along time. Six PNC configurations corresponding to the free energy minimum are extracted from the US trajectories to prepare the initial simulation systems. We remove the distance restraints between the phosphates and calcium within the cluster, which are originally applied in the free energy calculation as described above. All the six simulations are performed for 10 ns with the PNC immersed in the water box (5 × 5 × 5 nm3).

Decomposition of binding free energy

The ion association is dictated by a delicate balance between the direct cation–anion interaction and ion–water interaction. The overall binding free energy profile integrates the thermodynamic contributions from all the governing factors, and therefore, fails to distinguish between the individual energetic contributions toward the ion association. To resolve this problem, we perform the free energy decomposition calculations by combining the US simulations with the force integration,70 which allows us to precisely assess the roles of the water/ion during the association process. This method has been rigorously validated with a successful application in our recent work on a peptide–surface binding process.70 The contribution from the ion–ion interactions (ΔAions) and water-induced contribution (ΔAsol) toward the overall binding free energy could be assessed by integrating the mean forces from the direct cation–anion interaction and water–ion interaction, respectively. The forces exerted on the ions are easily collected as an average over a series of configurations extracted from the US simulations and the mean force effectively converges after 20 ns (Fig. S3, ESI). The convergence analysis within the framework of US could be implemented to monitor whether sufficient sampling is attained to collect the mean forces along the RC. Details about the current free energy decomposition method can be obtained from our earlier study.70

PNC aggregations

In order to address the question of how these highly charged PNCs further aggregate, we perform large-scale simulations designed at various experimental conditions, complemented with additional free energy calculations. Despite the promising development in supercomputers and parallel computing algorithms, nowadays, it is still extremely challenging to directly model the aggregation process of electrolyte aqueous solutions mainly due to the limitations of the applicable simulation timescale. Alternatively, a concentrated Ca–P solution is normally used to accelerate the aggregation process.58,71,72 The experimental study demonstrated that the ACP arising from the aggregation of calcium triphosphate is a fractal of Ca2(HPO4)32−, indicating the potential role of Ca2+ during PNC aggregation.13 We fully consider the available experimental data and prepare four simulation systems as described below.

(1) A simulation box of 10 × 10 × 10 nm3 contains more than 30[thin space (1/6-em)]000 water molecules and 60 PNCs (calcium triphosphate) with 120 free Ca2+ ions added to neutralize the system charge, and therefore, the concentration ratio of calcium/phosphate (Ca/P) is 1.0. The calculated concentration of Ca2+ is about 300 mM, similar to that in several previous computational studies with relatively high ion concentrations.58,71,72 Two different initial configurations are formulated with PNCs and free Ca2+ ions randomly scattered in the box. (2) We prepared a simulation system with excess free calcium ions (Ca/P = 2.0) to investigate the effect of Ca2+ concentration on PNC aggregations. The excess free calcium ions existing in the simulation box lead to a system having a positive charge. Chloride ions (Cl) are, therefore, added to maintain the entire system charge neutral. (3) A reference solution with only PNCs (Ca(HPO4)34−) is also investigated (without any counterions) to further reveal the role of calcium ions in PNC aggregation. It is noted that for such non-neutral systems, artifacts may arise due to the required introduction of a uniform background charge distribution for neutralizing the simulation system in the Ewald algorithm.73 With regard to the reference system, the highly charged PNCs are randomly distributed in the box during the simulation process due to strong repulsive forces (detailed discussions are provided below), resulting in a nearly homogeneous distribution of charges in the real simulation system.73 This could help minimize the artifacts for such a non-neutral system. (4) In order to verify the dependence of the observed aggregation behavior on the ion concentration, we further performed simulations at a much lower concentration, where the ion concentrations are an order of magnitude lower than that in the above systems (i.e., 30 mM vs. 300 mM). All the MD simulations described here are conducted in NPT ensembles for at least 40 ns by monitoring whether the size distributions of clusters are stable over time. If the distance between the calcium ion in one PNC and any phosphorus atom belonging to another PNC is shorter than 4.5 Å, then we define that these two PNCs are considered to be in the same cluster, where the shared solvent molecules have been squeezed out.

To further quantitatively investigate the driving force behind the aggregation of these highly charged species, namely, Ca(HPO4)34−, we compute the free energies, as shown in Fig. S4 (ESI), by the US technique combined with WHAM. The interaction between the two clusters in the solution is first investigated, which is compared to that for the aggregation of two clusters bridged by a free calcium ion (Fig. S4, ESI). An extra spring potential is also used to restrain the distance between calcium and each of the phosphate ions to minimize the possible effect of the conformation variation of each individual PNC on the free energy calculations along the cluster aggregation; however, these charged species are highly stable as discussed in the following sections.

Results and discussion

Formation and stability of PNC

In order to provide a comprehensive description of the Ca–P association in aqueous solutions, the PMF profiles calculated in a multistep fashion (Fig. S1a, ESI) are shown in Fig. 1a, from which the thermodynamic features of these ion association species could be quantitatively estimated. The three PMF profiles, respectively, associated with the formations of CaHPO4, Ca(HPO4)22−, and Ca(HPO4)34− uniformly exhibit two apparent free energy minima separated by a barrier at around 4.5 Å. The free energy minimum at about 5.2 Å corresponds to the SSP state,44 where the ions share two appropriately oriented water molecules, each of which forms a salt bridge with Ca2+ and a hydrogen bond with HPO42− (Fig. 1b). As the ions further approach each other, a free energy barrier of about 2.0–3.0 kcal mol−1 has to be overcome to reach the CIP state, along which the two shared water molecules are squeezed out. The PMF profile captures a split minimum (3.2 and 3.7 Å) for the CIP state in all the three consecutive association steps (Fig. 1b). The structure analysis clearly demonstrates that such a split minimum is due to the two possible coordination states of the incoming phosphate to the calcium ion, namely, monodentate and bidentate coordination (Fig. 1c and d),45 similar to that observed for calcium carbonate association.74 The coexistence of monodentate and bidentate coordination of the phosphate ions to calcium was detected by Zhang et al., when they investigated the Ca–P clusters with X-ray adsorption near-edge spectroscopy.33 The Ca–P PMF profile shown in Fig. 1a effectively reproduces the locations of the minima for the monodentate and bidentate coordination states, as well as the free energy barriers, in comparison with the simulations performed by Demichelis et al.45 The free energy for the CIP state from our present work is also in line with their value45 along with the results from Mancardi et al. with ab initio MD simulations in another work.51 These consistent results indicate that the applied CHARMM force field offers a decent parameterization to study the Ca–P association in an aqueous solution with our present setup.
image file: c9cp00919a-f1.tif
Fig. 1 (a) PMF profiles for three phosphate ions sequentially bound to calcium. Ca, CaP1, and CaP1P2 represent Ca2+, CaHPO4, and Ca(HPO4)22−, respectively, where P1, P2, and P3 refer to the phosphate ions coordinated to calcium at steps 1, 2, and 3, respectively; (b–d) representative atomic configurations of the ion association clusters captured from the free energy minima corresponding to the SSP state (b) and the CIP state (including both monodentate (c) and bidentate (d) coordination of HPO42− to calcium) for each PMF profile. Only the water molecules shared by the ions are shown for the purpose of clarity. Color designations for the atoms are as follows: blue for calcium; red for oxygen; gray for phosphorus; white for hydrogen. The oxygen atoms of the phosphate ion that approach calcium are shown in yellow for facilitating comparison; (e) probability distributions of the three distances between calcium and each of the phosphate ions for calcium triphosphate (Ca(HPO4)34−) from the six independent 10 ns simulations; (f–h) PMF profiles shown in (a) are decomposed into the water-induced contributions (ΔASol) and ion–ion interaction (ΔAions) terms via the force integration method. The error bar is estimated by calculating the 95% confidence interval of the mean value based on the divided four blocks of the simulation trajectories (block averaging). The error bars are smaller than the size of the symbols.

The features identified in the PMF profiles basically exhibit no dependence on the number of HPO42− coordinated to Ca2+, except for the free energy minima at the SSP state, which slightly reduces in magnitude (less negative) as Ca2+ sequentially recruits more HPO42− (Fig. 1a). This suggests that the formations of CaHPO4, Ca(HPO4)22−, and Ca(HPO4)34− experience a similar kinetic and thermodynamic process. In general, the calcium recruiting each phosphate leads to a free energy drop of about −6.5 kcal mol−1, and therefore, the stability of these Ca–P intermediate states are predicted to be Ca(HPO4)34− > Ca(HPO4)22− > CaHPO4. Calcium triphosphate (Ca(HPO4)34−) with each HPO42− exhibited bidentate coordination to Ca2+, which is presumably the most stable ion association complex, with an overall binding free energy of about −19.5 kcal mol−1 with respect to the isolated ions in solution. The dynamics and stability of Ca(HPO4)34− are further examined by additional MD simulations with the initial conformations for the CIP state extracted from the US simulations. Fig. 1e shows the probability distributions of the three Ca–P separation distances. The three phosphate ions equally surround the calcium center, as indicated by these effectively overlapped curves. Two peaks emerge at around 3.2 and 3.7 Å, representing the bidentate and monodentate coordination states of HPO42− to Ca2+, respectively, with a probability of about 90% to be in the bidentate coordination. Furthermore, the Ca–P distance fluctuates in the range of 3.0–3.9 Å, suggesting no rupture on the Ca(HPO4)34− structure over the six independent 10 ns simulations. These results confirm that Ca(HPO4)34− with three bidentate ligands could stably exist in the aqueous solution.

At first glance, it is counterintuitive that the formation of Ca(HPO4)34− from the association of two negatively charged units (Ca(HPO4)22− and HPO42−) is thermodynamically favorable with the same magnitude to CaHPO4 formation from the two oppositely charged ions. In principle, the electrostatic repulsive force should prevent ion association between Ca(HPO4)22− and HPO42−. It is well known that the ion association in an aqueous solution is jointly dominated by the ion–ion interaction and ion–water interaction. Free energy decomposition calculations are further carried out to quantify the roles of water phase and Ca–P interactions along the ion association process.

As shown in Fig. 1f for CaHPO4, the ion pair interaction provides a strong favorable contribution toward the overall binding free energy, as expected (electrostatic attractions). However, the solvent induces a monotonically increasing unfavorable free energy contribution as the two ions approach each other, which is in contrast to the case between the two hydrophobic solute molecules in the aqueous solution.75 In a few studies,48,76–79 it has been accepted that ion pairing in an aqueous solution is driven by an entropy increase due to the release of initially immobilized water molecules in the ion hydration shells into the bulk solution. However, our result may suggest that for a highly charged ionic solute (such as the ones in our study), the entropy gain from water release is completely offset by the unfavorable enthalpy change due to the loss of strong ion–dipole interactions between the released water and ionic solute, which leads to an overall unfavorable water-induced contribution to the PMF (Fig. 1f). This observation is in good agreement with that discussed in Yu et al.'s study involving NaCl pairs.48 Finally, huge cancellation is achieved between the favorable ion–ion attraction and unfavorable water-induced contribution, which results in an overall free energy profile with a relatively small magnitude.48 With regard to the formation of Ca(HPO4)22−, the decomposition study of the binding free energy exhibits a similar trend to that for CaHPO4 (Fig. 1g). However, the cation–anion interaction provides a smaller contribution (nearly halved) to the binding free energy, presumably due to the decrease in the electrostatic attraction between HPO42− and CaHPO4 because of the reduced charge. The water-induced contribution to the PMF, although still unfavorable, correspondingly decreases, too, implying increased competition between enthalpy loss and entropy gain upon the release of water molecules from the ion solvation shells.

The formation of Ca(HPO4)34−, however, exhibits drastically distinct features. The ion pair interaction (Ca(HPO4)22− and HPO42− repulsion), as expected, yields an unfavorable contribution, but the water-induced term alternatively presents a favorable contribution, serving as the driving force toward the formation of Ca(HPO4)34−, which is in striking contrast to the observations for the first two cases. It could be rationalized that the coordination of HPO42− ions around Ca2+ (Fig. 1d) significantly disrupts the cation's solvation shells, and therefore, weakens cation–water interactions. As a result, the entropy gain from the release of water molecules in the hydration shell is proposed to be the dominating factor that drives ion paring (Ca(HPO4)22− and HPO42−) association.79 The free energy decompositions clearly reveal that the ion pairing from two polyatomic anions (e.g., (Ca(HPO4)22− and HPO42−)) could be thermodynamically favorable with a cation as a bridge, and the role of the water phase closely correlates with the characteristics of the water hydration shell of ions. Based on our current study and consistent with intuition, the ion pairing is unexpected to occur between two monatomic ions with the same charge, since entropy gain is insufficient to compensate for the unfavorable enthalpy change.

Pathway for PNC formation

The identification of the most thermodynamically feasible pathway for ion association requires accurate evaluations of the free energy profiles for all the possible pathways, particularly emphasizing upon the free energy barriers between different metastable states. To fulfill this goal, we further conduct multidimensional free energy calculations with two or three separation distances between calcium and the three phosphate ions as the RCs for the pathways, as shown in Fig. S1b (ESI).

The 2-D free energy contours constructed from the MTD simulations with two and three RCs exhibit fairly similar results (Fig. 2a and b). The blue region approximately comprises four blocks representing the bidentate- and monodentate-binding modes between HPO42− and Ca2+; the regions with abrupt color changes (e.g., from blue to yellow/green) indicate the locations of free energy barriers along ion association, as shown in Fig. 2a and b. The MTD simulations confirm that the bidentate coordination of the phosphate ions to calcium is the most stable structure with a free energy of about −13.0 kcal mol−1 for Ca(HPO4)22−, as indicated by the dark-blue region (Fig. 2a). The obtained value is in good agreement with the corresponding US result ((−6.5) + (−6.5) = −13.0 kcal mol−1, as shown in Fig. 1a). This further validates that both the calculation methods are capable to characterize the free energy profiles for ion association involving strong electrostatic interactions and complicated water rearrangement as in our case if a delicate setup for the configurational sampling is properly designed.

image file: c9cp00919a-f2.tif
Fig. 2 Free energies of the Ca–P associations, obtained from MTD simulations performed with multiple RCs, as shown in Fig. S1b (ESI). The RCs are defined as the separation distances between calcium and each of the phosphate ions. (a and b) 2-D free energy contour maps from MTD simulations performed with two RCs and three RCs, respectively. The black triangles indicate the locations of the free energy barriers. The grey dashed lines indicate the pathways along which two or three phosphate ions simultaneously bind to one calcium ion (path 2 and path 3); (c) 1-D PMF profiles along the pathways, directed by the grey dashed lines in (a and b). The PMF profile for Ca2+ consecutively recruiting HPO42− (ion pairing of Ca2+ and HPO42− is taken as an example) is also provided for comparison (path 1).

Special attention needs to be paid to the thermodynamic properties for the simultaneous coordination of two or three phosphate ions to one calcium ion along the pathways labeled by the arrows in Fig. 2a and b. When the ions are dragged together, the ion association needs to cross a free energy barrier of about 5 kcal mol−1 and 9.5 kcal mol−1 for the two and three phosphate ions synchronously binding to calcium, respectively. However, for the pathway shown in Fig. 1a, the PNC formation only needs to overcome a series of small free energy barriers of about 2.5 kcal mol−1. These results suggest that the sequential calcium-recruiting phosphate ions (path 1 in Fig. S1a, ESI) is the most favorable pathway among the three currently proposed pathways for the formation of ion association complexes experimentally observed from a thermodynamics standpoint. As a special note, the three proposed pathways (Fig. S1, ESI) only represent extreme situations for the association mechanism; in reality, the actual pathway may be located somewhere in between, depending on the specific solution condition.

Mechanism of PNC aggregations

Next, we investigate the interesting question of how the negatively charged PNCs further aggregate to form hydrated ACP.13 Inspired from the observed favorable interaction between Ca(HPO4)22− and HPO42−, we perform the simulation for a reference system containing only PNCs and water molecules in a cubic box to investigate whether a similar thermodynamic trend holds true for the aggregation of Ca(HPO4)34−. It should be noted that the simulation system is not charge-neutral, and the details about the simulation setup can be found in the Methods section. The Ca/P ratio in the reference system is in complete accordance with that in the experimentally identified polymeric assemblies (0.3–0.4).13,33 As shown in Fig. 3a, no aggregation occurs for these highly charged PNCs, which maintain random distributions in the simulation box after 40 ns. The statistical analysis involving the formed aggregates reveals that all the PNCs are solvated in an isolated manner with zero collision probability during the entire simulation period (Fig. 3b). This could be rationalized by a nearly monotonic increase in the PMF profile (Fig. 3c) as the two PNCs are pushed together. These results indicate that the electrostatic repulsion dictates the overall free energy change as the PNCs approach each other, which is in contrast to that for the favorable pairing interaction between Ca(HPO4)22− and HPO42− discussed above.80 The shoulder (around 4.4 Å) on the PMF suggests a local intermediate state with free energy of about +3.4 kcal mol−1, where several salt bridges (Ca–O) are observed between calcium in one PNC and phosphate ions from the other PNC. However, such an intermediate state is thermodynamically unstable because it exists in the shallow free energy well. Our results suggest that the thermodynamically favorable factors, e.g., the release of water molecules or salt bridge formations, cannot compensate for the significant electrostatic repulsion between the PNCs.
image file: c9cp00919a-f3.tif
Fig. 3 (a) Snapshots taken at 0 and 40 ns, for the reference system containing only PNCs and water molecules. The color designation is the same as that in Fig. 1. Water molecules have been omitted for the purpose of clarity; (b) probabilities to observe PNCs in aggregations (blue line) and free PNCs in an aqueous solution (red line) as a function of time; (c) PMF profile for the PNC–PNC interaction as they approach each other. A local free energy well is marked by a red circle, where a representative structure of the PNC pairing is provided as the graphical inset. To distinguish between the two PNCs, the same color designation is used, except for the incoming Ca(HPO4)34−, where Ca2+ and oxygen atoms in HPO42− are displayed in green and yellow, respectively.

Therefore, we hypothesize that the neutralization of the highly charged species is the first and mandatory step toward PNC aggregation. Under the experimental conditions, the presence of counterions, such as H+, Na+, and Ca2+, presumably could alleviate electrostatic repulsions. Moreover, calcium has been proposed to participate in the formation of ACP as well as via PNC aggregations in solution.13,33

To verify the current proposal, we further perform two additional MD simulations with different Ca/P ratios (Ca/P = 1.0 and 2.0). The extracted snapshots illustrate that PNCs successfully aggregate into two assemblies by using Ca2+ within the system with Ca/P = 1.0 (Fig. 4a). For a calcium-rich system (Ca/P = 2.0), it is observed that although the self-assemblies of PNCs indeed occur, a number of small aggregates (each consisting of only a few PNCs) scatter in the box. This suggests that the excess calcium ions bound to the PNCs make the small aggregates have a positive charge, which prevents further aggregation.

image file: c9cp00919a-f4.tif
Fig. 4 (a) Snapshots of the PNC aggregation (Ca/P = 1.0) taken at 0 and 40 ns; (b) snapshots of the PNC aggregation (Ca/P = 2.0) taken at 0 and 40 ns. The color designation is the same as that in Fig. 3. Chloride as the counterion is shown in dark gray. Water molecules have been omitted for clarity; (c) the number of aggregates and hydrogen bonds formed during the self-assembly process as a function of time for the two investigated systems. Fitting lines are drawn to facilitate visualization. The numbers “1” and “2” refer to the simulation systems with Ca/P = 1.0 and 2.0, respectively; (d and e) probabilities to find free Ca2+ and free PNCs (i.e., not in aggregates) as a function of time.

The dynamic process of PNC aggregation under two different conditions could be elucidated by tracing the amount of aggregates along time. The analysis for the system with Ca/P = 1.0 demonstrates that during the first stage of rearrangement (∼6 ns), the PNCs rapidly cluster into about 12 individual aggregates (Fig. 4c). The fast first stage is followed by a slower process, during which these small aggregates merge to form two large ones within 30 ns (Fig. 4a and c). From the probability plot to observe free Ca2+/PNCs as a function of time (Fig. 4d), we see a faster decay for free Ca2+ in comparison with free PNCs, suggesting that the binding of Ca2+ to the PNCs occurs before PNC aggregation. These results indicate that calcium ions may promote the self-assembly of PNCs. Finally, all the PNCs and calcium ions are involved in the formation of two aggregates (Fig. 4a and d). The simulations performed at an ion concentration ten times lower than the above ones also display the aggregation of PNCs into small clusters (Fig. S5, ESI), confirming that the observed trend is independent of ion concentration.

In the presence of excess Ca2+, the PNCs quickly aggregate to form small clusters via the inclusion of Ca2+ at the first stage (Fig. 4c and e), exhibiting similar self-assembly dynamics as those at low Ca2+ concentrations. However, the descending trend regarding the number of aggregates is not observed after the fast stage, indicating no further merging of these small aggregates into large assemblies during the simulations (Fig. 4b). Our calculations show that 90% of the PNCs and about 60% of the calcium ions aggregate into 14 clusters (Fig. 4b and e), suggesting a much slower self-assembly process within the current simulation period, which is also consistent with the result from the hydrogen bond analysis within the aggregates (Fig. 4c). Further, we observe a much slower self-assembly process for the calcium-insufficient system (Ca/P = 0.5), where sodium ions are added to maintain the entire system charge neutral. This could be rationalized by the presence of insufficient calcium ions to reduce the free energy barrier associated with PNC aggregations. This also indicates that compared to calcium, sodium cannot efficiently initiate PNC aggregation.

The present results highlight the crucial role of Ca2+ in driving the aggregation of highly charged PNCs, which is in good agreement with the experimental data.13 It could be concluded that the aggregation of PNCs is triggered by the presence of free calcium ions in solution, which strongly bind to PNCs, as supported by the binding PMF profile (Fig. 5a). We find an overall downhill free energy change along Ca(HPO4)34− and Ca2+ association. Several shoulders on the free energy profile are observed, which are related to the formation of salt bridges between the oxygen atoms of the phosphate and incoming calcium (Fig. 5c). For example, a sudden free energy drop is observed from 7.0 to 6.7 Å, which is caused by the monodentate coordination between one of the phosphate ions in Ca(HPO4)34− to free calcium. Another significant free energy drop occurs in the range from 5.8 to 3.9 Å, which is mainly due to the formation of a fairly stable complex, where the two calcium ions symmetrically reside at opposite sides of the plane formed by the three phosphate ions (Fig. 5c). Such a binding mode maximizes the number of possible salt bridges (five between the phosphate ion and each of the two calcium ions with two bidentate and one monodentate coordination). The coexistence of two coordination models for the compact Ca–P structure was also experimentally and theoretically suggested for the aggregates formed by PO43− and Ca2+.33,81 In summary, the free energy profile clearly demonstrates that the ion pairing of Ca(HPO4)34− and Ca2+ through a sequence of intricate structural transitions is a thermodynamically spontaneous process with no apparent free energy barriers (Fig. 5a), leading to a rapid disappearance of free Ca2+ in solution (Fig. 4d).

image file: c9cp00919a-f5.tif
Fig. 5 (a) PMF profile for the formation of Ca2(HPO4)32− resulting from the coordination of Ca(HPO4)34− to Ca2+; (b) PMF profile for the formation of Ca3(HPO4)66− from the pair interaction between Ca2(HPO4)32− and Ca(HPO4)34−; (c and d) molecular structures of the complexes corresponding to the states indexed on the PMF curves. The color designation is the same as that in Fig. 1, except for calcium ions being pushed to PNC, where bridged Ca2+ and the one in the incoming Ca(HPO4)34− are displayed in green for the convenience of identification.

As an interesting observation, although the newly formed ion cluster, termed as post-nucleation aggregates13 with a composition of Ca2(HPO4)32−, still possesses a charge of −2; the free energy calculation reveals that it binds very favorably to Ca(HPO4)34− with a free energy drop of about −9.5 kcal mol−1 at 4.2 Å (Fig. 5b). A similar multistep binding process is noticed, which strongly correlates to the formations of different coordination models between the phosphate ions and three calcium ions, including the breaking of the pre-formed monodentate coordination corresponding to the free energy barriers with small magnitudes (Fig. 5b and d). The atomic conformation of the formed complex corresponding to the free energy minimum reveals that the three calcium ions of Ca3(HPO4)66− are distributed as an equilateral triangle topology to form salt bridges with phosphates to the maximum extent (Fig. 5b and d), which certainly stabilizes the aggregated structure by neutralizing the unfavorable electrostatic repulsions between the two PNCs.

Our simulations clearly demonstrate that PNC aggregation is initiated by the formation of Ca2(HPO4)32−via Ca(HPO4)34− and Ca2+ association. With the presence of Ca2+, the free energy of the self-assembly of two PNCs (compared to the respective isolated reference state) dramatically drops from +3.2 kcal mol−1 (Fig. 3c) to −9.5 kcal mol−1 (Fig. 5b) at the same separation distance of 4.2 Å. Under the relevant experimental conditions, more calcium ions presumably could bind with calcium triphosphate prior to aggregation, resulting in a varied free energy landscape as compared to that discussed above. These binding calcium ions help neutralize the negatively charged complexes and are in favor of the subsequent ACP formation. The present work theoretically demonstrates that the low calcium ion concentration (Ca/P = 1.0) is thermodynamically sufficient to trigger the aggregation of highly charged polyatomic complexes. Our calculation results successfully explain, from a free energy perspective, the experimental observation that Ca2(HPO4)32− acts as the nanometer-sized building block of ACP.2,13 In addition, our results also shed light on the formation of branched polymeric assemblies, as observed by Habraken et al., who proposed a reaction-limited aggregation process involving hydrogen bond formation between the two complexes and the release of hydration water.13 However, these two thermodynamically favorable contributions are not suggested to account for the stabilization of polymeric clusters in the current study. Alternatively, other cations, such as H+ and Na+, are expected to behave similarly as Ca2+ to induce PNC aggregation.


In the present study, MD simulations are applied to quantitatively identify PNC-based pathways for early-stage Ca–P nucleation. In particular, based on the available experimental data, the RC-specified strategies for the free energy calculations are employed to characterize the energetic and kinetic features of nonclassical nucleation from the initially free ions in solution to the eventually formed amorphous phase.

Our computational results have been extensively compared to the available experimental observations and they successfully yield quantitative explanations toward a series of dilemmas concerning the dynamics and stability of the so-called ion association clusters and the new phase formed by them, as well as their impact on nucleation barriers within the nonclassical nucleation theory framework. Based on our theoretical calculations, we propose the mechanism underlying the formation of the experimentally observed calcium triphosphate, and we demonstrate that the most favorable pathway is a sequential one, where the calcium ion recruits the phosphate ion one by one during Ca–P nucleation. Our study also presents thermodynamic evidence that calcium triphosphate with HPO42− exhibiting bidentate coordination with Ca2+ could exist in an aqueous solution in a very stable manner with binding free energy of around −19.5 kcal mol−1. In addition, the current study determines low free energy barriers of about 2.5 kcal mol−1 in the PNC-based nucleation mechanism, where the favorable associations of the charged species play a crucial role in PNC formations and their aggregation into the ACP. Moreover, we also demonstrate that the uptake of an extra calcium ion helps trigger the aggregation of highly charged PNCs, which significantly reduces the related free energy barriers from +3.2 to −9.5 kcal mol−1.

The developed techniques in our present work on the basis of the pathway-designed free energy calculations combined with the available experimental data offer a feasible and general solution to quantitatively and systematically study ion associations and crystal nucleation/growth in an aqueous solution at the atomic level, which is normally inaccessible to most of the existing experimental acquisitions. The current study also establishes a worthwhile foundation to further investigate the more complicated ion association clusters involving phosphates with different protonation states, such as the experimentally detected PO43− and H2PO4, which ultimately leads to a comprehensive understanding of the molecular mechanism underlying the nonclassical nucleation process and phase separation.

Conflicts of interest

There are no conflicts to declare.


We gratefully acknowledge the financial support provided by the National Natural Science Foundation of China (Grant 21771106 and 21676136), and the Natural Science Foundation of Jiangsu Province, P. R. China (Grant BK20160995). The computational resources generously provided by the High Performance Computing Center of Nanjing Tech University are greatly appreciated.

Notes and references

  1. J. J. De Yoreo, P. U. Gilbert, N. A. J. M. Sommerdijk, R. L. Penn, S. Whitelam, D. Joester, H. Zhang, J. D. Rimer, A. Navrotsky and J. F. Banfield, Science, 2015, 349, aaa6760 CrossRef PubMed .
  2. D. Gebauer, M. Kellermeier, J. D. Gale, L. Bergström and H. Cölfen, Chem. Soc. Rev., 2014, 43, 2348 RSC .
  3. J. J. De Yoreo, Nat. Mater., 2013, 12, 284 CrossRef CAS PubMed .
  4. M. Kellermeier, E. Melero-García, F. Glaab, R. Klein, M. Drechsler, R. Rachel, J. M. García-Ruiz and W. Kunz, J. Am. Chem. Soc., 2010, 132, 17859 CrossRef CAS PubMed .
  5. A. R. F. Herman Cho, J. P. K. Raluca Craciun and D. D. A. Neil Shah, J. Am. Chem. Soc., 2006, 128, 2324 CrossRef PubMed .
  6. D. Gebauer, A. Volkel and H. Colfen, Science, 2008, 322, 1819 CrossRef CAS PubMed .
  7. A. F. Wallace, L. O. Hedges, A. Fernandez-Martinez, P. Raiteri, J. D. Gale, G. A. Waychunas, S. Whitelam, J. F. Banfield and J. J. De Yoreo, Science, 2013, 341, 885 CrossRef CAS PubMed .
  8. E. M. Pouget, P. H. H. Bomans, J. A. C. M. Goos, P. M. Frederik, G. D. With and N. A. J. M. Sommerdijk, Science, 2009, 323, 1455 CrossRef CAS PubMed .
  9. D. Gebauer, P. N. Gunawidjaja, J. Y. Ko, Z. Bacsik, B. Aziz, L. Liu, Y. Hu, L. Bergström, C. W. Tai and T. K. Sham, Angew. Chem., Int. Ed., 2010, 49, 8889 CrossRef CAS PubMed .
  10. R. Demichelis, P. Raiteri, J. D. Gale, D. Quigley and D. Gebauer, Nat. Commun., 2011, 2, 590 CrossRef PubMed .
  11. D. Gebauer and H. Cölfen, Nano Today, 2011, 6, 564 CrossRef CAS .
  12. A. Dey, P. H. Bomans, F. A. Muller, J. Will, P. M. Frederik, G. D. With and N. A. J. M. Sommerdijk, Nat. Commun., 2010, 9, 1010 CAS .
  13. W. J. E. M. Habraken, J. Tao, L. J. Brylka, H. Friedrich, L. Bertinetti, A. S. Schenk, A. Verch, V. Dmitrovic, P. H. H. Bomans, P. M. Frederik, J. Laven, P. van der Schoot, B. Aichmayer, G. de With, J. J. DeYoreo and N. A. J. M. Sommerdijk, Nat. Commun., 2013, 4, 1507 CrossRef PubMed .
  14. D. Kim, B. Lee, S. Thomopoulos and Y.-S. Jun, Cryst. Growth Des., 2016, 16, 5359 CrossRef CAS PubMed .
  15. D. Kim, B. Lee, S. Thomopoulos and Y.-S. Jun, Nat. Commun., 2018, 9, 962 CrossRef PubMed .
  16. G. Bitan, A. Lomakin and D. B. Teplow, J. Biol. Chem., 2001, 276, 35176 CrossRef CAS PubMed .
  17. S. Fermani, C. Vettraino, I. Bonacini, M. Marcaccio, G. Falini, J. A. Gavira and J. M. G. Ruiz, Cryst. Growth Des., 2013, 13, 3110 CrossRef CAS .
  18. P. Smeets, A. R. Finney, W. J. E. M. Habraken, F. Nudelman, H. Friedrich, J. Laven, J. J. DeYoreo, P. M. Rodger and N. A. J. M. Sommerdijk, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E7882 CrossRef CAS PubMed .
  19. J. Baumgartner, A. Dey, P. H. H. Bomans, C. Le Coadou, P. Fratzl, N. A. J. M. Sommerdijk and D. Faivre, Nat. Mater., 2013, 12, 310 CrossRef CAS PubMed .
  20. W. A. House, Environ. Technol., 1999, 20, 727 CrossRef CAS .
  21. Y. Wang, T. Azaïs, M. Robin, A. Vallée, C. Catania, P. Legriel, G. Pehau-Arnaudet, F. Babonneau, M.-M. Giraud-Guille and N. Nassif, Nat. Mater., 2012, 11, 724 CrossRef CAS PubMed .
  22. S. Elsharkawy, M. Al-Jawad, M. F. Pantano, E. Tejeda-Montes, K. Mehta, H. Jamal, S. Agarwal, K. Shuturminska, A. Rice, N. V. Tarakina, R. M. Wilson, A. J. Bushby, M. Alonso, J. C. Rodriguez-Cabello, E. Barbieri, A. del Río Hernández, M. M. Stevens, N. M. Pugno, P. Anderson and A. Mata, Nat. Commun., 2018, 9, 2145 CrossRef PubMed .
  23. W. Habraken, P. Habibovic, M. Epple and M. Bohner, Mater. Today, 2016, 19, 69 CrossRef CAS .
  24. A. Sato and A. Iwasaki, Biomaterials, 2004, 25, 3569 CrossRef PubMed .
  25. G. Daculsi, J. M. Bouler and R. Z. LeGeros, Int. Rev. Cytol., 1997, 172, 129 CAS .
  26. M. Bohner, S. Tadier, N. V. Garderen, A. D. Gasparo, N. Döbelin and G. Baroud, Biomatter, 2013, 3, e25103 CrossRef PubMed .
  27. L.-N. Niu, S. E. Jee, K. Jiao, L. Tonggu, M. Li, L.-G. Wang, Y.-D. Yang, J.-H. Bian, L. Breschi, S. S. Jang, J.-H. Chen, D. H. Pashley and F. R. Tay, Nat. Mater., 2016, 16, 370 CrossRef PubMed .
  28. L. Wang and G. H. Nancollas, Dalton Trans., 2009, 2665 RSC .
  29. D. Erdemir, A. Y. Lee and A. S. Myerson, ChemInform, 2009, 40, 621 CrossRef .
  30. F. Betts and A. S. Posner, Mater. Res. Bull., 1974, 9, 353 CrossRef CAS .
  31. A. S. Posner and F. Betts, Acc. Chem. Res., 1974, 8, 273 CrossRef .
  32. L.-W. Du, S. Bian, B.-D. Gou, Y. Jiang, J. Huang, Y.-X. Gao, Y.-D. Zhao, W. Wen, T.-L. Zhang and K. Wang, Cryst. Growth Des., 2013, 13, 3103 CrossRef CAS .
  33. Q. Zhang, Y. Jiang, B.-D. Gou, J. Huang, Y.-X. Gao, J.-T. Zhao, L. Zheng, Y.-D. Zhao, T.-L. Zhang and K. Wang, Cryst. Growth Des., 2015, 15, 2204 CrossRef CAS .
  34. L. Wang, S. Li, E. Ruiz-Agudo, C. V. Putnis and A. Putnis, CrystEngComm, 2012, 14, 6252 RSC .
  35. K. Onuma and A. Ito, Chem. Mater., 1998, 10, 3346 CrossRef CAS .
  36. F. Nudelman, K. Pieterse, A. George, P. H. H. Bomans, H. Friedrich, L. J. Brylka, P. A. J. Hilbers, G. de With and N. A. J. M. Sommerdijk, Nat. Mater., 2010, 9, 1004 CrossRef CAS PubMed .
  37. G. Mancardi, C. E. Hernandez Tamargo, D. Di Tommaso and N. H. de Leeuw, J. Mater. Chem. B, 2017, 5, 7274 RSC .
  38. Y. Wu, M. J. Glimcher, C. Rey and J. L. Ackerman, J. Mol. Biol., 1994, 244, 423 CrossRef CAS PubMed .
  39. S. V. Dorozhkin, Acta Biomater., 2010, 6, 4457 CrossRef CAS PubMed .
  40. Y. Wang, S. Von Euw, F. M. Fernandes, S. Cassaignon, M. Selmane, G. Laurent, G. Pehau-Arnaudet, C. Coelho, L. Bonhomme-Coury, M.-M. Giraud-Guille, F. Babonneau, T. Azaïs and N. Nassif, Nat. Mater., 2013, 12, 1144 CrossRef CAS PubMed .
  41. J. D. Termine and E. D. Eanes, Calcif. Tissue Res., 1972, 10, 171 CrossRef CAS PubMed .
  42. J. Vincze, M. Valisko and D. Boda, J. Chem. Phys., 2010, 133, 154507 CrossRef PubMed .
  43. N. F. A. van der Vegt, K. Haldrup, S. Roke, J.-R. Zheng, M. Lund and H. J. Bakker, Chem. Rev., 2016, 116, 7626 CrossRef CAS PubMed .
  44. Y. Marcus and G. Hefter, Chem. Rev., 2006, 106, 4585 CrossRef CAS PubMed .
  45. R. Demichelis, N. A. Garcia, P. Raiteri, M. R. Innocenti, C. L. Freeman, J. H. Harding and J. D. Gale, J. Phys. Chem. B, 2017, 122, 1471 CrossRef PubMed .
  46. J. Hack, D. C. Grills, J. R. Miller and T. Mani, J. Phys. Chem. B, 2016, 120, 1149 CrossRef CAS PubMed .
  47. L. Mikael, V. Robert and J. Pavel, Langmuir, 2008, 24, 3387 CrossRef PubMed .
  48. S. Yu and T. L. Beck, J. Phys. Chem. B, 2017, 121, 2189 CrossRef PubMed .
  49. W. Kunz, Specific ion effects, World Scientific Pub. Co., Singapore, Hackensack, NJ, 2010 Search PubMed .
  50. E. Andreas, N. Andreas, H. Glenn and B. Richard, J. Phys. Chem. B, 2015, 119, 5270 CrossRef PubMed .
  51. G. Mancardi, U. Terranova and N. H. de Leeuw, Cryst. Growth Des., 2016, 16, 3353 CrossRef CAS .
  52. B. Hess, C. Kutzner, D. V. D. Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435 CrossRef CAS PubMed .
  53. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926 CrossRef CAS .
  54. B. R. Brooks, C. L. Brooks, A. D. Mackerell, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J.-Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York and M. Karplus, J. Comput. Chem., 2009, 30, 1545 CrossRef CAS PubMed .
  55. W.-L. Zhao, Z.-Q. Wang, Z.-J. Xu and N. Sahai, Phys. Chem. Chem. Phys., 2018, 20, 13047 RSC .
  56. R. Pokhrel, B. S. Gerstman, J. D. Hutcheson and P. P. Chapagain, J. Phys. Chem. B, 2018, 122, 3782 CrossRef CAS PubMed .
  57. W.-L. Zhao, Z.-J. Xu, Y. Yang and N. Sahai, Langmuir, 2014, 30, 13283 CrossRef CAS PubMed .
  58. J. Ma, J. Mater. Sci., 2014, 49, 3099 CrossRef CAS .
  59. Z.-J. Xu, Y. Yang, W.-L. Zhao, Z.-Q. Wang, W. J. Landis, Q. Cui and N. Sahai, Biomaterials, 2015, 39, 59 CrossRef CAS PubMed .
  60. N. V. Lukasheva and D. A. Tolmachev, Langmuir, 2016, 32, 125 CrossRef CAS PubMed .
  61. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089 CrossRef CAS .
  62. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695 CrossRef .
  63. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182 CrossRef CAS .
  64. Q.-C. Wei, W.-L. Zhao, Y. Yang, B.-L. Cui, Z.-J. Xu and X.-N. Yang, ChemPhysChem, 2018, 19, 668 CrossRef CAS .
  65. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187 CrossRef .
  66. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011 CrossRef CAS .
  67. Z.-J. Xu, Q.-C. Wei, W.-L. Zhao, C. Qiang and N. Sahai, J. Phys. Chem. C, 2018, 122, 4372 CrossRef CAS .
  68. A. Barducci, G. Bussi and M. Parrinello, Phys. Rev. Lett., 2008, 100, 020603 CrossRef PubMed .
  69. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci and R. A. Broglia, Comput. Phys. Commun., 2009, 180, 1961 CrossRef CAS .
  70. Z.-J. Xu, X. Yang, Q.-C. Wei, W.-L. Zhao, B.-L. Cui, X.-N. Yang and N. Sahai, Langmuir, 2018, 34, 7932 CrossRef CAS PubMed .
  71. Y. Yang, Q. Cui and S. Nita, Langmuir, 2010, 26, 9848 CrossRef CAS PubMed .
  72. N. Almora-Barrios and N. H. de Leeuw, Cryst. Growth Des., 2012, 12, 756 CrossRef CAS .
  73. J. S. Hub, B. L. de Groot, H. Grubmüller and G. Groenhof, J. Chem. Theory Comput., 2014, 10, 381 CrossRef CAS PubMed .
  74. K. Henzler, E. O. Fetisov, M. Galib, M. D. Baer, B. A. Legg, C. Borca, J. M. Xto, S. Pin, J. L. Fulton and G. K. Schenter, Sci. Adv., 2018, 4, eaao6283 CrossRef PubMed .
  75. E. Sobolewski, M. Makowski, C. Czaplewski, A. Liwo, S. Ołdziej and H. A. Scheraga, J. Phys. Chem. B, 2007, 111, 10765 CrossRef CAS PubMed .
  76. Y. Marcus, J. Solution Chem., 1987, 16, 735 CrossRef CAS .
  77. J. Iwahara, A. Esadze and L. Zandarashvili, Biomolecules, 2015, 5, 2435 CrossRef CAS PubMed .
  78. K. D. Collins, Methods, 2004, 34, 300 CrossRef CAS PubMed .
  79. Q. Zhang, B.-B. Zhang, L. Jiang and W. Zhuang, Chin. J. Chem. Phys., 2013, 26, 694 CrossRef CAS .
  80. B. Cantaert, Y. Y. Kim, H. Ludwig, F. Nudelman, N. A. J. M. Sommerdijk and F. C. Meldrum, Adv. Funct. Mater., 2012, 22, 907 CrossRef CAS .
  81. T.-J. Lin and C.-C. Chiu, Phys. Chem. Chem. Phys., 2018, 20, 345 RSC .


Electronic supplementary information (ESI) available. See DOI: 10.1039/c9cp00919a

This journal is © the Owner Societies 2019