ARTICLE

Received 6 Jul 2011 | Accepted 18 Nov 2011 | Published 20 Dec 2011 DOI: 10.1038/ncomms1604

Stable prenucleation mineral clusters are

liquid-like ionic polymers
Raffaella Demichelis1, Paolo Raiteri1, Julian D. Gale1, David Quigley2 & Denis Gebauer3

Calcium carbonate is an abundant substance that can be created in several mineral forms by
the reaction of dissolved carbon dioxide in water with calcium ions. Through biomineralization,
organisms can harness and control this process to form various functional materials that can
act as anything from shells through to lenses. The early stages of calcium carbonate formation
have recently attracted attention as stable prenucleation clusters have been observed, contrary
to classical models. Here we show, using computer simulations combined with the analysis
of experimental data, that these mineral clusters are made of an ionic polymer, composed of
alternating calcium and carbonate ions, with a dynamic topology consisting of chains, branches
and rings. The existence of a disordered, flexible and strongly hydrated precursor provides
a basis for explaining the formation of other liquid-like amorphous states of calcium carbonate,
in addition to the non-classical behaviour during growth of amorphous calcium carbonate.

1 Department of Chemistry, Nanochemistry Research Institute, Curtin University, PO Box U1987, Perth, Western Australia 6845, Australia. 2 Department of
Physics and Centre for Scientific Computing, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK. 3 Department of Physical Chemistry, University
of Konstanz, Konstanz D-78457, Germany. Correspondence and requests for materials should be addressed to J.D.G. (email: J.Gale@curtin.edu.au).

NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications 1

© 2011 Macmillan Publishers Limited. All rights reserved.
ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604

C
alcium carbonate is a ubiquitous mineral, often created as a therefore explored their impact for the specific case of calcium
result of biomineralization, that is used by nature to perform carbonate. Concentrations between 0.5 and 0.06 M are explored
many diverse functions in marine organisms, from skeletons using extensive simulations of up to 70 ns. Here, the lower bound is
and shells1 to even optical devices2. It also represents a common approximately an order of magnitude greater than the experimental
scale forming from hard water leading to technological problems3, values at which it is possible to maintain solutions containing pre-
though deposition of carbonates can be put to beneficial use as nucleation clusters without them undergoing nucleation. In order
a means of sequestering carbon dioxide4. Despite its economic, to make direct contact with the actual experimental conditions, we
ecological and scientific relevance, our knowledge is far from com- have also performed more limited simulations for systems contain-
plete in regard to the early stages of the complex process that can ing up to 6.4 million atoms that allow the bicarbonate and calcium
ultimately produce this mineral. concentrations to be reduced to 10 and 0.4 mM, respectively. By var-
Recently, there have been several key advances in our under- ying the relative concentrations of carbonate and bicarbonate ions,
standing of the nucleation and growth of calcium carbonate. It has we have also simulated pH values in the range of 8.5–11.5, thereby
been established that formation of amorphous calcium carbonate spanning the experimental conditions.
(ACC) and its subsequent transformation into crystalline phases In the present work, we show that the stable prenucleation clus-
provides a competing pathway to direct creation of the minerals ters of calcium carbonate are in fact ionic polymers consisting of
calcite, aragonite and vaterite5–7. There is growing evidence that chains of cations and anions held together by only ionic interac-
organisms exploit this alternative route for biomineralization as tions. These chains can be linear or branched, with a dynamic
the amorphous phase can be stored until needed, at which point structure that is constantly evolving, yet stable with respect to the
crystallization can be directed to yield the required crystalline separated solvated ions.
polymorph. Furthermore, this process may not follow a conven-
tional mechanism of the type envisaged within classical nucleation Results
theory, in which an activation barrier must be overcome before sig- Atomistic simulation of speciation. Starting from free ions in
nificant association of ions can occur beyond simple ion pairing. solution, the initial stages of speciation involve the expected formation
The existence of stable prenucleation clusters of calcium carbon- of the ion pair, CaCO30, and to a lesser extent, CaHCO3 + (ref 14).
ate before nucleation was initially shown by ion potential meas- At the lower end of the pH range studied, where there is an excess of
urements in combination with analytical ultracentrifugation8, and bicarbonate anions, further association is observed to form species
later confirmed by cryogenic transmission electron microscopy9. of calcium coordinated to two or even three bicarbonate anions at
The size of these clusters has been estimated to be either ~2 nm 0.5 M, though bicarbonate ion pairing becomes unimportant as the
or 0.6–1.1 nm from the two experimental techniques, respectively; concentration decreases. In contrast, calcium and carbonate begin
the discrepancy may be explained depending on the differing sen- to assemble themselves into a new structural motif consisting of
sitivity of the methods to the surrounding hydration layer. While chains of ions resembling a polymer (Fig. 1). These clusters can also
evidence grows for the presence of such stable clusters, the nature include bicarbonate ions at low pH, but normally only at the end
of these species in terms of composition and structure remains of the chain. As time progresses these chains lengthen by further
‘open to speculation’10. collisions with ion pairs and/or ions.
Following the above experimental work, there have been attempts Where the above chains differ from a conventional organic poly-
to use computer simulation to understand the formation of ACC11. mer is in the dynamic nature of their structure. Being held together
Based on a force field model that was accurately calibrated against only by ionic interactions, the chains can frequently break and
experimental free energies of solvation12, we have previously shown reform allowing them to explore a range of configurations involving
that both the free energy of ACC should monotonically decrease as linear chains, rings and branched structures (Fig. 2). Full details of
more ion pairs of calcium carbonate add from solution and that the the coordination numbers of ions within clusters can be found in the
initial association of ion pairs also lowers the free energy, suggest- Supplementary Figures S2–S7, including as a function of concentra-
ing that nucleation of this material may be barrierless and therefore tion and pH. Here, we focus on the high concentration case of 0.5 M
completely non-classical13. This finding unfortunately appears to where the majority of carbonate ions are part of clusters instead of
contradict experiment. These earlier simulations were for solutions ion pairs. In this case, the most common coordination number of
of pure calcium (Ca2 + ) and carbonate (CO32 − ) ions, that is, at high carbonate by calcium is 2, which is consistent with the chain-like
pH, whereas at experimental conditions, the pH is lower and bicar- model. The next most common coordination numbers are 1 and 3,
bonate ions (HCO3 − ) dominate the equilibrium. Experimental corresponding to an ion at the end of the chain and a branching
measurements of the drop in free calcium at the point of nucleation point, respectively. Higher coordination numbers, which are char-
(Δn) show that the magnitude of this event decreases with increas- acteristic of crystalline or ACC, are rarely observed. By forming a
ing pH8, as shown in Supplementary Figure S1, apparently consist- chain, only two waters are removed from the solvation sphere per
ent with the simulations in the limit of high pH. ion, thereby retaining much of the enthalpy of solvation. Further-
Despite the above investigations, the details of the prenucleation more, these clusters have a dynamic structure, which includes the
form of calcium carbonate remain unclear. The challenge is to con- ability to fold and coil like a polymer. This conformational freedom
ceive of a structural form that can exist with a stability intermedi- and the associated entropic contribution to the free energy is the key
ate between ion pairs in solution and the amorphous phase. Any to the stability of these dynamic clusters.
such structure must retain a high degree of hydration otherwise its Formation of chain-like structures is observed in all simulations,
enthalpy would be less exothermic than that of the ion pairs, which regardless of pH and composition. What varies with these condi-
are strongly solvated in water. Conversely, such clusters must also tions is the size distribution, the lifetime of a given length and the
exhibit a high level of disorder to compete entropically with ACC. degree of branching along the chain. At low pH, the chain length
To resolve this issue, we have performed molecular dynamics is limited by competition between carbonate and bicarbonate, with
simulations of ions in solution to probe the details of their asso- the latter species generally acting as a chain terminator as it forms
ciation. While it would normally be highly improbable to observe a weaker link between two calcium ions. At low concentration,
significant levels of aggregation using unbiased simulations, in this growth is limited by the total number of ions available within the
case a direct approach is feasible due to the spontaneous nature of simulation cell. Furthermore, the growth of clusters becomes diffu-
cluster formation. In general, both the concentration of ions and sion limited and so as the concentration decreases, the time taken to
pH can influence the balance of equilibria in solution, and we have grow larger aggregates increases.

