Phase-Field Modeling of Biomineralization in Mollusks and Corals: Microstructure vs Formation Mechanism.

While biological crystallization processes have been studied on the microscale extensively, there is a general lack of models addressing the mesoscale aspects of such phenomena. In this work, we investigate whether the phase-field theory developed in materials' science for describing complex polycrystalline structures on the mesoscale can be meaningfully adapted to model crystallization in biological systems. We demonstrate the abilities of the phase-field technique by modeling a range of microstructures observed in mollusk shells and coral skeletons, including granular, prismatic, sheet/columnar nacre, and sprinkled spherulitic structures. We also compare two possible micromechanisms of calcification: the classical route, via ion-by-ion addition from a fluid state, and a nonclassical route, crystallization of an amorphous precursor deposited at the solidification front. We show that with an appropriate choice of the model parameters, microstructures similar to those found in biomineralized systems can be obtained along both routes, though the time-scale of the nonclassical route appears to be more realistic. The resemblance of the simulated and natural biominerals suggests that, underneath the immense biological complexity observed in living organisms, the underlying design principles for biological structures may be understood with simple math and simulated by phase-field theory.

Introduction

Crystalline materials formed by solidification from the liquid state play an essential role in our civilization.[1,2] This class of matter incorporates most of the technical alloys, polymers, minerals, drugs, food products, and so on. Owing to their importance, mathematical models describing the process of crystallization under the respective conditions were and are being developed. Relying on the statistical physical description of phase transitions, the evolving numerical methods, and the ever-increasing computational power, computational materials science reached the level where knowledge-based design of crystalline matter is possible for certain classes of materials (see, e.g., refs (2−4)). The models that address the behavior of matter during crystalline solidification range from the molecular time and length scales to the engineering scales. They include ab initio computations; particle-based methods like molecular dynamics (MD), Monte Carlo, or population dynamics simulations and different types of continuum models ranging from the density functional theory of classical particles, via coarse-grained models (such as the time-dependent Ginzburg–Landau, Cahn–Hilliard, and phase-field type order parameter theories that belong to the family of classical field theoretical models widely used in modeling phase transitions of various complexity), to the macroscopic continuum models applicable on engineering time- and length-scales. While this inventory allows the modeling of a substantial range of crystallization phenomena, there are complex cases, for which its use is not straightforward. Such examples are the biomorphic (inorganic) materials[5−10] that form worm-shape or arboresque morphologies by aggregation of crystalline particles, and the process of biomineralization;[11−19] that is, the formation of hierarchically structured organic–inorganic composites in biological systems. Examples of biomineralization include the formation of mollusk shells,[13,14] skeletons of corals[15] and cell walls of diatoms,[16] kidney stones,[17] bones and teeth,[18] and magnetite crystals in the magnetosomes of magnetotactic bacteria,[19] to name a few. The materials formed by biomineralization often have surprisingly good mechanical properties owing to their hierarchical microstructure (see e.g., refs (13,14)). Recent imaging and analytic methods provide detailed information on the respective microstructures, which in turn may give clues to the formation mechanism: many of these microstructures are well-known from materials science (such as dendrites, spherulites, cellular, columnar shapes, etc.).[11−15,17−19] This raises the possibility that with some adjustment/further development, the models developed in materials science can be used to reverse engineer the biomineralization process and learn the pathways used by nature to create these complex structures, which may inspire new technologies for creating novel composite materials.[20−25]

Recently, we explored the possibility of developing predictive mathematical models for biomorphic crystallization and for relatively simple biomineralization processes by adopting well-established methods of computational materials science and adjusting them to the circumstances as necessary.[26−28] The research done so far is confined yet to relatively simple cases of extracellular biomineralization such as mollusk shell formation[26,27] or microstructure evolution of spherulitic structures in coral skeletons[28] but is expected to deepen the general understanding in the field, and the tools developed in the course of this research might open the way for modeling more complex cases of crystallization in biological systems such as formation of bones, kidney stones, and so on.

In the present paper, we concentrate on the modeling aspects of such an approach, outlining possible minimum requirements for phase-field modeling of biological crystallization processes, and demonstrate that with appropriate choice of the model parameters and boundary conditions phase-field models can approximate the polycrystalline microstructure formed in simple cases of biomineralization (shell formation in mollusks such as bivalves, gastropods and cephalopods, and sprinkle formation in coral skeletons).

Microstructures Formed during Biomineralization

Before outlining the phase-field models we used in the present research, we give a short account of the experimental results on the observed microstructures. Polycrystalline microstructures formed in biomineralization processes have been investigated by a variety of experimental methods, including optical microscopy (OM), scanning electron microscopy (SEM), electron back scattering diffraction (EBSD), X-ray tomography, and polarization-dependent imaging contrast (PIC,[29,30]), among others. Here we concentrate on experimental results obtained on the microstructure of molluscan shells and coral skeletons.

Mollusk Shells

Mollusk shells are complex organo-mineral biocomposites with a broad range of species-dependent microstructures.[13,14,30−39] A schematic view of bivalve anatomy having a nacre-prismatic shell is shown in Figure . Moving inward, the sequence of the individual layers is as follows: the organic periostracum, a leathery “skin”, that encloses the domain where biominerlization takes place, a dominantly mineral (calcium carbonate, CC) prismatic layer, and the nacre composed of CC tablets and organic interlamellar membranes, the submicron thick extrapallial liquid,[31,34−36] and then the outer calcifying epithelum layer of the mantle. Images showing typical microstructures of mollusk shells of similar type are displayed in Figure .[13,26,37,38]

Figure 1

Schematic view of bivalve molluscan anatomy with successive magnification of the mantle–nacre interface.[39] A thin liquid-filled extrapallial space is indicated (its thickness and content is open to debate). The interlamellar membrane is made of a viscoelastic chitin-based organic substance, whereas the mineral constituent is crystalline CC (aragonite). (Reproduced with permission from ref (39). Copyright 2009 United States National Academy of Sciences). Here “growth direction” indicates the lateral growth of the shell. However, this is based on layerwise growth perpendicular to the surface via formation of new layers that also grow laterally, a process which is responsible for both the thickening and the sidewise spreading of the shell. In the present work, we model the local thickening of the shell.

Figure 2

Some typical mircrostructures observed in mollusk shells: (a) cross-sectional EBSD orientation map for the shell of Katelysia rhytiphora,[37] the outer side is on the left (note the transitions between layers of different crystallite morphologies); (b) outer randomly oriented granular domain in the shell of Unio pictorum (SEM);[26] (c) columnar prismatic domain of Pinna nobilis (X-ray tomography reconstruction),[13] and (d) plate-like structure of the nacre of Unio pictorum (SEM);[26] (e) spherulitic layer, section perpendicular to growth in the shell of Haliotis rufescens (SEM).[38] [(a) Reproduced from ref (37) under Creative Commons Attribution 4.0 License. Copyright Authors 2019. (b), (d) Reproduced with permission from ref (26) . Copyright 2018 Wiley. (c) Reproduced with permission from ref (13) . Copyright 2014 Springer Nature. (e) Reproduced from ref (38) . Copyright 2002 American Chemical Society.].

A recent study shows that in members of three classes of mollusks Unio pictorum (bivalve), Nautilus pompilius (cephalopod), and Haliotis asinine (gastropod), the shell displays a common sequence of ultrastructures: a granular domain composed of randomly oriented crystallites, a prismatic domain of columnar crystallites, and the nacre[26,27] (Figure ). It has been shown that the layered structure of nacre may contain screw dislocation-like defects (see Figure ).[14,39−44]

Figure 3

Hierarchy of ultrastructures in the shell of some of the mollusks. Left: for bivalve Unio pictorum (growth direction: upward); (A) The shell, (B) schematic drawing of microstructure, (C–F) scanning electronmicroscopy (SEM) images. Right top: for cephalopod Nautilus pompilius (growth direction: to the right). Right bottom: for gastropod Haliotis asinina (growth direction: to the right); (A–C) SEM image, (C,D) magnified views of the granular domain and the prismatic → nacre transition. Note the similar sequence of ultrastructures during growth: granular → prismatic (columnar) → nacre (alternating mineral/organic layers). (Left: Reproduced with permission from ref (26). Copyright 2018 Wiley. Right: Reproduced with permission from ref (27). Copyright 2019 United States National Academy of Sciences).

Figure 4

