Simulation of Calcium Phosphate Prenucleation Clusters in Aqueous Solution: Association beyond Ion Pairing

Classical molecular dynamics simulations and free energy methods have been used to obtain a better understanding of the molecular processes occurring prior to the first nucleation event for calcium phosphate biominerals. The association constants for the formation of negatively charged complexes containing calcium and phosphate ions in aqueous solution have been computed, and these results suggest that the previously proposed calcium phosphate building unit, [Ca(HPO4)3]4–, should only be present in small amounts under normal experimental conditions. However, the presence of an activation barrier for the removal of an HPO42– ion from this complex indicates that this species could be kinetically trapped. Aggregation pathways involving CaHPO4, [Ca(HPO4)2]2–, and [Ca(HPO4)3]4– complexes have been explored with the finding that dimerization is favorable up to a Ca/HPO4 ratio of 1:2.


Determination of the lowest free energy structure
In the cases where an ion association involved only one Ca, the lowest free energy structure was determined as simply the lowest free energy as a function of both the distance and Ca coordination number (CN). For associations involving two Ca, the free energy was determined as a function of three collective variables (CVs). Therefore, the free energy was compared as a function of each CV alone or as a combination of each pair of 2 CVs (where the remaining CV(s) are integrated out). For example, in the case of [Ca(HPO4)2] -2 + Ca +2 → [(CaHPO4)2], the three 3-dimensional plots (free energy vs 2 CVs) in Figure S3 corresponding to all combinations of the three CVs yield a consistent picture of the free energy landscape and minimum. The free energy as a function of the ion Ca CN and distance indicates that the minimum is around 3.5 Å with an ion Ca CN of 4. As a function of the cluster Ca CN and distance, it indicates that the minimum is around 3.5 Å with a cluster Ca CN of 4. As a function of the ion Ca CN and cluster Ca CN, it indicates that both CNs should be 4. All of this leads to a structure which should have 4 water oxygens coordinated to each Ca and a Ca-Ca distance of approximately 3.5 Å.
In Figures S1-S20 the free energy profiles are given as a function of 1 and 2 CVs for all of the ion association processes studied. In the case of the 1-D free energy profiles versus distance then a dashed black line is included to indicate the asymptotic long-range limit used to align the free energy and check convergence. Table S1. The association reaction performed and the corresponding standard free energy, time simulated, biased CV(s) and walled CV(s) for each process, where "d" refers to distance, "CN" refers to coordination number and "COM" refers to the center of mass of a cluster. In the "Walled CV" column, all distance walls are upper walls, all coordination number walls are lower walls and all subscripts denote which atom the wall occurs between, corresponding to numerical appearance in the association reaction as labeled. The bias factor used was 5 except in the cases of 8, 11, 12, 17 and 18 where the bias factor was 15.   Table S2. The standard free energies for formation of each ion cluster. Where multiple pathways are available, then the free energy for each association mechanism is given and the uncertainty is reported as the standard error.