2 NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications

© 2011 Macmillan Publishers Limited. All rights reserved.
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604 ARTICLE

a 0.6
Chain pH = 11.5
pH = 10.5
0.5

Conditional probability
pH = 10.0
pH = 9.5
0.4 pH = 8.5

0.3

0.2 Branch

0.1 Terminal

0
0 1 2 3 4 5
Coordination number

Figure 2 | Cluster coordination environment. Probability of finding a

given coordination number of carbonate by calcium as a function of pH,
conditional on the ion being part of a cluster. Points are labelled in the
figure to indicate the underlying cluster structure corresponding to the
Low pH coordination number. Here ‘terminal’ indicates a carbonate either at the
end of a chain or in an ion pair; ‘chain’ represents a carbonate that bridges
two calciums in a chain or ring; ‘branch’ denotes anions where the chain
b bifurcates. Error bars represent the s.d. of the average for successive 1-ns
data windows.

a 0.16
10 ns
0.14 20 ns
0.12 30 ns
40 ns
0.10
P(n)

0.08
0.06
0.04
0.02
0.00
0 10 20 30
n
b 0.07
10 ns
0.06 20 ns
High pH 30 ns
0.05 40 ns

0.04
P(n)

Figure 1 | Species observed in simulations of calcium (bi)carbonate 0.03

solutions. (a,b) Illustration of the species observed at low and high pH,
respectively, with snapshots of the simulation box for 0.5 M at the centre 0.02
(animations of the time–evolution of the simulation are available as 0.01
Supplementary Movies 1 & 2 for pH 10 and 9.5, respectively). Here atoms
are represented as spheres where calcium, carbon of carbonate, carbon 0.00
0 50 100 150 200
of bicarbonate, oxygen and hydrogen are coloured green, blue, yellow, red n
and white, respectively. Ca–C and C–C distances below the cutoffs of 3.9
and 5.0 Å, respectively, are highlighted by purple bonds. For the largest Figure 3 | Cluster size distribution. Probability, P(n), of finding a cluster
cluster (top of b), the oxygen atoms are omitted for clarity so that the containing n ions averaged over successive time intervals of 10 ns for
connectivity is apparent. simulations at a calcium concentration of 0.5 M. Figures show probabilities
for (a) pH 8.5 and (b) pH 11.5. In b, dots represent the true data points while
The differing stability of clusters involving carbonate and bicar- the lines represent smoothed fits to these values as a guide to the eye.
bonate can be determined from the cluster size distribution and its
time evolution (Fig. 3 and also Supplementary Fig. S8 for interme-
diate pH values). In the presence of excess HCO3 − , an exponential size reached in the simulations is a consequence of the number of
decay in cluster size is observed. Here, the simulations are limited by ions present, whereas at lower concentrations the size would be lim-
the amount of carbonate and the timescale for species to encounter ited by the balance between the frequency of collisions between ions
each other in the solution. For carbonate-rich systems, we find a size and the cluster, versus the rate of ion loss.
distribution with multiple peaks that shift to increasing values with
time until a limiting size is reached at higher concentration. This Calcium coordination number. Support for the above struc-
indicates that the clusters we observe are stable, rather than metasta- tural model of a chain-like structure comes from experiment
ble, with respect to solvated ions or ion pairs in solution. The limiting (Supplementary Table S1). By using a speciation model fitted to the

NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications 3