(a–h) Screw dislocation-like defects at the growth front of the nacreous layer of various species.[14,39] (a) Aragonitic nacre structure in the shell of the bivalve Pteria avicula. Reproduced with permission from ref (40). Copyright 2007 The Royal Society (U.K.). (b) Calcitic seminacre structure in the shell of brachiopod Novocrania anomala. Reproduced with permission from ref (41). Copyright 2010 The Paleontological Association. (c) Calcitic foliated structure in the shell of Ostrea edulis. Reproduced with permission from ref (42). Copyright 2007 Elsevier. (d) Aragonitic foliated structure in the shell of monoplacophoran Rokopella euglypta. Reproduced with permission from ref (43). Copyright 2009 American Chemical Society. Scale bars indicate 5 μm. SEM images showing the growth front of the nacreous layer for bivalves: (e–g) Pteria avicula and (h) Pteria hirundo. Reproduced with permission ref (40). Copyright 2007 The Royal Society (U.K.). (i,j) Target and spiral patterns formed in Belousov–Zhabotinsky reaction. Reproduced with permission from ref (44). Copyright 1996 Royal Society of Chemistry.

Of these structures, the prismatic layers show mechanical flexibility, whereas the nacre (also called “mother of pearl”) is fairly rigid but hard; the combination of the two yields a surprisingly strong yet flexible biocomposite (see e.g., ref (22)).

Microstructure of Coral Skeletons

The multiply branched shapes of coral skeletons are covered by a large number of coral polyps[45] and the connecting living tissue, which secretes calcium carbonate to create a hard shelter (the corallite, a tubular hollow structure on which the polyp sits),[46] into which the polyp can retreat if danger is detected (see Figure a,b). The polyps are transparent, their color originates from photosynthesizing algae (zooxanthallea) that live in symbiosis with the polyp and feed the polyp sugars and oxygen. The surface of the skeleton is intricately structured,[46] depressions, ridges, cavities are arranged into complex patterns reflecting the radial symmetry of the polyp (Figure b). The CC crystal (aragonite) building the porous skeleton (Figure c) has a spherulitic microstructure, that is, a radial arrangement of crystals radiating from a common center. The center does not have to be a point; it can be a line or a plane in plumose spherulites. In the case of coral skeletons, the centers are curved planes, termed centers of calcification (CoCs) (Figure d).[15,47] From these CoCs, acicular fibers grow radially and then arrange into fan-like bundles that finally group into a feather-duster-like shape termed “trabecula”.[15] Besides the plumose spherulitic structure, randomly oriented nanoscale crystallites “sprinkles” are also present[28,48,49] (Figure e). Recent PIC mapping experiments performed at synchrotron indicate that the amount of the submicron-sized sprinkle crystallites varies considerably among different coral species: in some of them, the sprinkles are missing (Figure f), whereas in others, submicron size crystallites appear at the perimeter of the skeleton, along grain boundaries, at the growth front of spherulitic trabeculae, or in bands formed in the interior of the skeleton.[28] The discovery of sprinkles inspired a refined model for spherulitic growth in corals, described in detail in ref (28). Randomly oriented sprinkles are the first nucleated crystals at the growth front. With further growth, those oriented radially have space and thus continue to grow, and those oriented tangentially run into each other and stop growing. This is why they stay small. A coarsening process then makes the larger crystal grow larger at the expense of the smaller ones, which disappear. In most mature spherulites, therefore, no sprinkles remain. In the skeleton of some coral species, however, some sprinkles do not disappear, presumably because they are kinetically trapped.

Figure 5

From coral polyps to the CC (aragonite) skeleton they form, highly structured from the cm to the nm scale. (a) Schematic drawing of the coral polyp (colored) sitting on the porous skeleton (white).[45] Reproduced with permission from Britannica eReader.com, a service of Encyclopædia Britannica, Inc. (b) SEM image of the surface formed below coral polyps (corallites) in the case of Porites sp. Reproduced with permission from ref (46). Copyright 2011 Elsevier. (c) Cross-section of corallite, Pollicipora damicornis (linear size of image is ∼3 mm, “d” stands for mineral bridges termed distally convex dissepiments). Reproduced with permission from ref (47). Copyright 2006 Springer Nature. (d) Schematic drawing of a “plumose” spherulite. Reproduced with permission from ref (15). Copyright 2017 American Chemical Society. (e) PIC map of Acropora pharaonis coral skeleton showing spherulitic microstructure of crystals radiating from a band of randomly oriented sprinkles (linear center of calcification, the bar indicates 10 μm). Reproduced with permission from ref (28) under Creative Commons CC-BY-NC-ND 4.0 license. Copyright 2020 Elsevier. (f) PIC map of Phyllangia coral skeleton showing spherulitic microstructure of crystals radiating from centers of calcification, however, with no sprinkles.

In our previous work, we modeled this process by phase-field simulations[28] and raised the possibility that other spherulites may grow this way, including aspirin, chocolate, and geologic crystals. However, the formation mechanism of sprinkle bands and the origin of different amount of sprinkles in the skeleton of different coral species is not yet fully understood. Mathematical modeling is expected to help to identify the governing factors. Molecular/ionic mobility in the calcifying fluid is expected to be orders of magnitude higher than in the solid. As a result, in the case of growth via ion-by-ion addition, either a slow supply of the ions or a kinetic barrier of ion deposition can keep the growth rate sufficiently low. This offers limitations to the possible mechanisms, as will be discussed later.

Biomineralization on the Nano- and Macroscale

The complexity of biomineralization stems mainly from the fact that the fluid and solids incorporate organic molecules, the role of which is largely unknown.[35,36,50,51]

For example, nanoscale amorphous calcium carbonate (ACC) globules are essential for the formation of mollusk shells and coral skeletons.[13,14,52−54] Possible pathways from free ions to crystalline CC are reviewed in ref (55). Evidently, all the possible processes cannot be explicitly incorporated into an orientation field based phase-field (OF-PF) approach designed to model polycrystalline microstructures on the mesoscale. In this work, we investigate two specific cases: (i) formation of crystalline CC via the classical route of ion-by-ion addition of Ca2+ and CO32–, and (ii) the nonclassical route via an amorphous precursor (ACC) deposited at the solidification front. In the latter case, we hypothesize that the crystallization rate is determined by the velocity of crystal growth into the ACC layer and not by the rate of supplementing ACC; that is, growth is controlled by the self-diffusion in ACC. This hypothesis can account for the typically months’ time scale of shell/skeleton formation, which would be difficult to interpret, for example, on the basis of ion-by-ion deposition directly from the extrapallial fluid or other aqueous solutions. A specific realization of mechanism (ii) is presented in ref (56). The extrapallial fluid contains various ions and organic molecules as shown by in vivo studies.[35,36]

Growth Rate of Mollusk Shells and Coral Skeletons

The shell of bivalves grows typically by 100–300 μm lunar-day increments,[57,58] corresponding to a growth rate of about v ≈ 4.2 × 10–11–1.3 × 10–10 m/s, which decreases with age.[57] The thickening rate of the shell of Tridacna deresa was estimated to be v ≈ 1.6 × 10–10–4.9 × 10–10 m/s in its early life, which decreases to v ≈ 3.2 × 10–11–2.3 × 10–10 m/s in the later life.[59] Comparable growth rates were reported for freshwater gastropods v ≈ 9.5 × 10–11 m/s.[60] In contrast, the coral skeleton growth rates range between about 1 and 37 cm/year, corresponding to 3.2 × 10–10 to 1.2 × 10–8 m/s.[61−63]

Modeling Section

There are two main categories of the phase-field (PF) models developed to address polycrystalline freezing: (a) the multiphase-field (MPF) models that assign a separate phase field for every crystal grain[64−70] and (b) the orientation-field based phase-field (OF-PF) approaches, in which the local crystallographic orientation is monitored by a scalar field (2D)[71−78] or quaternion/rotation matrix fields (3D).[78−85] Recent developments in these areas were reviewed in refs (70) and (78), respectively.

Both the MPF and OF-PF approaches have their advantages, yet complex polycrystalline growth forms (such as disordered dendrites, crystal sheaves, various types of spherulites, and fractal-like polycrystalline aggregates, etc.) were so far modeled exclusively using the OF-PF approach. A further advantage of these models is that they allow a continuous variation of the orientation, a feature particularly useful in modeling biomineralization.

In the OF-PF models, the local phase state of the matter is characterized by a coarse grained structural order parameter, the “phase field”, that is time- and space-dependent and monitors the transition between the liquid and solid states. This field is usually coupled to other slowly changing fields, such as the concentration field of the constituents, and the orientation field. The free energy contains bulk free energy and gradient terms for these fields. The equations of motion can be obtained via a variational route. The thermal fluctuations are represented by adding noise (nonconserved fields) or flux noise (conserved fields) to the equations of motion that satisfy the fluctuation–dissipation theorem, allowing for homogeneous nucleation[72]