Species
Association  Table S3. The free energies for formation of each ion cluster as reported by Yang et al. 1 in their published work and recalculated in this work via alignment to the asymptotic longrange limit as seen in Figure S21.  Figure S1. Reaction 1 (Ca 2+ + HPO4 2-→ CaHPO4) using the distance between Ca 2+ and the P of HPO4 2and the Ca-OW coordination number as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the calcium water oxygen CN, (c) free energy as a function of both the distance between species and calcium water oxygen CN, and (d) an example geometry at the minimum. Figure S2. Reaction 2 (CaHPO4 + HPO4 2-→ Ca(HPO4)2 2-) using the distance between Ca 2+ , the Ca of CaHPO4 and the Ca-OW coordination number as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the calcium water oxygen CN, (c) free energy as a function of both the distance between species and calcium water oxygen CN, and (d) an example geometry at the minimum. Figure S3. Reaction 3 (Ca(HPO4)2 -2 + Ca +2 → (CaHPO4)2) using the distance between the Ca atoms and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two distinct calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S4. Reaction 4 (CaHPO4 + Ca 2+ → Ca2HPO4 2+ ) using the distance between Ca 2+ and COM of CaHPO4 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two distinct calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S5. Reaction 5 (Ca2HPO4 2+ + HPO4 2-→ (CaHPO4)2) using the distance between Ca 2+ and COM of Ca2HPO4 2+ and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S6. Reaction 6 (CaHPO4 + CaHPO4 → (CaHPO4)2) using the distance between Ca 2+ atoms and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two nominally symmetry equivalent calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S7. Reaction 7 (HPO4 -2 + (CaHPO4)2 → Ca2(HPO4)3 -2 ) using the distance between the COM of HPO4 -2 and the COM of (CaHPO4)2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two nominally symmetry equivalent calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs.     Figure S13. Reaction 13 (Ca2(HPO4)3 -2 + HPO4 -2 → Ca2(HPO4)4 -4 ) using the distance between the COM of Ca2(HPO4)3 -2 and the COM of HPO4 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S14. Reaction 14 (Ca2(HPO4)4 -4 + HPO4 -2 → Ca2(HPO4)5 -6 ) using the distance between the COM of Ca2(HPO4)4 -4 and the COM of HPO4 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the two calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S15. Reaction 15 (Ca(HPO4)3 -4 + Ca(HPO4)2 -2 → Ca2(HPO4)5 -6 ) using the distance between the Ca 2+ of Ca(HPO4)3 -4 and the Ca 2+ of Ca(HPO4)2 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of each of the two distinct calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S16. Reaction 16 (Ca(HPO4)3 -4 + HPO4 -2 → Ca(HPO4)4 -6 ) using the distance between the COM of Ca(HPO4)3 -4 and the COM of HPO4 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of the calcium water oxygen CN, (c) free energy as a function of both the distance between species and the calcium water oxygen CN, and (d) an example geometry at the minimum. Figure S17. Reaction 17 (Ca(HPO4)4 -6 + CaHPO4 → Ca2(HPO4)5 -6 ) using the distance between the Ca 2+ of Ca(HPO4)4 -6 and the Ca 2+ of CaHPO4 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of each of the two distinct calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S18. Reaction 18 (Ca2(HPO4)5 -6 + HPO4 -2 → Ca2(HPO4)6 -8 ) using the distance between the COM of Ca2(HPO4)5 -6 and the COM of HPO4 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of one of the two calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each of the calcium water oxygen CNs. Figure S19. Reaction 19 (Ca(HPO4)3 -4 + Ca(HPO4)3 -4 → Ca2(HPO4)6 -8 ) using the distance between the Ca 2+ of Ca(HPO4)3 -4 and the Ca 2+ of Ca(HPO4)3 -4 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of the distance between species, (b) free energy as a function of either of the two calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S20. Reaction (Ca(HPO4)4 -6 + Ca(HPO4)2 -2 → Ca2(HPO4)6 -8 ) using the distance between the Ca 2+ of Ca(HPO4)4 -6 and the Ca 2+ of Ca(HPO4)2 -2 and the two Ca-OW coordination numbers as the biased CVs. We show (a) the free energy as a function of distance between species, (b) free energy as a function of one of the two distinct calcium water oxygen CNs, (c) an example geometry at the minimum, (d) free energy as a function of both calcium water oxygen CNs, and (e,f) free energy as a function of the distance between species and each one of the two calcium water oxygen CNs. Figure S21. Figure 1 from the main manuscript, but here recreated using the free energy curves reported by Yang et al. 1 plotted as a function of the distance between Ca 2+ and P within a ( ) ( group. The red, green and blue lines are the free energy profiles for ,

Association / Species
] ( + ( à[ ( ) ] ( , respectively. The corresponding dashed lines represent the analytic solution in the long-range limit, which is used to align the absolute free energy. It should be noted that there is significant uncertainty in the alignment of the last of these curves (blue) since the data is only provided out to a shorter distance than would normally be used and the slope has yet to converge to the asymptote.

Supplementary Discussion
On the alignment of free energy profiles As discussed in the Methods section of the main manuscript, the free energy profiles calculated using metadynamics in this work are aligned to an asymptotic limit described by the isolated ions in solution. The pairing free energy for two point particles interacting via a screened electrostatic potential is known analytically: where q is the charge, e0 is the vacuum permittivity, er is the relative permittivity and r is the distance between particles. This can be used to align the free energies obtained from the metadynamics (or umbrella sampling) calculations, whose long-range limit depends on the parameters of the simulations, such as the bias factor or the number of Gaussians that have been deposited. Besides being a useful method of determining convergence of the calculations at the long-range limit, this alignment is necessary to properly normalize the free energies relative to standard state conditions.
The study of Yang et al. 1 does not employ alignment to the asymptotic limit of the free ions in solution in their umbrella sampling free energy profiles. Therefore, the values originally quoted are not at standard conditions and cannot be directly compared to those given in the present manuscript. In order to facilitate a comparison we have attempted to align the data of Yang et al, as shown in Figure S21, resulting in the values quoted in Table S3.