© 2011 Macmillan Publishers Limited. All rights reserved.
ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604

Dry NP
Wet NP
25 DOLLOP

P(Rgyr)
20

15
ΔA (kBT )

10 7 8 9 10 11 12
Rgyr (Å)
5
Figure 5 | Probability distributions as a function of radius of gyration.
0 Probability distributions, P, are compared for dry and wet ACC
4 5 6 7 8 9 10 nanoparticles (NP) containing 36 formula units of CaCO3, as well as
R gyr (Å) for a DOLLOP cluster of the same size. The distributions for the NP are
obtained by direct simulation, while that of DOLLOP comes from the free
Figure 4 | Free energy as a function of radius for a six formula unit energy profile obtained by US. It should be noted that the total probability
cluster. The Helmholtz free energy relative to the lowest energy structure, is unnormalized as the three regions come from separate simulations and
ΔA, is shown as a function of the radius of gyration, Rgyr, based on the cannot be collectively reweighted.
calcium and carbon atoms. Here, the units of energy are Boltzmann’s
constant multiplied by temperature, kBT. Each atom has an ambient
thermal energy of 1.5 in these units. Four sample atomic configurations (up to three formula units) or to both wet and dry ACC nanopar-
are illustrated, along with their first solvation shell, for different radii to ticles led to a very similar free energy profile. As all additions were
demonstrate the wide variety of structures accessible within ambient exothermic, there was no evidence for a barrier to the growth of
thermal energy in the case of the three largest radii. calcium carbonate from ion pairs into an amorphous phase. Here
we have computed the radial free energy profile for a larger 36 for-
mula unit piece of DOLLOP to compare against previous results
experimental data it is possible to compute the average coordina- for an amorphous nanoparticle of the same composition. This gives
tion number of calcium ions that are bound within stable prenu- a similar radial free energy profile as that discussed above for six
cleation clusters. In the pH range of 9.25 to 10.0, we compute that formula units of DOLLOP, though with a narrower radius range.
this coordination number is 2 ± 0.2 (average ± maximum devia- Using data from previous simulations13, we can also determine the
tion). Although this represents an average over the entire cluster, distribution of radii corresponding to ACC, both in its dry state and
the only simple structural models that are consistent with this are when containing structural water. From a comparison of the result-
a ring or long chain. Direct in situ structural information is diffi- ing radial probability distributions, shown in Figure 5, it is clear
cult to obtain. However, ex situ EXAFS studies of rapidly quenched that ACC and DOLLOP represent separate minima on the energy
precursor species also point to twofold coordination of calcium by landscape. To compress the radius of DOLLOP towards that of the
carbonate15. Furthermore, mass spectrometry indicates a series of denser ACC nanoparticles it must undergo significant dehydra-
clusters containing up to eight carbonate ions formed by successive tion. Because of this, there must exist a significant activation barrier
ion pair additions16. between DOLLOP and ACC that prevents spontaneous transitions
from being observed between the two structures. The discrepancy
Free energy sampling. To quantify the thermodynamics of the pol- between the findings of earlier simulations and experiment can
ymer-like chains, we have computed the free energy as a function of therefore be resolved by the observation that DOLLOP and ACC
the gyration radius of the calcium and carbon atoms, an approach have quite distinct densities and levels of hydration. Although the
similar to that used to understand the energetics of biomolecules17. free energy profile for addition of a further ion pair to both cluster
As can be seen from Figure 4, where the free energy landscape for types is similar, it requires a nucleation event to occur for DOLLOP
a six formula unit cluster is shown, this is very revealing. At small to transform to ACC.
radii, corresponding to compressing the cluster towards a bulk-like
structure, the free energy rises steeply. No stable state is found by Speciation analysis of simulation data. A speciation model has
compressing the radius alone, indicating that a significant activation been fitted to the populations observed in the above models, allow-
barrier exists to dehydration. Outside of this, the free energy sur- ing the equilibrium constants and corresponding free energies to
face is incredibly flat and the radius of gyration can be changed by be determined at 298.15 K. Here, we make the same assumption as
almost a factor of two at an energetic cost that is less than ambient used in the experimental analysis that the equilibrium constant for
thermal energy per degree of freedom. This remarkable flexibility is ion pair association is independent of cluster size, as supported by
like distorting the shape of a liquid droplet. For this reason, we have previous simulations13. Three separate fits were performed: (1) a fit
named this stable cluster form of calcium carbonate dynamically to the data for pH 9.5, (2) a fit to the data for pH 10 and (3) a fit to
ordered liquid-like oxyanion polymer (DOLLOP). We have deliber- the complete set of data points for all pH values simultaneously. The
ately not incorporated calcium carbonate into this name as we wish results are contained in Table 1. Although the model contains some
to stress that such species may also be observed for other oxyanions. approximations (primarily that only the most numerous species
Furthermore, given that Pouget et al.9 have demonstrated that such are accounted for and not all possible reactions), three clear points
clusters also exist after nucleation, in addition to being the first emerge. First, the association of ion pairs to form larger clusters has
formed stable species, it seems inappropriate to continue using the an exothermic free energy. Second, the free energies obtained do
phrase ‘prenucleation clusters’. not vary strongly with the pH of the system in the range examined.
In an earlier simulation study13, it was found that the addition Third, the free energies obtained are close to those obtained from
of calcium carbonate ion pairs to either small clusters of ion pairs umbrella sampling (US) techniques thereby providing separate

4 NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications

© 2011 Macmillan Publishers Limited. All rights reserved.
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604 ARTICLE