The inclusion of foreign surfaces/particles, represented by appropriate boundary conditions that set the wetting properties, allows the modeling of heterogeneous nucleation[83−85] and solidification in confined space.[77,86] The addition of phase-field noise makes the modeling of Brownian motion of solid particles possible.[87] The OF-PF models may also be coupled to hydrodynamic flow, which can be represented by either the Navier–Stokes equation[88,89] or the lattice Boltzmann technique. With an appropriate overlapping grid technique, flow coupled motion of growing solid particles can also be treated.[90]

The main virtue of PF modeling is that the growth morphology can be computed on the basis of the freezing conditions and the thermophysical properties. Images obtained by OF-PF modeling in two- and three dimensions (2D and 3D)[74,78,84,86] are compared with experiments[74,86,91−101] in Figure .

Figure 6

Complex crystallization structures in the experiments (1st and 3rd rows) and in the corresponding simulations performed using OF-PF models at the Wigner Research Centre for Physics (2nd and 4th rows). Upper block: From left to right: Rhombic dodecahedron crystals (upper panel reproduced from ref (91), lower panel reproduced with permission from ref (84). Copyright 2008 IOP Publishing Ltd.). Columnar dendrites (upper panel reproduced with permission from ref (92). Copyright 2016 Elsevier). Equiaxed dendrites (upper panel reproduced with permission from ref (93). Copyright 2016 Elsevier, lower panel reproduced with permission from ref (78). Copyright 2013 Springer Nature). “Dizzy” dendrites (reproduced with permission from ref (74). Copyright 2003 Springer Nature). Crystal sheaves (upper panel reproduced with permission from ref (94). Copyright 1993 American Chemical Society, lower panel reproduced with permission from ref (78). Copyright 2013 Springer Nature). Orientation field in spherulites (upper panel reproduced with permission from ref (95). Copyright 2007 American Chemical Society). Arboresque spherulites (upper panel reproduced with permission from ref (96). Copyright 2003 American Chemical Society, lower panel reproduced with permission from ref (78). Copyright 2013 Springer Nature). Lower block: From left to right: Scratch induced dendritic crystallization in polymer film (reproduced with permission from ref (86). Copyright 2003 Springer Nature). Dumbbell-shape spherulites (upper panel reproduced with permission from ref (97). Copyright 2006 Wiley, lower panel reproduced with permission from ref (84). Copyright 2008 IOP Publishing Ltd.). Spherulites in temperature gradient (upper panel reproduced with permission from ref (98). Copyright 2003 Springer Nature, lower panel reproduced with permission from and ref (86). Copyright 2013 Springer Nature). Freezing soap bubble vs dendrite growing in spherical shell (upper panel reproduced with permission from ref (99)). Effect of oscillating temperature (reproduced with permission from ref (78). Copyright 2013 Springer Nature). “Quadrites” formed by nearly 90° branching (upper panel by courtesy of B. Okerberg,[100] lower panel reproduced with permission from ref (78). Copyright 2013 Springer Nature). “One eyed” spherulites (upper panel reproduced with permission from ref (101). Copyright 2002 Elsevier).

Recently, PF modeling has been extended for microstructure formation in mollusk shells and coral skeletons, and promising agreement was seen between experiment and the predicted microstructures.[26−28] We give below a detailed account of the modeling efforts and present new results.

Herein, we continue further the quest for a minimum phase-field model of the biomineralization process during the formation of mollusk shells and coral skeletons. Evidently, we cannot model the living organism in the framework of this approach; their functions are represented by appropriate boundary conditions. Unfortunately, usually little is known of the thermodynamics of the multicomponent fluids involved, of the interfacial free energies and their anisotropies, and the respective diffusion coefficients. Therefore, we aim at identifying the main components that a minimal theory needs to incorporate for qualitatively reproducing the microstructures seen in the experiments.

In the present study, we rely on three specific formulations of the phase-field theory. In two of them,[26,27] the local state is characterized by the phase field, ϕ(r, t) that monitors the process of crystallization, the concentration field c(r, t) specifying the local composition, and a scalar orientation field θ (r, t), which represents the local crystallographic orientation in 2D. The latter field is made to fluctuate in time and space in the liquid, a feature that represents the short-range order present in the liquid state (as done in refs (72,74−79) and (83−86)). The third model does not include an orientation field. It has been developed to handle two-phase spiraling dendrites during eutectic solidification in ternary alloys in 3D,[102] assuming a fixed relative orientation of the solid phases.

In Figure , we summarize the problems addressed herein, the models used, the assumed micromechanisms, and the corresponding phase transitions.

Figure 7

Summary of biomineralization-related problems that are addressed within this work using phase-field methods. In each case, the applied model, the respective micromechanism, and the relevant phase transition are specified. Phase-field models PF1, PF2, and PF3 are defined in the text.

Phase-Field Model 1 (PF1)

The first model will be used to address four cases of biomineralization: (i) mollusk shell formation from aqueous solution by ion-by-ion addition (ions from the extrapallial fluid attach to the surface of the crystalline phase), (ii) mollusk shell formation via amorphous precursor (ACC → CCC transition), (iii) formation of columnar nacre via ion-by-ion attachment, and (iv) formation of coral skeletons via ion-by-ion attachment. In all these cases, the same mathematical model will be used, however, with specific input parameters and initial- and boundary conditions.

Model PF1 is defined by eqs −6 shown below. This model is similar to the standard binary PF theory by Warren and Boettinger;[103] however, it is supplemented with an orientation field as done in refs (72,74−78) Accordingly, the free energy of the heterogeneous system is expressed aswhere T is the temperature. Parameters ε2 = (12/√2) γiδi /Ti and w(c) = (1 – c) wA + c wB, where w = (12/√2) γi/(δiTi) are expressed in terms of the free energy γi, the thickness δi of the crystal–liquid interface, and the melting point Ti of the ith pure component (i = A or B that stand for the organic and mineral components, respectively). In model PF1 ε2 = 0 is chosen. s = s(ϑ, θ) = 1 + s0 cos{kϑ – 2πθ} is an anisotropy function corresponding to an interfacial free energy of k-fold symmetry and strength s0, whereas ϑ = arctan(ϕ/ϕ) is the angle of the normal of the interface in the laboratory frame, while ∇ϕ = [ϕ, ϕ]. The angular (circular) variables ϑ and θ are normalized so that they vary between 0 and 1. The bulk free energy density reads as,and varies between the free energy densities of the crystal and mother phases (fC and fM, respectively) as prescribed by the interpolation function p(ϕ) = ϕ3 (10–15ϕ + 6ϕ2). Here fC and fM were taken from the ideal solution model. The orientation free energy density is as follows:where the parameter H can be used to tune the magnitude of the grain boundary energy. Owing to the scalar orientation field, model PF1 is applicable exclusively in 2D.

The time evolution of the heterogeneous system is described by variational equations of motion (EOMs):where M, M, and M are mobilities that determine the time scale of the evolution of the individual fields, and are related to coefficients of the self-diffusion, interdiffusion, and rotational diffusion.[72,74−78] The chemical and orientation mobilities are made phase-dependent as M = M + [1 − p(ϕ)]{MM – MC}, where i = c or θ, and indices M and C denote values for the mother and crystalline phases. The corresponding dimensionless mobilities are defined as m = MεT/D, m = M/D, where M= (vm/RT)c(1 – c)D, and m = MξHT/D. Here vm is the average molar volume of the components, R the gas constant, ξ the length scale, whereas HT is the energy scale of the grain boundary energy. Gaussian white noise terms ζ are added to the EOMs to represent the thermal fluctuations (here i = ϕ, c, and θ). 3D generalizations of model PF1 can be found elsewhere,[78−82] which, however, require quaternion or rotation matrix representation of the crystallographic orientation, as opposed to the scalar field used here.

Phase-Field Model 2 (PF2)

The second model will be used to provide a refined model of mollusk shell formation including the nacreous structures in the case of (i) ion-by-ion attachment and (ii) amorphous precursor mediated process. This OF-PF model was originally developed to describe eutectic solidification, while keeping a fixed orientational relationship between the two solid phases inside the crystal grains.[104] To realize this, the square-gradient term, (1/2)εc2T (∇c)2, was retained in the free energy density (choosing εc2 = 2ε2), and a more complex orientational free energy term was usedthat realizes a fixed orientational relationship at the solid–solid phase boundaries. Here h(c) = (1/2){1 + cos[2π (c – cα)/(cβ – cα)]}, cα and cβ are the CC concentrations in the two solid solution phases, whereas F1(|∇θ|) = |∇θ| and F2(|∇θ|) = a + b|cos(2mπd|∇θ|)|. Here a, b, m, are constants, and d is the characteristic thickness of the sold-solid phase boundary. The EOMs are derived the same way as in the case of model PF1. Accordingly, model PF2 is defined by eqs , 2, 4–6, and 7.

Phase-Field Model 3 (PF3)

This model will be used to address the formation of 3D topological defects in sheet nacre via crystal growth into hACC (here the third component plays the role of water that has smaller solubility in CCC than hACC). Modeling of sheet nacre by PF3 is a 3D analogue of the amorphous precursor mediated case in PF2. The free energy of the crystallizing system is given by the expression[102]where c = [c1,c2,c3] and Σc = 1, whereas w and ε2 are constants. The bulk free energies of the solid and liquid phases are taken from the regular and ideal solution model:andwhere Ω are the binary interaction coefficients in the solid.

The respective EOMs obtained variationally are as followsandwhere Mc is the 3 × 3 mobility matrix. With the choice of 1 and −0.5 for the diagonal and off-diagonal elements, the criterion Σc = 1 is automatically satisfied. The diffusion is switched off in the bulk solid. Further information on model PF3 is given in ref (102).

Numerical Solutions

The EOMs were solved numerically in a dimensionless form, on rectangular uniform grids, using finite difference discretization with forward Euler time stepping. The PF1 and PF2 codes were run on a CPU cluster of 480 CPU cores using MPI protocol. Typical runs on a 1000 × 2000 grid took between about 8 to 15 h, depending on the number of time steps that varied from 2.4 × 105 to 5 × 105, as required by the velocity of crystallization. The code for PF3 was run on high-end graphics processing units (GPUs) and was solved on a 3D rectangular grid.

Materials’ Parameters

We review here the present status of input data required for a quantitative modeling. The OF-PF models require a fairly detailed information on the systems studied. This incorporates the free energy of all the relevant phases as a function of temperature and composition; all the interface energies; and the translational, chemical, and rotational diffusion coefficients. Since in the biomineralization problems, we address here, the dominant CC polymorph is the metastable aragonite, we use the thermophysical data available for this polymorph, as much as possible. In case, where no information is available, we use values for another CC polymorph (calcite). The input data are collected in Tables and 2 for Model PF1, in Tables and 4 for Model PF2, a few common ones are presented in Table . In these Tables, “aq. sol. → CCC” indicates data relevant to the ion-by-ion mechanism, whereas “ACC → CCC” denote those for the amorphous precursor mediated case.

Table 1

Dimensionless Mobility Coefficients for Model PF1a

	mφ	mc,M	mc,C	mθ,M	mθ,C
aq. sol. → CCC	3.75	1.0	10–20	120	120 × 10–20
ACC → CCC	3.75	1.0	10–14	120	120 × 10–14

The subscripts M and C stand for the mother and crystalline phases. The chemical mobility of the former was used as reference, as its chemical diffusion coefficient was used in making the EOMs dimensionless.

Table 2

Dimensionless Thermodynamic Data Used in Model PF1 (Ideal Solution Thermodynamics[26])a

quantity		value
Tr	= T/TA	0.911
Tr,B	= TB/TA	0.786
ΔgA	= ΔGA/RT	–0.1184
ΔgB	= ΔGB/RT	0.1554

Here ΔGA,B = ΔHA,B(T–TA,B)/TA,B (A stands for CC and B for the organic component).

Table 3

Dimensionless Mobility Coefficients for Model PF2a

	mφ	mc,M	mc,C	mθ,M	mθ,C
aq. sol. → CCC	0.0144	1.0	10–20	12	12 × 10–20
ACC → CCC	0.0144	1.0	10–14	12	12 × 10–14

The subscripts M and C stand for the mother and crystalline phases. The chemical mobility of the former was used as reference, as its chemical diffusion coefficient was used in making the EOMs dimensionless.

Table 4

Dimensionless Thermodynamic Data Used in Model PF2 (Regular Solution Thermodynamics[27])a

quantity		value
Tr	= T/TE	0.720
Tr,A	= TA/TE	1.169
Tr,B	= TB/TE	1.286
ΔgA	= ΔGA/RT	– 0.5802
ΔgB	= ΔGB/RT	– 1.0477
ωM	= ΩM/RT	2.0510
ωC	= ΩC/RT	3.6335

Here ΩM,C = Ω0,M,C – Ω1,M,CT.

Table 5

Materials and Computational Data Used in PF1/PF2

quantity	value	unit	ref
γA (CCC – aq. sol.)	150	mJ/m2	(123,124)
γB (organic–aq. sol.)	118	mJ/m2	this work
γA (CCC – ACC)	87	mJ/m2	this work
γB (organic – ACC)	68	mJ/m2	this work
vm (CCC – aq. sol.)	26.7	cm3/mol
vm (CCC – ACC)	32.4	cm3/mol
ξ	2.1 × 10–6	m
δ	4.15 × 10–8	m
Δx	6.25 × 10–3	dimensionless
Δt	4.75 × 10–6	dimensionless

Thermodynamics

Unfortunately, only limited thermophysical information is available even for the pure CC system from experiment and MD simulations, such as the phase diagrams[105,106] and the equilibrium shapes reflecting the anisotropy of the interface energy.[107,108] During biomineralization, however, a variety of ions and organic macromolecules are present that may influence/control the crystallization process.[109−115] Accordingly, it is a nontrivial task to obtain accurate input data for mesoscale modeling; for example, selective adsorption of ions or organic molecules on different crystal faces may change growth morphology[111,112] or influence the formation of polymorphs of CC.[114,115]

Owing to this lack of information, we present here generic approaches that are based on simplified hypothetical model systems of properties similar to those used in refs (26−28).

Diffusion Coefficients

As noted above, we address here two scenarios for the formation of crystalline CC (CCC): (i) diffusion controlled growth of crystalline CC directly from aqueous solution via ion-by-ion addition; (ii) diffusion controlled growth of crystalline CC into a hydrated ACC (hACC) layer that is assumed to be deposited on the solidification front (by vesicles or the-ion-by-ion process) with a sufficient rate so that deposition is not the rate-limiting process. These scenarios differ in the diffusion coefficient we assign to the “mother phase” (either aqueous solution or ACC) that crystallizes. Since the equations of motion were made dimensionless using the diffusion coefficient of the mother phase, and we assume that the relative magnitudes of M, M, and M remain the same in the mother phase independently whether it is liquid or amorphous, the two scenarios differ in only the dimensionless mobilities assigned to the crystalline phase.

Aqueous Solutions

The coefficient of ion diffusion in aqueous solutions at room temperature is in the order of DL ≈10–9 m2/s.[116]

Amorphous CC

The ion diffusion in ACC at 300 K is in the order of DACC,ion ≈ 10–15m2/s.[117] However, the diffusion coefficient of the water molecules from MD simulations is typically DhACC,H ≈ 10–14–10–13m2/s for the slow H2O molecules, although a few percent of H2O molecules that have orders of magnitude faster diffusion (DhACC,H ≈ 10–11 m2/s) are also present.[117,118] However, biogenic ACC is almost anhydrous.[119] Therefore, water diffusion is expected to play a negligible role. The rate limiting factor for the structural transition is expected to be the slowest of these processes; accordingly, we use the diffusion coefficient for the ions, DACC,ion ≈ 10–15m2/s.[117]

Crystalline CC

We are unaware of self-diffusion data for aragonite. There are, however, experimental data for calcite. The Mg diffusion coefficient in calcite is about Dcalcite ≈ 10–21 m2/s at 823 K, which, extrapolates to Dcalcite ≈ 10–53 m2/s at room temperature provided that the diffusion mechanism does not change.[120] A different estimate is obtained via extrapolating the Mg–Ca interdiffusion data of ref (121) to room temperature, which yields Dcalcite ≈ 10–29 m2/s at 300 K for Mg diffusion in calcite. An even higher value emerges from the radioactive tracer method Dcalcite ≈ 10–23 m2/s.[122]

Summarizing, herein, we opt for DACC,ion ≈ 10–15m2/s and Dcalcite ≈ 10–29 m2/s for ion diffusion in ACC and calcite, respectively; assuming thus that the diffusion data for calcite can be viewed as a reasonable order-of-magnitude estimate for the other polymorphs, including aragonite. The corresponding dimensionless mobility data are presented in Tables and 3. The data of magnitude 10–20 to 10–12 indicate that in the crystal, independently of the mother phase, the reduced mobilities are essentially zero for the concentration and orientation fields; that is, these fields do not show time evolution within the crystal on the time scale of the simulations.

Interfacial Free Energies

The experimental and theoretical results for the water-aragonite interfacial free energy are about 150 mJ/m2.[123,124] A recent ab initio theoretical treatment provides a considerably larger value (280 mJ/m2), and information on its anisotropy for small aragonite clusters.[125] Herein, we use 150 mJ/m2. We are unaware of data for the free energy of the ACC-aragonite interface. Using Turnbull’s relationship[126] for the interfacial free energy, γ = αΔH/(Nvmc2)1/3, where α is a constant, ΔH heat of transformation, N0 the Avogadro-number, and vmc is the molar volume of the crystalline phase, a crude estimate can be made on the basis of the enthalpy difference between ACC and calcite: ΔHcalcite-ACC = (14.3 ± 1.0) kJ/mol.[119] A similar value may be expected for aragonite.[119] Considering α = 0.55 from MD simulations,[127] one obtains γaragonite-ACC ≈ 87 mJ/m2, which result, however, needs independent confirmation by other experimental/theoretical methods.