confirmation that the speciation model is valid. As the statisti- is known to form similar extended chains in both the gas/liquid
cally best speciation data comes from high concentrations, the free phases and in ionic liquids20,21, but the weaker dipolar interactions
energy values relate to these conditions. between molecules are insufficient to preserve these oligomers in
To demonstrate the consequences of the free energies in Table 1, aqueous solution. Arguably a closer analogue to the present system
we have computed the speciation at the experimental conditions, as is that of supramolecular polymers, where typically hydrogen bonds
given in Table 2. Here, the fractions of calcium bound in ion pairs hold the fragments together, rather than covalent linkages22,23.
with both carbonate and bicarbonate are quoted, along with the The key to forming ionic polymers that are stable in solution is the
fraction bound in larger clusters (DOLLOP). The calculations are balance between the solvation free energies of the cation/anion
performed at the same pH values considered experimentally with versus the strength of the cation–anion interaction. Given that other
the same concentration of bicarbonate/carbonate buffer. A linear ions possess similar solvation properties to calcium and carbonate,
relationship is found between the amounts of added and bound cal- there is no reason suppose that this phenomenon is unique to the
cium ions within the low concentration regime, in accord with the present material.
experimental findings. Therefore, we quote specific values only for With any real or virtual experiment, it is important to discuss
one concentration value (0.4 mM) as an illustration. any potential limitations. In the present simulations every care has
The fraction of bound ions in Table 2 can be compared with the been taken to ensure that the thermodynamics of all interactions
experimental values in Supplementary Table S1. This shows that are calibrated against experiment and ab initio quantum mechanics
there is good agreement between the theoretical and experimental where possible. However, no model for such complex systems can
speciation, especially when considering the exponential depend- be perfect. To demonstrate that our observations are not an arte-
ence of the equilibrium constant on any error in the free energy. The fact of the underlying parameterization, we have also performed
largest discrepancy in the model is that the free energy of ion pair- simulations with a different force field that includes reactivity24.
ing between calcium bicarbonate is too exothermic by ~5 kJ mol − 1. Again the formation of flexible chains of ion pairs was found to
To examine the influence of this, the experimental equilibrium con- occur, with the chains being, if anything, more stable than those
stant for this reaction (Kbi) of Plummer and Busenberg18 can be observed with the present force field. Therefore the observation of
substituted into the speciation model. Results for this scenario are stable dynamic chain-like structures is a robust finding, regardless
also given in Table 2. As can be seen, the effect of correcting this of any small quantitative deviations that might exist. On a quantita-
equilibrium constant is to greatly diminish the role of bicarbonate tive level, the finding that the free energy profiles for association of
ion pairing to the point where it is negligible. clusters involve differences at the level of ambient thermal energies
is consistent with the experimental fraction of bound calcium being
Discussion in the range 0.38 to 0.76. A driving force substantially greater than
Association of ionic species in solution has long been known for thermal energy, kBT, would lead to very different ratios of calcium
the case of ion pairs. Even the formation of triple ions is known to concentration.
occur, for example between metal ions and the sulphate anion19. The influence of concentration on the observations also deserves
Beyond this, it is largely assumed that formation of higher oligom- further discussion. Because of the limitations on what is feasible
ers corresponds to the creation of unstable species until the size of with current computers, it is only possible to obtain detailed sta-
the critical nucleus is reached, after which the loss of translational tistical information from simulations at higher concentrations than
entropy of the ions is compensated for by the cohesive energy of those used in the experimental characterization of prenucleation
the proto-solid phase. Our results demonstrate that much larger clusters. To overcome this, a speciation model has been fitted to the
stable chain-structured oligomers are possible for calcium carbon- simulation results at high concentrations. The equilibrium constants
ate. Unlike conventional organic polymers, these chains are held obtained can then be used to extrapolate to the lower concentra-
together by the ionic bond between a calcium cation and carbon- tions employed experimentally. The theoretical speciation model is
ate anion rather than a covalent linkage, something that is argu- built on the assumption that DOLLOP grows by adding successive
ably unprecedented for this type of material. Hydrogen fluoride ion pairs with a constant exothermic free energy. At low concentra-
tions, this yields an approximately linear relationship between the
calcium added and the number of bound metal ions, which agrees
Table 1 | Thermodynamics of speciation. almost quantitatively with the experimental observations, thereby
supporting the validity of this extrapolation.
Reaction pH 9.5 pH 10 Fit to A further argument relating to the influence of concentrations
only only all pHs
is that direct formation of a crystalline phase by a classical pathway
Ca2 + + HCO3 − ⇒CaHCO3 + − 11.4 − 11.4 − 11.3 may occur under the conditions simulated that bypasses the aggre-
Ca2 + + CO32 − ⇒CaCO30 − 18.6 − 20.5 − 20.3 gation of ions to give stable DOLLOP clusters. As we are presently
CaCO30 + (CaCO30)n⇒(CaCO30)n + 1 − 22.9 − 22.0 − 21.7
unable to simulate the nucleation event, we cannot comment defini-
Free energies (kJ mol − 1) for three reactions obtained by fitting a speciation model to the tively on what happens at this stage. However, it is possible to argue
simulation data at three concentrations of 0.5, 0.28 and 0.06 M.
that prenucleation clusters of the type described are extremely likely

Table 2 | Speciation of aqueous calcium carbonate/bicarbonate.

pH 9.00 pH 9.25 pH 9.50 pH 9.75 pH 10.0

Fraction of bound Ca2 + in CaCO30 ion pairs for Kbitheory 0.157 0.196 0.205 0.209 0.212
Fraction of bound Ca2 + in DOLLOP for Kbitheory 0.417 0.581 0.683 0.731 0.755
Fraction of bound Ca2 + in CaHCO3 + ion pairs for Kbitheory 0.202 0.104 0.051 0.026 0.013
Fraction of bound Ca2 + in CaCO30 ion pairs for Kbiexp 0.187 0.202 0.208 0.209 0.212
Fraction of bound Ca2 + in DOLLOP for Kbiexp 0.484 0.649 0.721 0.753 0.766
Fraction of bound Ca2 + in CaHCO3 + ion pairs for Kbiexp 0.041 0.018 0.008 0.004 0.002
Speciation is computed as a function of pH using the free energies from Table 1 based on the fit to all pH values. Calculations are performed for a carbonate/bicarbonate buffer concentration of 10 mM.
Values are quoted for a calcium concentration of 0.4 mM, as an example, given that the bound fraction is linear over the experimental range.

NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications 5

© 2011 Macmillan Publishers Limited. All rights reserved.
ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604

to maximize the distance between them. During this simulation, the

same variety of structures are explored as for simulations under dif-
ferent conditions; linear and branched chains are observed, while
ions occasionally become solvent separated by one or more shells,
though often return to the cluster again, as can be seen in Figure 6.
Hence, this behaviour is qualitatively no different from that observed
for the same cluster either at high concentrations or when simulated
in water alone (Supplementary Movie 3). Changing the concentra-
tion alters the frequency of collisions between ions/ion pairs and
clusters, and so the maximum size will vary according to the bal-
ance of association/dissociation events. Indeed, the lifetime of a
particular cluster becomes longer, on average, as the concentration
decreases due to the reduced frequency of collision. Despite this
variation in cluster size and lifetime with conditions, the underlying
structural characteristics of DOLLOP remain the same.
In the process of biomineralization, organisms control the nucle-
ation and growth of calcium carbonate to create complex-functional
materials to serve a variety of purposes. While it still remains prob-
able that much of this control is exerted through the ACC phase
into which DOLLOP can ultimately transform by partial dehydra-
tion, there is also the possibility that this new form may also have
a role under particular conditions. For example, it has been pro-
posed26 that certain peptides may bind chains of calcium carbon-
ate, significantly altering the structure of the resulting material. Our
Figure 6 | Structures of the four formula unit prenucleation clusters. work indicates that these chains are likely to be longer than the 3–4
Structures shown represent the configurations of four separate clusters ions previously speculated to exist as a cluster26, and is consistent
after 1 ns of simulation at experimental conditions ([Ca] = 0.4 mM, with the observation of a long-range effect of the peptide. Other
[HCO3 − ] = 10 mM, pH = 10). Atoms are coloured green, blue and red for simulation studies have indicated that organic molecules (specifi-
calcium, carbon and oxygen, respectively. Note that very similar structures cally those containing Asp residues in this case) can have a different
can be seen in the Supplementary Movie 3 where four ion pairs are started influence by modifying the solvation environment of Ca2 + even at
in a solvent-separated state in aqueous solution. Surrounding water and a separation of ~4.9 Å27.
remote sodium, bicarbonate and carbonate species are omitted for clarity. Liquid-like phases of calcium carbonate have also been prepared
by addition of acidic polypeptides to generate polymer-induced liq-
uid precursor28. Recent work has demonstrated that liquid emul-
to still exist. Under genuinely homogeneous conditions (that is, sions can be formed without the use of organic additives if the
where there is no pre-existing seed or interface to grow from), ions materials are levitated to prevent nucleation16,29, while addition of
have to assemble in solution before the nucleation of the crystalline proteins can either stabilize or destabilize this liquid ACC30. The
phase. Previous simulations have already shown that calcite nano- concept of a dynamic polymer structure involving highly hydrated
particles are only the most stable configuration beyond sizes of a calcium carbonate provides a unifying framework that connects all
few nanometres13. Hence, it is hard to explain why ions would form of these observations, while the existence of liquid precursors for
metastable ordered structures when there is a stable disordered con- materials beyond calcium carbonate indicates that this may well
figuration. From an entropic perspective, the direct creation of an be a more general phenomenon. Finally, manipulation of this inor-
ordered structure with a low volume in configuration space when ganic polymer offers new pathways by which organic molecules can
there is an ensemble of structures with a large volume and lower influence biomineralization and biomimetic synthesis at the earliest
free energy is improbable. The key point to emphasize here is that stages of crystal growth31,32.
we make no claim that the present stable clusters described lie on
all nucleation pathways that might occur under different condi- Methods
tions. However, they are very likely to be present during the early Force field. Here, existing parameters13 for aqueous calcium carbonate systems
were extended to include bicarbonate and sodium species. The intramolecular
stages of speciation under most homogenous conditions, though carbonate parameters were also refined to accurately describe the vibrational
either nucleation to ACC or dissolution and growth of a competing spectrum of this anion33, as described in Supplementary Table S2, using the
species may subsequently occur. relaxed fitting algorithm within GULP34. As for carbonate, the water–bicarbonate
Finally on the issue of concentration, we should consider the parameters were fitted to the experimental free energy of solvation. Given there are
results obtained for the largest simulations, performed at the experi- limited experimental data to fit for bicarbonate, we have used quantum mechanical
information to supplement this at the M06/6-31 + G** to M06/6-311 + + G** level35
mental conditions, in comparison with those conducted at higher in NWChem36. Quantum mechanical information is primarily used to fit repulsive
concentrations. As simulations at this scale will be limited by the interactions and to validate the energies/structures of complexes observed during
rate of diffusion and collision frequency of ions, it is not possible molecular dynamics. Final parameters for bicarbonate inter- and intramolecular
to observe the spontaneous formation of clusters in a statistically interactions are given in Supplementary Tables S3 and S4, respectively.
meaningful way. However, by placing clusters taken from simula- As association between calcium and carbonate/bicarbonate anions is an
important feature of the speciation, we have computed the binding free energies of
tions at higher concentrations into this simulation cell it is possible the ion pairs by US, yielding energies of − 27 and − 10.5 kJ mol − 1 for CaCO30 and
to determine whether there is any significant change in the structure CaHCO3 + , respectively. When corrected to standard concentrations, these values
and lifetime of the clusters at this lower concentration. As the times- are shifted to approximately − 22 and − 5.5 kJ mol − 1. The corresponding experi-
cale for water exchange in the calcium coordination environment mental values are − 18 and − 6.5 kJ mol − 1 (ref. 18), though the value for CaCO30
is ~100 ps25, several nanoseconds of simulation is sufficient for a includes the contribution of prenucleation clusters.
To model the experimental conditions, it is also necessary to include Na +
cluster to undergo a significant level of rearrangement and even dis- whose interaction with carbonate was fitted to natrite (Lv/298 K)37 under the
sociate completely if unstable. Four clusters containing four formula approximations of full sodium occupancy and no modulation. The parameters
units of calcium carbonate were placed in the simulation cell so as were tested for eitelite38, with magnesium–carbonate parameters obtained from