Once the thermodynamic data are fixed for components A and B, and the interfacial free energy is given for one of the components, models PF1 and PF2 predict the interfacial free energy for the other component, provided that the interface thicknesses are similar.[103] For materials of comparable entropy of transformation, this realizes γ ∝ Ttrans, where Ttrans is the temperature of the phase transition, a relationship that works well for the solid–liquid interfacial free energy of metals.[128]

Qualitative Modeling

Despite our efforts to collect a full set of the required materials parameters, owing to uncertainties of the thermodynamic diving force of crystallization and of the interface energy estimates, the simulations we present can only be regarded as qualitative. They are aimed at demonstrating that phase-field modeling has the potential to capture various microstructural/morphological aspects of biomineralization. This summary of the present status of input data may give hints where further experiments and microscopic theory can help mesoscale modeling.

Results and Discussion

Modeling of Microstructures Mimicking Mollusk Shells

Herein, we hypothesize that (1) the granular domain is produced by heterogeneous nucleation (either on the surface of periostracum or on organic particles); and (2) the physical background of the columnar → nacre morphological transition is the observation that decreasing the driving force of solidification, an initially diffusionless process (full solute trapping) is replaced by partitioning, which first appears in the form of alternating layers rich in one or the other component (see Figure ).[129,130]

Figure 8

Band formation as predicted by the binary Warren-Boettinger[103] (WB) type phase-field models: (a) transition from diffusionless solidification to partitioning with decreasing driving force (supersaturation decreases from left to right) in the presence of a single solid phase (ideal solution); and (b) transition from lamellar eutectics toward bands parallel to the growth front with increasing driving force (undercooling) in a WB model supplemented with square gradient term for concentration,[130] a case in which two solids of different composition form (regular solution). Lower panel was reproduced with permission from ref (130). Copyright 2017 Springer.

The latter mechanism is the one that causes band formation in Figures and 10. Further assumptions made during the application of these models (e.g., boundary conditions specific to the studied problems) are recapitulated below. In the case of PF1, reduction of the driving force leads to a transition of a chemically homogeneous solid to alternating solid–liquid layers, whereas in the case of PF2 partitioning appears via alternating α–β bands. For even smaller driving forces, PF1 would produce seaweed/dendritic structures, whereas PF2 would yield lamellae perpendicular to the growth front.

Figure 9

Comparison of the microstructure of (a) the shell of mollusk Unio pictorum(26) as shown by electron microscopy images (see also Figure left block; reproduced with permission from ref (26). Copyright 2018 Wiley) with simulations (b–e) obtained by model PF1 assuming (b, c) ion-by-ion attachment from extrapallial fluid (D = 10–9 m2/s) or (d, e) crystallization from an ACC layer (D = 10–15 m2/s). Growth direction is upward. Orientation (b, d) and grain boundary (c, e) maps are shown. In panels (b) and (d), colors denote different crystallographic orientations. In the experimental images of (a), the bars correspond to 3, 3, 1, and 1 μm, respectively, from top to bottom. Note the presence of the three characteristic domains in the experiment and in both types of simulations: granular, columnar prismatic, and sheet nacre structures. The qualitative phase-field simulations were performed on 1000 × 2000 grids (corresponding to 13.125 μm × 26.25 μm with the present choice of model parameters). (In the high driving force (lower) part of both these simulations lDv /D > 1; that is, diffusionless crystallization takes place, whereas in the upper domain, alternating mineral and organic layers mimicking the sheet nacre form.).

Figure 10

Three stages of microstructure evolution in Model PF2 obtained assuming (a–d) ion-by-ion attachment from extrapallial fluid and (e–h) by crystallization from ACC precursor(orientation (a–c) and (e–g), and composition maps (d) and (h) are shown.): (a,e) Formation of granular structure via heterogeneous nucleation dominated equiaxed solidification; (b, f) columnar growth via directional solidification in concentration gradient yielding the prismatic structure, and (c, g) layerwise formation of alternating CCC and organic layers (sheet nacre). In panels (a–c) the mother phase is an aqueous solution (D = 10–9 m2/s) in the simulation shown; whereas an amorphous precursor (D = 10–15 m2/s) is assumed in the simulation shown in panels (d–f). The wavelength of the layered structure is roughly proportional to the free energy of the mother phase – CCC interface. The thickness of the mother phase (extrapallial domain) is assumed to be constant. In (a–c) and (e–g) different colors stand for different crystallographic orientation, while colors white and black stand for the mantle of the mollusk and the mother phase, respectively. In (d) and (h), gray, yellow, and blue indicate the mother phase, the organic phase, and CCC. The qualitative phase-field simulations were performed on 2000 × 200 grids (corresponding to an area of 26.25 μm × 2.625 μm with the present choice of model parameters). In (a–c) and (e–g) time elapses downward, whereas (c) and (d) and (g) and (h) display snapshots taken at the same time.

Shell-Like Microstructure in Model PF1

In a recent OF-PF study,[26] we made the following assumptions that define the conditions under which eqs –6 were solved, when modeling the formation of mollusk shells within model PF1:

The CC crystals grow into the extrapallial fluid by the molecule/ion attachment mechanism.

Binary ideal solution thermodynamics (CC and organic component) is applied. Evidently, treating the extrapallial fluid as a quasi-binary solution is a gross simplification. During crystallization of the CC-rich crystal, and an organic-component-rich “fluid” forms from the original homogeneous mixture. This construction was used as a simple means to provide thermodynamic driving force for CCC precipitation.

CC-supersaturation of the extrapallial fluid decreases exponentially with the distance x from the periostracum, owing to a spatially dependent amount of the organic component: c(x) = cmin + (cmax – cmin){1–exp(− 9x/L)}, where L is the thickness of the extrapallial space. (This is a hypothesis. We are unaware of any experimental information pro or contra.)

Crystallization of CC starts via heterogeneous nucleation on the periostracum.

The anisotropy of the CCC-mother phase interfacial free energy is neglected.

In the present work, besides this, we explore a different scenario shown in Figure , in which the CCC crystal grows into an ACC precursor that forms continuously ahead of the crystallization front. The CCC front propagates into this ACC layer and has no direct contact with the extrapallial fluid. It is also assumed that the formation of the ACC layer is fast enough, so that it is not the rate-limiting process. The respective amorphous → crystal transition has in principle a reduced driving force, while orders of magnitude smaller diffusion coefficients prevail in the amorphous phase, yielding a far longer time scale for crystallization, when compared with crystallization from an aqueous solution. In this scenario, the phase field monitors the amorphous → crystal transition, rather than the crystallization of a liquid. We furthermore assume that the coefficients of the translational, chemical, and rotational diffusion decrease proportionally by orders of magnitude during this process, retaining the same relative magnitudes of the mobilities as in the liquid (see Table ).

Since the EOMs are solved in dimensionless form, where time is dedimensionalized using the chemical diffusion coefficient of the mother phase, the dimensionless chemical mobilities of the mother phase remain unchanged. What differs between the present computation and the previous one in ref (26) is the magnitude of the individual mobilities in the crystalline and the mother phase (see Table ). Following the general principles of phase-field modeling (see e.g. ref (131)), the phase-field mobility is assumed to be independent of the phase. In contrast, in agreement with the experimental diffusion coefficients, the chemical and orientational mobilities are assumed to be 20 orders of magnitude larger in the extrapallial fluid then in the crystal, and 14 orders of magnitude larger in ACC than in CCC.

We retain the assumptions made in ref (26), as listed above with the difference that in the nonclassical mechanism it is the CC content of the ACC layer that decreases exponentially (while the organic content increases) with the distance from the periostracum. For the sake of simplicity, we employ the same dimensionless driving force as in ref (26) for both the classical and the nonclassical cases.

The microstructures that evolved in the two cases are compared in Figure . The characteristic microstructural transitions are present in both simulations. Whether from the fluid or the amorphous phase, first small randomly oriented CCC grains form in the neighborhood of the periostracum, of which crystal grains grow further inward yielding elongated crystals of random orientation that compete with each other. With increasing distance from the periostracum, that is, with decreasing supersaturation, growth slows down, and the separation of the two constituents becomes possible, which leads to the formation of alternating CCC and organic-rich layers, closely resembling the nacre. The mineral layers in the “nacre” are composed of segments of different orientations (see Figure b,d), which can be viewed as a 2D analogue of the usual 3D mineral platelets. The predicted sequence of the morphological transitions is similar for both mechanisms (i.e., for aq. sol. → CCC and ACC → CCC); however, there are minor differences in the relative thicknesses of the granular, columnar prismatic, and layered nacre structures.

While the respective microstructures are rather similar, the typical size scales for the ACC → CCC transition is smaller than for the aq. sol. → CCC, roughly proportionally with the respective interfacial free energies. Alternating CCC and organic-rich layers akin to sheet nacre form here due to a process described in ref (129), with the difference that under the present conditions a roughly flat “banded structure” forms, and thermal diffusion is replaced by chemical diffusion. Apparently, the orientational information is only partly transferred through the organic layers. Mineral bridges (discontinuities of the organic layers) are also observed. The granular → columnar → layered morphological transitions occur here because of changes in the growth velocity.

At high supersaturations nucleation and diffusionless solidification takes place forming the granular domain. At medium supersaturations nucleation ceases, only competing growth of the existing particles takes place yielding the columnar domain, whereas at small supersaturations alternating CCC and organic layers occur, forming a layered structure that closely resembles the sheet nacre.

Diffusionless crystallization is possible, when the diffusion length lD = D /v is comparable to the thickness of liquid/crystal or amorphous/crystal interfaces (d ≈ 10–9–10–8 m[132,133]), where v is the growth rate. The transition from diffusion controlled to diffusionless growth takes place in the regimes 10–4 < dv/D < 1 or 10–2 < dv/D < 10, depending on the model.[134] Considering lD ≈ 10–9 m and a typical experimental growth rate of v ≈ 10–10 m/s, one finds that the mother phase needs to have a diffusion coefficient of D ≈ 10–15–10–20 m2/s to show transition toward diffusionless growth on the time scale required. This clearly rules out the possibility that the CC crystals form dominantly by direct ion-by-ion attachment from the extrapallial fluid as D ≈ 10–9 m2/s applies for the latter process. In turn, this range of D is consistent with crystallization from amorphous CC, as the magnitude of D falls in the range diffusion coefficient takes in the amorphous state.[135]

Summarizing, while the two mechanisms considered here lead to similar microstructures, crystallization via the ACC precursor seems preferable to direct solidification via ion-by-ion addition from the extrapallial fluid, as in the latter case the diffusion coefficients are rather high, that is, crystallization is expected to be fast, unless the “fluid” is of high viscosity (not realistic for the extrapallial fluid). Another problem of direct precipitation from the aqueous solution is that due to the high diffusion coefficient, the assumed initial exponential spatial dependence of CC supersaturation is only temporary on the time scale of shell growth, unless crystal growth is so fast that diffusional equilibration cannot take place. This is, however, at odds with the experimental growth rates.

Note that modeling of the experimentally observed orientational ordering in the columnar (prismatic) and nacre domains that yields coalignment of the c’ axis of the crystallites, requires a 3D orientation field. Work is underway in this direction.

Finally, we note that the thickness of the organic and mineral layers in the nacre are typically 10–40 nm and 300–600 nm,[136−140] respectively. Our qualitative simulations give a considerably larger relative thickness for the organic layer. This is partly because we intended to model the whole granular → columnar → nacre sequence, and because of limitations of available computational power, we cannot have sufficient spatial resolution to realize a more realistic thickness ratio. If modeling is limited to the nacre, one is expected to achieve a better agreement.

Shell-Like Microstructures in Model PF2

To relax some of the simplifying assumptions made in PF1, a refined model (PF2) was proposed for modeling the formation of mollusk shells in ref (27). In this model, two solid phases form simultaneously from the liquid state, a mineral-rich and an organic-rich, while a fixed relative orientational relationship is forced between the solid phases formed inside the same crystal grain. This realizes a strong orientational coupling between the solid phases. Simultaneous formation of two solids occurs, for example, in eutectic or peritectic systems. In the refined approach, we opted for the former case. The main assumptions that set the conditions, under which eqs , 2, 4, 6, and 7 defining model FP2 were solved, are as follows:[27]

Crystal growth of CC happens via molecule/ion attachment.

A binary eutectic model thermodynamics (regular solution) applies.

The mineral content of the extrapallial fluid emitted at the surface of the mantle decreases exponentially with time.

Formation of granular CC crystals starts by heterogeneous nucleation on organic heterogeneities, whose number density is assumed to decrease exponentially with the distance from the periostracum.

The thickness of the extrapallial domain (distance between the mantle and the solidification front) remains constant. (In the simulation, the position of the mantle surface varies in accord with the solidification rate.)

In this approach, the assumption that CC supersaturation decreases toward the mantle is removed and is replaced by the more natural assumption that the CC supersaturation at the mantle decreases exponentially with time that leads to an analogous result, however, in a more natural way. The simulation presented in ref (27) shows a good qualitative agreement with the experimental microstructures observed in various types of mollusk shells. It is reassuring that model PF2 recovers the experimental microstructure though in a somewhat more ordered form, which could however be made more random by varying the noise added to the equations of motion. At the present state of affairs, it is difficult to decide which of models PF1 or PF2 should be considered superior to the other, yet the two-solid model is probably closer to the reality.

Herein, we explore whether the microstructure remains similar, when assuming ACC mediated crystallization in the framework of model PF2. The results obtained by the two mechanisms of crystallization are compared Figure . Apparently, as in the case of model PF1, the two micromechanisms for crystalline CC formation yield similar microstructures.

Note that in model PF2, the nacre is composed of two alternating solid phases. Under the conditions used in our simulation the organic component remains in an amorphous or nanocrystalline state. Remarkably, the predicted “nacre” structure recovers such details of the experimentally observed microstructure as the “mineral bridges” and the line defects across the organic layers termed “aligned holes”.[33,34,138,140] We will show, however, in the next section that despite this close similarity, the formation mechanism of the nacre can be considerably more complex than predicted here.

Discussion of Results from Models PF1 and PF2

First, we compare the solidification rates obtained from the simulations performed assuming DL ≈10–9 m2/s (aqueous solution) for the mother phase, with those using DACC,ion ≈ 10–15m2/s (ionic diffusion in ACC). We wish to stress that the velocities evaluated from the simulations can only be considered qualitative, as the thermodynamic driving force we used in this study might be well away from the true ones, which in turn may differ from the value obtained for the pure aqueous solution–CCC system, or the ACC–CCC system. Work is underway to perform more quantitative phase-field simulations to determine the growth and dissolution rates in pure systems, for which reasonably accurate experimental data are available.

The growth rate results are presented in Figure . As one expects the average growth rate is roughly proportional to the diffusion coefficient of the mother phase. The growth velocities predicted for crystallization from ACC by models PF1 and PF2 are about two and one orders of magnitude higher (v ≈ 10–8 and 10–9 m/s, respectively) than the experimental ones[57−60] (≈10–10 m/s, see gray zone in Figure ), whereas the rates predicted for the ion-by-ion addition are about v ≈ 10–1 and 10–2 m/s, which are about 8–9 orders of magnitude too high. On this ground, the mechanism based on the fast ion-by-ion addition can be excluded. One may perhaps argue that the production of calcium and carbonate ions in the surface layer of the mantle (outer epithelium) may be the rate-limiting process, which may be then taken so slow as to match the experimental growth rate. However, in that case, it is not the diffusion in the mother phase that controls the time scale of the process, and thus, the mechanisms models PF1 and PF2 rely on would not be present.

Figure 11

Growth rates predicted by models PF1 (red lines) and PF2 (blue lines) for the ion-by-ion mechanism from aqueous solution and for crystallization from ACC precursor. Note the ∼7 orders of magnitude difference in growth velocity predicted for the two mechanisms. For comparison, the range of experimental data is also shown (gray domain). Note the oscillating growth rate in the layered domain.

The velocity vs time relationships show characteristic differences for the two models: PF1predicts a steeply increasing velocity, followed by a plateau decreasing slowly with time for the domain of alternating solid–liquid layers. In the present simulations, the early stage behavior of the two models is different due to reasons different of the growth mode: while model PF1 starts with surface induced heterogeneous nucleation on the periostracum, PF2 relies on volumetric heterogeneous nucleation on organic impurities (note that both heterogeneous nucleation mechanisms can be adopted in both models). As a result, in model PF2, fast initial crystallization is observed during the formation of the granular layer via volumetric heterogeneous nucleation. This is followed by steady-state growth (roughly constant growth velocity) in both the prismatic and the nacreous domains, due to the lack of long-range diffusion during fast eutectic solidification. In contrast, in PF1 a continuously decelerating solidification is observed, which is combined with oscillating growth rate in the nacre. The hypothesized mechanism (volumetric heterogeneous nucleation in a thick layer at the periostracum) that creates protection for the mollusk fast in the early stage of shell formation is advantageous from the viewpoint of survival.