6 NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications

© 2011 Macmillan Publishers Limited. All rights reserved.
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604 ARTICLE
magnesite. Here, the largest error in the cell parameters was only 0.77%. The an ion pair adds. Recognizing the potential influence of bicarbonate ion pairing
sodium–bicarbonate interactions were obtained by scaling the Buckingham A on the equilibria, we also include this as a competing process with a separate equi-
parameter for the carbonate case by the ratio of oxygen charges. These values were librium constant. A least squares fit was then performed to the cluster populations
tested for sodium bicarbonate where errors in cell parameters were between 1.13 obtained from the simulations.
and 2.24%. Finally, the sodium–water interaction was taken from Grossfield et al.39
and checked to ensure that the change of water model has negligible effect. Experiment. The experimental details are described elsewhere for potential meas-
urements, titration and speciation8, as well as for structural analyses of quenched
Atomistic molecular dynamics simulation. Molecular dynamics simulations precursors15. In this work, data from these studies have been reanalysed and used
were performed for solutions containing calcium, carbonate and bicarbonate ions for discussion. Specifically, the multiple-binding equilibrium introduced in previ-
under conditions of ambient temperature (300 K) and pressure (1 atmosphere). ous work8 has been examined in more depth. Evaluation of this physicochemical
Concentrations of 0.5, 0.28 and 0.06 M were examined based on the molarity of model gives the parameters, x and K, which are pH dependent. Here, x represents
the calcium ions, while the ratio of carbonate to bicarbonate was varied to yield the number of binding sites for calcium ions on a carbonate ion, and the cor-
conditions equivalent to pH values of 8.5, 9.5, 10, 10.5 and 11.5 based on the responding equilibrium constant is given by K. As shown by titration and calcium
experimental equilibrium constant between the species. Cubic simulation cells were potential measurements8, equal amounts of calcium and carbonate ions are bound
used of initial dimension ~10 nm that were varied within an NPT ensemble using a in prenucleation clusters. Hence, it is reasonable to assume that calcium ions have
timestep of 1 fs. Typically this volume contains more than 30,000 water molecules. the same number of binding sites as carbonate ions, at the respective pH values.
Single ions are initially distributed randomly throughout the solution while avoiding In solution, during the prenucleation stage, there are free ions and ions bound in
close contacts between atoms. All simulations were run for between 50 and 70 ns clusters, and the number of binding sites is an averaged value that includes those
of real time using LAMMPS40. A chain of five Nosé–Hoover thermostats was ions with empty binding sites, that is, ions that are free in solution. In order to
employed with thermostat and barostat relaxation times of 0.1 and 1 ps, respectively. derive the macroscopic coordination number, Ncalcium, in prenucleation clusters
The electrostatics were treated using the Particle–Particle Particle–Mesh Ewald from x, the fraction of bound ions, f, has to be taken into account;
algorithm with a precision of 10 − 4.
Additional simulations were performed using a larger cubic simulation cell of nbound (ca2 + )
f = ,
side ~40 nm, containing just over 2.1 million water molecules for a bicarbonate nbound (ca2 + ) + nfree (ca2 + )
concentration of 10 mM and a carbonate/bicarbonate ratio consistent with a pH of
10. To match the experimental conditions just before the nucleation event at this
pH, 16 Ca ions were placed in this cell to give a concentration of ~4×10 − 4 M, while where nbound(Ca2 + ) and nfree(Ca2 + ) are the amounts of bound and free calcium,
charge balance was maintained through the addition of sodium cations, as per the respectively. Ncalcium, which is valid only for ions bound in clusters, is then given by:
real system. Owing to the magnitude, these simulations were run for 1 ns. As simu- x
lations at this scale will be limited by the rate of diffusion and collision frequency of N calcium = .
ions, it is not possible to observe the spontaneous formation of clusters. However, f
by placing clusters taken from simulations at higher concentrations into this simu-
lation cell it is possible to determine whether there is any significant change in the
structure and lifetime of the clusters at this lower concentration. References
US simulations41 were performed using the PLUMED plug-in42 in the NVT 1. Belcher, A. M. et al. Control of crystal phase switching and orientation by
ensemble on a single cluster within a cubic cell containing ~6,300 water molecules soluble mollusc-shell proteins. Nature 381, 56–58 (1996).
(6 formula unit cluster) or 15,200 water molecules (36 formula unit cluster). Here, 2. Aizenberg, J., Tkachenko, A., Weiner, S., Addadi, L. & Hendler, G. Calcitic
the radius of gyration was employed as the collective variable for US. The calcula- microlenses as part of the photoreceptor system in brittlestars. Nature 412,
tion for the 6 formula unit DOLLOP was performed with 16 windows, each run 819–822 (2001).
for 20 ns including equilibration, equally spaced in the interval 4.2–9.9 Å using a 3. Chen, T., Neville, A. & Yuan, M. Calcium carbonate scale formation—assessing
spring constant of 0.95 eV/Å. For the 36 formula unit case the interval spanned the initial stages of precipitation and deposition. J. Petrol. Sci. Eng. 46, 185–194
was 6.6–14.5 Å and the spring constant used was 0.5 eV/Å. The strength of each (2005).
restraining potential was chosen so that the equilibrium probability distribution, 4. Matter, J. M. & Kelemen, P. B. Permanent storage of carbon dioxide in
P(Rgyr), in each window overlapped with its neighbour at no more than 1.2 s.d.’s geological reservoirs by mineral carbonation. Nat. Geosci. 2, 837–841
from the mean. Expansion by dissociation of the cluster was allowed, however, to (2009).
suppress periodicity artefacts, the maximum distance between any two solute ions 5. Beniash, E., Aizenberg, J., Addadi, L. & Weiner, S. Amorphous calcium
was prevented from exceeding half of the simulation cell by a smooth and continu- carbonate transforms into calcite during sea urchin larval spicule growth.
ous repulsive function. Probability distributions were computed only in regions Proc. Roy. Soc. Lond. 264B, 461–465 (1997).
where this repulsive function was negligible. The free energy was constructed using 6. Radha, A. V., Forbes, T. Z., Killian, C. E., Gilbert, P. U. P. A. & Navrotsky,
weighted histogram analysis within the WHAM software43. A. Transformation and crystallization energetics of synthetic and biogenic
amorphous calcium carbonate. Proc. Natl Acad. Sci. USA 107, 16438–16443
Analysis of simulation data. Analysis of the unbiased speciation simulations is (2010).
based on snapshots of the entire system taken at 5 ps intervals, and of the solute 7. Rodriguez-Blanco, J. D., Shaw, S. & Benning, L. G. The kinetics and
ions at intervals of 0.1 ps. The analysis of clustering was based on a geometric crite- mechanisms of amorphous calcium carbonate (ACC) crystallization to calcite,
rion for defining connections between ions. This required the distance connecting via vaterite. Nanoscale 3, 265–271 (2011).
centres of two ions (taken as the coordinates of the carbon for molecular anions) 8. Gebauer, D., Völkel, A. & Cölfen, H. Stable prenucleation calcium carbonate
to be shorter than 3.9 Å for Ca–C and 5.0 Å for C–C. These values were chosen to clusters. Science 322, 1819–1822 (2008).
be close to the transition state separating directly coordinated species from those 9. Pouget, E. M. et al. The initial stages of template-controlled CaCO3 formation
that are solvent separated. A recursive search over all connections was used to revealed by Cryo-TEM. Science 323, 1455–1458 (2009).
identify discrete instantaneous clusters. Cluster size statistics were accumulated 10. Meldrum, F. C. & Sears, R. P. Now you see them. Science 322, 1802–1803
over time windows of 1 ns to eliminate noise arising from sporadic ion–cluster (2008).
and cluster–cluster proximity. 11. Tribello, G. A., Bruneval, F., Liew, C. & Parrinello, M. A molecular dynamics
Each cluster of non-trivial size was characterized by the number of ions with study of the early stages of calcium carbonate growth. J. Phys. Chem. B 113,
each coordination number in the range 1–5 (higher values not being observed). 11680–11687 (2009).
These quantities were averaged over the lifetime of the cluster. The cluster charge 12. Raiteri, P., Gale, J. D., Quigley, D. & Rodger, P. M. Derivation of an accurate
and composition were also recorded. Here, we adopted a more stringent cluster force-field for simulating the growth of calcium carbonate from aqueous
definition, requiring persistence with the same ions for >2 ps. Correspondingly, solution: a new model for the calcite-water interface. J. Phys. Chem. C 114,
any cluster that breaks and re-forms within a time of 2 ps was considered to have 5997–6010 (2010).
persisted over that interval. Hence, the total number of ions assigned to these 13. Raiteri, P. & Gale, J. D. Water is the key to non-classical nucleation of
persistent clusters is a conserved quantity. At each frame we recorded the number amorphous calcium carbonate. J. Am. Chem. Soc. 132, 17623–17634
and maximum size of both persistent and instantaneous clusters. We also recorded (2010).
the number of ions of each type with zero connectivity, allowing direct connection 14. Greenwald, I. The dissociation of calcium and magnesium carbonates and
to the experimental free ion concentration. bicarbonates. J. Biol. Chem. 141, 789–796 (1941).
The cluster size distribution information was used to perform a simple specia- 15. Gebauer, D. et al. Proto-calcite and proto-vaterite in amorphous calcium
tion analysis to provide a validation of the free energies for ion association. Here, carbonates. Angew. Chem. Int. Ed. 49, 8889–8891 (2010).
we assume that equilibria exist between calcium and carbonate ions to form ion 16. Wolf, S. E. et al. Carbonate-coordinated metal complexes precede the formation
pairs, which are then able to associate further to form clusters with two separate of liquid amorphous mineral emulsions of divalent metal carbonates. Nanoscale
equilibrium constants (that is, one for ion pairing and one for addition of an ion 3, 1158–1165 (2011).
pair to a cluster). In the spirit of the experimental speciation model8, we assume 17. Kohn, J. E. et al. Random-coil behavior and the dimensions of chemically
that the equilibrium constant does not depend on the size of the cluster to which unfolded proteins. Proc. Natl Acad. Sci. USA 101, 12491–12496 (2004).

NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications 7

© 2011 Macmillan Publishers Limited. All rights reserved.
ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms1604

18. Plummer, L. N. & Busenberg, E. The solubilities of calcite, aragonite and 36. Bylaska, E. J. et al. NWChem, A Computational Chemistry Package for Parallel
vaterite in CO2-H2O solutions between 0 and 90 °C, and an evaluation of Computers, Version 5.1.1 (Pacific Northwest National Laboratory, 2009).
the aqueous model for the system CaCO3-CO2-H2O. Geochim. Cosmochim. 37. Arakcheeva, A. et al. The incommensurately modulated structures of natural
Acta 46, 1011–1040 (1982). natrite at 120 and 293 K from synchrotron X-ray data. Am. Miner. 95, 574–581
19. Schröder, D. et al. Direct observation of triple ions in aqueous solutions of (2010).
nickel(II) sulfate: a molecular link between the gas phase and bulk behavior. 38. Knobloch, D., Pertlik, F. & Zemann, J. Crystal structure refinements of
J. Am. Chem. Soc. 133, 2444–2451 (2011). buetschlite and eitelite: a contribution to the stereochemistry of trigonal
20. Redington, R. L. Infrared absorbance of hydrogen fluoride oligomers. carbonate minerals. Neues Jahrb. Mineral. 230–236 (1980).
J. Phys. Chem. 86, 561–563 (1982). 39. Grossfield, A., Ren, P. & Ponder, J. W. Ion solvation thermodynamics from
21. Enomoto, T., Nakamori, Y., Matsumoto, K. & Hagiwara, R. Ion-ion interactions simulation with a polarisable force field. J. Am. Chem. Soc. 125, 15671–15682
and conduction mechanisms of highly conductive fluorohydrogenate ionic (2003).
liquids. J. Phys. Chem. C 115, 4324–4332 (2011). 40. Plimpton, S. J. Fast parallel algorithms for short-range molecule dynamics.
22. Lehn, J.- M. Supramolecular polymer chemistry—scope and perspectives. J. Comp. Phys. 117, 1–19 (1995).
Polymer Intl. 51, 825–839 (2002). 41. Torrie, G. M. & Valleau, J. P. Non-physical sampling distribution in Monte
23. de Greef, T. F. A. & Meijer, E. W. Supramolecular polymers. Nature 453, Carlo free-energy estimation: umbrella sampling. J. Comp. Phys. 23, 187–199
171–173 (2008). (1977).
24. Gale, J. D., Raiteri, P. & van Duin, A. C. T. A reactive force field for aqueous- 42. Bonomi, M. et al. PLUMED: a portable plugin for free-energy calculations
calcium carbonate systems. Phys. Chem. Chem. Phys. 13, 16666–16679 (2011). with molecular dynamics. Comput. Phys. Commun. 180, 1961–1972
25. Ohtaki, H. & Radnai, T. Structure and dynamics of hydrated ions. Chem. Rev. (2009).
93, 1157–1204 (1993). 43. Grossfield, A. WHAM: the Weighted Histogram Analysis Method, version
26. Metzler, R. A., Tribello, G. A., Parrinello, M. & Gilbert, P. U. P. A. Asprich 2.04,http://membrane.urmc.rochester.edu/content/wham.
peptides are occluded in calcite and permanently disorder biomineral crystals.
J. Am. Chem. Soc. 132, 11585–11591 (2010).
27. Hamm, L. M., Wallace, A. F. & Dove, P. M. Molecular dynamics of ion Acknowledgements
hydration in the presence of small carboxylated molecules and implications for This work was supported by Australian Research Council grant DP986999 and EPSRC
calcification. J. Phys. Chem. B 114, 10488–10495 (2010). grant EP/H00341X/1. We thank the Pawsey Centre, iVEC, NCI and HECToR for
28. Gower, L. B. & Odom, D. J. Deposition of calcium carbonate films by a computing facilities. D.G. thanks Markus Antonietti for discussions regarding the
polymer-induced liquid-precursor (PILP) process. J. Crystal Growth 210, structure of prenucleation clusters.
719–734 (2000).
29. Olszta, M. J., Odom, D. J., Douglas, E. P. & Gower, L. B. A new paradigm for
biomineral formation: mineralization via an amorphous liquid-precursor.
Author contributions
R.D., P.R., D.Q. and J.D.G. performed and analysed the simulations. D.G. conducted the
Connect. Tissue Res. 44, 326–334 (2003).
experimental work and analysis. All authors contributed to the writing and editing of the
30. Wolf, S. E., Leiterer, J., Pipich, V., Barrea, R., Emmerling, F. & Tremel, W. Strong
manuscript.
stabilization of amorphous calcium carbonate emulsion by ovalbumin: gaining
insight into the mechanism of ‘polymer-induced liquid precursor’ processes.
J. Am. Chem. Soc. 133, 12642–12649 (2011). Additional information
31. Meldrum, F. C. & Cölfen, H. Controlling mineral morphologies and structures Supplementary Information accompanies this paper at http://www.nature.com/
in biological and synthetic systems. Chem. Rev. 108, 4332–4432 (2008). naturecommunications
32. Sommerdijk, N. A. J. M. & de With, G. Biomimetic CaCO3 mineralization
using designer molecules and interfaces. Chem. Rev. 108, 4499–4550 (2008). Competing financial interests: The authors declare no competing financial interests.
33. Coleyshaw, E. E., Crump, G. & Griffith, W. P. Vibrational spectra of the Reprints and permission information is available online at http://npg.nature.com/
hydrated carbonate minerals ikaite, monohydrocalcite, lansfordite and reprintsandpermissions/
nesquehonite. Spectrochim. Acta A 59, 2231–2239 (2003).
34. Gale, J. D. GULP: capabilities and prospects. Z. Krist. 220, 552–554 (2005). How to cite this article: Demichelis, R. et al. Stable prenucleation mineral clusters are
35. Zhao, Y. & Truhlar, D. G. The M06 suite of density functionals for main group liquid-like ionic polymers. Nat. Commun. 2:590 doi: 10.1038/ncomms1604 (2011).
thermochemistry, thermochemical kinetics, noncovalent interactions, excited
states, and transition elements: two new functionals and systematic testing of License: This work is licensed under a Creative Commons Attribution-NonCommercial-
four M06-class functionals and 12 other functionals. Theoret. Chem. Acc. 120, Share Alike 3.0 Unported License. To view a copy of this license, visit http://
215–241 (2008). creativecommons.org/licenses/by-nc-sa/3.0/

8 NATURE COMMUNICATIONS | 2:590 | DOI: 10.1038/ncomms1604 | www.nature.com/naturecommunications

© 2011 Macmillan Publishers Limited. All rights reserved.