We note that because of computational limitations the maximum linear size of the computational domain we used was about 26.4 μm both in the PF1 and PF2 simulations. However, analogous structures can be produced on a larger size scale via reducing the rate by which the driving force of crystallization decreases with position/time and with an extended initial domain filled with heterogeneous nucleation centers.

It is appropriate to mention that while our models describe the formation of the granular and prismatic layers and the sheet nacre reasonably well within the framework of directional solidification, the predicted mechanism for the formation of the nacre via alternately precipitating mineral and organic layers may be oversimplified.

This is especially true for columnar nacre: experiments indicate that the formation of a quasi-periodic network of organic membranes precedes the formation of the CCC layers, which fill the space between the organic membranes later, as illustrated in Figure .[136−140] Although one could imagine that the organic membranes form periodically in space via an oscillating chemical reaction of the diffusion-reaction type, while the extrapallial fluid fills the space between the membranes, apparently, the real mechanism is more complex: first a multilayer outer membrane forms, of which the individual layers are exfoliated by nucleating CC crystals.[136,137] Furthermore, recent experimental works show a thin ACC layer on the surface of the CCC tablets of the nacre,[141,142] implying that (h)ACC particles may play an essential role in the formation of nacre. This, in turn, indicates that a particle-by-particle addition is more appropriate for solidification, as in the case of coral skeletons (see discussion in section ).

Figure 12

Electron microscopy images displaying the formation of columnar nacre in the case of (a) Haliotis rufescens (reproduced with permission from ref (136), under Creative Commons Attribution 4.0 License, Copyright 1982 Authors; the linear size is about 7.1 μm), and (b) Phorcus turbinatus (reproduced with permission from ref (137). Copyright 2015 Trans Tech). Here ss and s stand for “surface sheets” and “organic sheet”, whereas sm and ilm indicate “surface membrane” and “interlamellar membrane”, respectively.

Summarizing, while models PF1 and PF2 cannot capture all details of the formation mechanism of nacre, remarkably similar microstructures are generated. This raises the possibility that on the basis of the mechanisms these models realize (diffusion controlled solidification at high driving forces), one may design and prepare artificial mollusk shell structures that inherit the mechanical excellence emerging from the hierarchical sequence of the granular, prismatic, and nacre-like ultrastructures. Thus, the present work opens up the way toward a novel design strategy for creating biomimetic/bioinspired composite materials.

Finally, to address the evolution of columnar nacre within the phase-field theory, one can incorporate preexisting organic membranes into the simulations “by hand”, including the mineral bridges/aligned holes seen in the experiments.[33,34,138,140] We used PF1 to model the formation of CCC stacks shown in Figure . The thin organic walls were assumed to have amorphous structure (local orientation varied randomly pixelwise) and we applied a boundary condition that ensured a contact angle of 100° with the solid–liquid interface. Simulations of this kind yield 2D “pyramid-like” stacks of CC layers as shown in Figure . Such simulations can be used to explore the effect of such parameters as growth anisotropy, contact angle, hole size/position, and so on. For example, Figure b implies that the growth velocity of stack height increases with increasing hole width.

Figure 13

Snapshots of 2D phase-field simulations (model PF1), showing the formation of columnar nacre-like CCC tile stacks (gray/yellow - different colors denote different crystallographic orientations) via solidification between preexisting organic membranes (brown horizontal lines) with mineral bridges provided by (a) aligned holes of uniform width (20 Δx), and (b) aligned holes of three different widths (from left to right, 20, 5, and 10 Δx, respectively). The simulations were performed on 2000 × 500 and 2000 × 1000 grids (a) with anisotropic interfacial free energy and with (b) kinetic and interfacial free energy anisotropy. The horizontal size of the simulations is 26.25 μm.

Helical Structures Predicted by Model PF3

Spectacular screw dislocation-like helical structures have been observed in mollusk shells akin to patterns formed in oscillating chemical reactions (Figure ). These 3D structures cannot be addressed in models PF1 and PF2 as the scalar orientation field used in them is valid in only 2D (as in 3D minimum the three Euler angle are needed to define the crystallographic orientation). Therefore, we use model PF3 to explore the possibility of forming such objects within the framework of the phase-field theory. Here two solid phases (α and β) precipitate from a homogeneous ternary liquid. Owing to the lack of relevant information, we retain the materials parameters used in ref (102). Under appropriate conditions, shown by the green diamonds in Figure , a layerwise structure composed of alternating α (mineral) and β (organic) layers form (the third component is water, the crystal grows into hACC), an analogue of the “nacre” observed in model PF2. Note that the layer-by-layer growth mechanism, by which the layered structure forms is spinodal nucleation of one phase on the other, and thus can be viewed as an extreme case of island growth. In this regime, we see the formation of helical structures that emerge in pairs of opposite chirality. As the spiraling eutectic dendrites in ref (102), these defects originate from an instability associated with diffusion of the third component.

Figure 14

Domain in an idealized ternary phase diagram,[102] in which screw dislocation like structures form (green diamonds). For comparison, domains of spiraling eutectic dendrites (red dots), and unordered eutectic structures (blue crosses) are also shown.

Different views and sections of such structures are shown in Figure . This chiral structure closely resembles the screw dislocation-like defects reported in experiments.[14,39] The similarity could be enhanced by incorporating kinetic/interfacial energy anisotropies yielding faceted growth perimeters. It cannot, however, be excluded that the screw dislocation-like helical structures have here a different origin than in the reality. For example, our dislocation pairs do not recombine even in the long time limit, presumably because of the lack of mechanical stresses that are not incorporated into the phase-field models used in this study. Further investigations are yet needed to clarify whether phase-field models of island growth[143] or anisotropic pinning[144] could provide more realistic dynamics.

Figure 15

Screw dislocation pair formed in model PF3 at composition c = (0.43, 0.55, 0.02). (a) cross section through the axes of the pair. (b) Top view. (c) Front view of the central part of the section shown in panel (a). (d) Multiple screw dislocation pairs formed at c = (0.44, 0.54, 0.02). Red and yellow colors indicate solid solution phases rich in the mineral and organic components. These are generic dimensionless computations (performed on 256 × 256 × 272 grids moving with the solidification front), illustrating the ability of model PF3 to capture the formation of helical structures, but they are not optimized for modeling screw dislocations in sheet nacre.

Simultaneous formation of alternating organic and mineral layers in the vicinity of topological defects predicted by our simulations is a working hypothesis. Previous work suggests that the formation of spiraling organic membranes precedes mineral deposition.[39,145] Biological systems often use liquid crystals based on chitin as a template.[146] It is hypothesized[39,40,145] in this case that self-organization of the liquid crystal phase leads to the formation of a helical organic scaffold, which serves as a template for the mineralized structures of the same configuration. This mechanism of scaffold formation may be analogous to the formation of helical structures in Liesegang systems,[147,148] which raises the possibility of using the phase-field inventory developed for such systems.[149]

Note that the present simulations were not optimized for the case of nacre. Further work is underway to characterize these structures and the dynamics of their formation within the phase-field theory.

Modeling of Coral Skeletons in Model PF1

Next, using model PF1, we try to find a qualitative answer to the question of why the skeleton of some corals species contain small randomly oriented crystallites, “sprinkles” (see Figure e), which occur at the perimeter and along grain boundaries and even form bands, whereas other species do not display this behavior, an observation discussed in some detail in ref (28). In our previous work, we demonstrated that conditions of mineralization can influence the amount of sprinkles; however, we have addressed only tangentially the formation of sprinkle bands.

In this section, we address the latter phenomenon. We hypothesize that the coral polyp sits on the corallite (top of the skeleton) and emits a supersaturated extracellular calcifying fluid (CF), which is not in direct contact with the seawater but fills in the cavities of the porous skeleton. Recent work indicates that the coral skeletons are deposited biologically and actively via attachment of hACC nanoparticles (of diameter 50 to 400 nm), while ion-by-ion addition fills the interstitial space among them.[150] This leads to the formation of a thin (<1 μm) amorphous surface layer: diffusion of solid hACC nanoparticles and ion diffusion in CF contribute to the deposition. The hACC surface layer crystallizes into aragonite plumose spherulites after dehydration.

As the hACC/ACC surface layer remains thin, crystallization is apparently fast enough to keep pace with hACC deposition. Thus, diffusion of the hACC nanoparticles is expected to be the rate-limiting phenomenon. Thus, although particle diffusion does not influence the thermodynamic driving force of the ACC → CCC transition, it does influence the effective crystal growth rate. The effective rate is small on the sides of the fingers due to the low concentration of nanoparticles and allows a longer time the consolidation of the crystal layer yielding larger grain sizes. At the tips, the effective crystal growth rate is larger because of a higher concentration of the nanoparticles; thus, a shorter time is available for consolidation, and the crystal grains remain smaller. This is observed in the experimental images (Figure a–c).

Figure 16

Comparison of the cross sectional microstructure of top part (corallite) of the skeleton of coral species Balanophyllia europea (a–c) to a simulation performed using by model PF1. The 2D simulation was performed on a 1000 × 2000 grid. Orientation maps obtained by PIC mapping are shown in (a–c). Reproduced with permission from ref (28) under Creative Commons license CC-BY-NC-ND 4.0. Copyright 2020 Elsevier. At high magnification, an abundance of small crystallites, sprinkles, can be seen around the centerline of the finger shown in panel (c). For comparison, the computed orientation maps are displayed in (d) and (e), whereas the computed composition map is presented in panel (f). In panels (b–e) different colors correspond to different crystallographic orientation. The white rectangles in (a) and (b) denote areas magnified in (b) and (c), respectively, whereas the white rectangles in panels (e) and (f) indicate the area shown magnified in panel (d). Panels (d–f) show results of a generic dimensionless computations that incorporate several dimensionless combinations of the relevant physical properties; accordingly, the corresponding physical size scale depends on the choice of these parameters. With appropriate choice of these parameters, the length scale of the experiments can be recovered.

The diffusion coefficient for hACC particles of radius R = 150 nm can be estimated by the Stokes–Einstein relationship Dp = kBT/(6πRpη), which leads to Dp ≈ 2 × 10–12 m2/s, for T = 300 K and a viscosity of η = 1 mPa·s. Taking a medium growth rate[61−63] of v = 5 cm/yr ≈ 1.6 × 10–10 m/s for the coral skeleton, one obtains a diffusion length of lD ≈ 0.92 mm, whereas for fast-growing corals (37 cm/yr)[63]lD ≈ 0.12 mm, which do not rule out diffusion-controlled solidification.

Accordingly, model PF1 is applied here as follows: the phase-field and the coupled (particle-) concentration field control solidification on the time scale of particle diffusion and lead to a Mullins–Sekerka-type diffusional instability[151] of the interface. The magnitude of the orientational mobility, in turn, determines whether the forming solid becomes orientationally homogeneous (single crystal), partly ordered (polycrystalline), or fully disordered (amorphous).

Since in the case of coral skeletons we do not need to produce nacre-like alternating organic–inorganic layers via a transition from diffusionless to diffusive growth, in model PF1 we can use conditions, where dv/D ≪ 1; that is, the system is far from the diffusionless growth regime. The dimensionless model parameters and boundary conditions are the same as those we used in ref (28) for coral skeletons.

The polycrystalline structure emerging under these conditions due to the diffusional instability is shown in Figure d–f. Note the rough surface and the liquid channels, with small crystallites at the tips and along the spine of the branches, and larger crystallites at the lateral surfaces. This distribution of grain size is the result of the fact that the growth velocity at the tips is larger because of the larger supersaturation (the tip meets fresh, nondepleted liquid, see Figure f), whereas in the interarm channels, the fluid is depleted, so growth is slow, yielding larger crystallites. This phenomenon is the result of growth front nucleation (GFN) that leads to more frequent GFN events with increasing growth velocity as discussed in refs (75−78,152). Whether the combination of diffusional instability with GFN or some other biology directed mechanisms is responsible for the appearance of the rough surface of the skeleton is unclear at present. It is thus desirable to investigate further consequences of the hypothesized control of the grain size distribution via diffusional instabilities and GFN.

Along these lines, we make predictions using model PF1 that may be tested experimentally and validate or disprove our assumptions. For this reason, we investigate the effect of model parameters that influence the amount of sprinkles forming in the model and thus may offer explanation why it varies from one species to the other.

Previous work indicated that the orientational mobility (related to the rotational diffusion coefficient), the thermodynamic driving force (influenced by supersaturation or temperature) may influence the intensity of GFN (formation of new grains at the growth front).[75−78,153] In the framework of this study, we varied individually these parameters, and in all cases, we obtained a transition from microstructures dominated by sprinkles to microstructures with few or no sprinkles. Of these parameters, only the temperature can be controlled with relative ease. Variation of particle concentration in the calcifying fluid below the coral polyp, or controlling the rotational diffusion coefficient of the particles in the fluid are probably beyond the reach of the experimenter. Therefore, we present only the microstructural/morphological changes predicted as a function of temperature (see Figure ). Apparently, according to our model, the coral skeletons grown at low temperatures should have larger amount of sprinkles than those grown at higher temperatures. This finding raises the question whether the amount of sprinkles is indeed a characteristic feature of the individual coral species (i.e., determined biologically) or some differences in the circumstances of mineralization (temperature and/or supersaturation) are responsible for the deviations. In any event, these simulations may offer a natural explanation for the differences seen in the amount of sprinkles in the skeleton of different coral species.

Figure 17

Morphology and microstructure evolution as a function of temperature as predicted by model PF1. The reference computation is shown in panel (d). The relative temperatures from left to right are ΔT = −3, – 2, – 1, 0, 1, and 2 K. The simulations were performed on 1000 × 2000 grids. Note the decreasing amount of sprinkles with increasing temperature (from left to right). These are generic dimensionless computations that incorporate several dimensionless combinations of the relevant physical properties; accordingly, the corresponding physical size scale depends on the choice of these parameters. With an appropriate choice of them, the length scale of the experiments can be recovered.

Conclusions

We have demonstrated that qualitative coarse-grained phase-field modeling offers a methodology to address specific mesoscale aspects of nonclassical crystallization phenomena taking place during biomineralization. In particular, we investigated to what extent qualitative phase-field modeling can contribute to the understanding of microstructure evolution during the formation of mollusk shells and coral skeletons for various species. We applied three different phase-field models: PF1, PF2, and PF3.

Our present findings can be summarized as follows:

Ultrastructure specific to the shells of mollusks Unio pictorum, Nautilus pompilius, and Haliotis asinina: Driving the solidification process from solute trapping toward partitioning via decreasing the thermodynamic driving force, binary phase-field models PF1 and PF2 recover the common sequence of granular → prismatic → nacre ultrastructures on a reasonable time scale, if CC crystallization takes place via an amorphous precursor. In contrast, within this scenario, CC crystallization via ion-by-ion deposition from aqueous solution appears to be orders of magnitude too fast when compared to experiments.

Nacre formation in mollusk shells: Models PF1 and PF2 describe reasonably well the formation of not only the granular and prismatic domains, but the appearance of sheet nacre as well, in which case the models indicate alternating precipitation of the organic and mineral components. The models seem to reproduce even such details as mineral bridges and aligned holes. Yet, for obvious reasons, they cannot predict the formation mechanism of columnar nacre, in which the formation of organic membranes precedes CC precipitation. However, representing the preexisting organic membranes via appropriate boundary conditions, a reasonable description can be obtained even in this case.

Screw dislocations in mollusk shells: Ternary phase field model PF3 predicts the formation of screw dislocations pairs in 3D, a phenomenon analogous to the experimental findings. Inclusion of elasticity into the model is needed to capture the proper dynamic behavior during growth.

Sprinkle formation in coral skeletons: model PF1 was used to explore the possible mechanism for the formation of nanoscale crystallites “sprinkles”, whose presence was reported recently in the skeletons of certain coral species. Assuming a diffusion controlled mechanism in confined space, we observe the formation of sprinkle bands at the spine of the arms of the corallite as a trace of fast solidification at the arm tips, whereas larger crystallites form at the sides of the arms. The simulations show that varying the orientation mobility (proportional to the rotational diffusion coefficient of the molecules/ions) or the driving force of crystallization (via changing either the supersaturation or the temperature), one can control the amount of sprinkles between essentially no sprinkle and dominantly sprinkled microstructures.

Unquestionably, the applied models should be viewed as only minimum models of the processes taking place during biomineralization. Yet, we believe that phase-field modeling complemented with biochemical/biological information has the potential to contribute to a better qualitative or even quantitative understanding of morphogenesis in simple cases. Introduction of more complex models can certainly improve the mathematical representation of the associated phenomena. Ultimate limitations of such approaches stem from the fact that living organisms cannot be modeled within this framework: they can only be represented by boundary conditions of different complexity. Despite these, in specific cases, we recovered structures closely resembling their biogenic counterpart. The resemblance of the simulated and natural biominerals suggests that, underneath the immense biological complexity observed in living organisms, the underlying design principles for biological structures may be so simple that they can be understood with simple math and simulated by phase-field theory.

Finally, we note that our simulations outline conditions, under which standard materials science processes can be used to create inorganic substances that mimic the microstructures observed to form in living organisms, a knowledge that may open up ways for creating new biomimetic/bioinspired composite materials.