Theoretical characterization of a model of aragonite crystal orientation in red abalone nacre 9 0 0 2 S.N. Coppersmith, P.U.P.A. Gilbert, and R.A. Metzler n Department of Physics, University of Wisconsin, Madison, WI 53706 a J 8 Abstract. Nacre, commonly known as mother-of-pearl, is a remarkable biomineral 2 that in red abalone consists of layers of 400-nm thick aragonite crystalline tablets confined by organicmatrix sheets, with the (001) crystalaxes of the aragonitetablets ] t oriented to within ±12o from the normal to the layer planes. Recent experiments f o demonstrate that this orientational order develops over a distance of tens of layers s from the prismatic boundary at which nacre formation begins. . t a Our previous simulations of a model in which the order develops because of m differential tablet growth rates (oriented tablets growing faster than misoriented - ones) yield patterns of tablets that agree qualitatively and quantitatively with the d n experimentalmeasurements. Thispaperpresentsananalyticaltreatmentofthismodel, o focusing on how the dynamical development and eventual degree of order depend on c model parameters. Dynamical equations for the probability distributions governing [ tablet orientations are introduced whose form can be determined from symmetry 1 considerations and for which substantial analytic progress can be made. Numerical v simulations are performed to relate the parameters used in the analytic theory to 2 8 those in the microscopic growth model. The analytic theory demonstrates that the 4 dynamical mechanism is able to achieve a much higher degree of order than naive 4 estimates would indicate. . 1 0 9 0 : 1. Introduction v i X Nacre, or mother-of-pearl, is a biomineral that attracts the attention of materials r a scientists, biologists, and mineralogists as well as physicists because of its remarkable mechanical properties and its incompletely elucidated formation mechanisms [1–5]. Aragonite, a hard but brittle orthorhombic CaCO polymorph, accounts for 95% of 3 nacre’smass, yetnacreis3000timestougherthanaragonite[6]. Nosyntheticcomposites outperform their components by such large factors. It is therefore of great interest to understand the mechanisms governing nacre formation. Nacre is a layered composite in which organic matrix (OM) sheets alternate with aragonite layers, each of which consists of tablets of irregular polygonal shape that completely fill the space between preformed OM layers. In red abalone, the OM layers are 30-nm thick, and the aragonite tablets are of thickness 400-500nm and width 5-6 microns [7,8]. The aragonite tablets are crystalline and oriented with their (001) crystal axeswithin±12◦ fromthenormaltothelayerplane[8,9]. TEMandAFMmeasurements Model of crystal orientation in nacre 2 have shown that the OM has pores through which aragonite can grow [10,11], and x-ray and x-ray photoelectron emission spectromicroscopy (X-PEEM) measurements demonstrate that nacre has stacks of co-oriented tablets [7,8,12], consistent with the hypothesis that aragonite crystals grow through pores in the OM sheets. Growth through pores explains how crystal orientation of the aragonite tablets is maintained between layers in the material, but it does not necessarily explain the physical mechanism giving rise to the orientational alignment in the first place. It is commonbeliefinthebiomineralizationfieldthatthealignmentofthearagonitecrystalc- axes is due to microscopic chemical templation by the OM[1,13], with organicmolecules providing surfaces that promote aragonite nucleation with preferred orientations. However, it is not obvious why such a mechanism would lead to a high degree of orientational order of the c-axes but not of the a or b axes of the aragonite tablets. Moreover, Ref. [8] reports x-ray photoelectron emission spectromicrosopy (X-PEEM) and microbeam X-ray diffraction that probe the degree of orientational alignment of the aragonitec-axes[7,14]anddemonstratethattheorientationalalignmentofthearagonite tablets increases systematically over a length scale of tens of microns starting from the prismatic boundary at which nacre growth originates. As discussed in Ref. [8], this observation suggests that the aragonite crystal orientation is the result of a dynamical process. In Refs. [8,15] a dynamical model is proposed for the development of orientational order of the aragonite tablets in nacre in which the presence of co-oriented stacks of tablets plays an essential role, and the ordering arises because oriented tablets grow faster than misoriented ones. Refs. [8,15] present numerical simulations of the model that successfully reproduce several different aspects of the pattern of tablet orientations in red abalone nacre measured using X-PEEM. This paper presents a more detailed investigation of the model and its behavior, including closed-form analytic results that apply in the limit that the probability of nucleating a tablet with orientation different from the one immediately below is small, which, based on comparison to the X-PEEM results, is the physically relevant regime. We find that the dynamical ordering mechanism can lead to a remarkably high degree of orientation of the tablet c-axes, because the width of the distribution of tablet orientations decreases very strongly as the fraction of tablets that nucleate with the “wrong” orientation decreases. The analytic characterization of the model uses methods similar to those used to study mutation-selection models in population biology [16–18]. However, the analytic formulation involves some parameters whose relationship to those of the growth model is not determined in Ref. [8]. These parameters are examined here, and it is shown that this relationship is not trivial. The relationship between the parameter sets depends on the spatial arrangement of the nucleation sites, and some aspects are quite insensitive to changes in the values of the microscopic parameters. Numerical simulations and mean- field arguments are used to relate the parameters in the analytic probabilistic model to the parameters in the original growth model for the case when the nucleation sites are chosen randomly with uniform probability on each layer. Model of crystal orientation in nacre 3 The paper is organized as follows. Section 2 presents the model, while section 3 presents the analysis of the probabilistic model that enables one to understand qualitatively some features of the behavior, using methods that have been developed to study models relevant to population biology [16–22]. The fixed point behavior is discussed in subsection 3.1, while the dynamical evolution is discussed in subsection 3.2. Section 4 discusses the relationships between the parameters of the growth model to those of the model for the probability distribution for tablet orientations. Sec. 5 is a discussion, and the conclusions are presented in Sec. 6. Appendix A presents additional details of the arguments justifying the functional forms of the equations used in the main text. 2. The model Figure 1 illustrates the basic mechanism of nacre growth [1,10,23]. First, organic matrix sheets spaced by approximately 0.4 microns are created. Crystalline aragonite tablets thennucleateonthefirstlayer atuncorrelatedrandomlocationsandgrowwhileconfined by the organic matrix sheets. The crystals in this layer continue to grow and fill out the space in the layer [1,24,25]. Each nacre tablet in a given layer nucleates and then grows until it reaches confluence with a neighboring tablets, so that the resulting tablet pattern resembles a Voronoi construction [2,5], with tablet in-plane width of order 5 microns. The crystal orientation of each tablet is highly probable to be the same as that of the tablet directly below its nucleation site [10,11,26], which reflects the presence of pores in the organic matrix (typically with diameter ∼ 5 − 50 nm) through which aragonite crystals can grow. As discussed in [2,8], there is one nucleation site per tablet, with an identifiable structure in the organic matrix that is large enough to have one or more pores going through it. The model that we examine here assumes that growth in a given layer is completed before tablets in the succeeding layer are nucleated, and that the positions of nucleation sites inallsuccessive layers arenucleated at locationsthat areuncorrelated withthose in the preceding layers. The growth rates in the first layer are chosen uniformly at random in the interval [1−δ/2,1+δ/2], and the tablets in each layer grow to confluence. With probability 1 − ǫ a tablet has the same growth rate as the tablet below its nucleation site, while with probability ǫ the tablet is assigned a growth rate chosen uniformly at random from the range [1 − δ/2,1 + δ/2]. Appendix A presents the arguments that these choices of the probability distributions are natural, given the geometry of the experimental system. The version of the model used here and in Ref. [15] is simpler than that examined in Ref. [8], which explicitly models columnar nacre by assuming that growth in a given layer is not completed before tablets in the succeeding layer are nucleated, so that nucleation sites in successive layers must be within a certain in-plane distance and are thus correlated; it is theoretically more attractive because it requires the specification of one fewer parameter. Ref. [15] presents numerical simulations of the model and experimental data on the pattern of tablet orientations in red abalone Model of crystal orientation in nacre 4 Figure 1. Schematic of growth of sheet nacre, which exhibits layer-by-layer growth. (a)Scaffoldingoforganicmatrixsheetsiscreated. (b)Crystallinearagonitetabletsare nucleated at random positions in the first layer of pre-existing scaffolding of organic matrix sheets. (c-e) Aragonite tablets are confined by organic matrix sheet in the next layer, but grow within the layer until they reach confluence. (f) Aragonite tablets are nucleated at random positions in the next layer. With high probability, the nucleated tablet has the same crystal orientation as that of the tablet directly below the nucleation site [10,11,26]. nacre,anddemonstratesthatnumericalsimulationsofthismodelusingparametervalues ǫ = 0.015 and δ = 0.25 yields good agreement with the experimental measurements on nacre from red abalone, Haliotis rufescens. The next section presents an analytic treatment of the behavior of the probability distribution function governing the density of tablets of different orientations as a function of the distance from the prismatic boundary. 2.1. Growth model used for simulations This paper reports simulation of a model in which nucleation sites are placed randomly on a square calculational domain with open boundary conditions. In each layer, the rate of growth of each tablet is chosen at random from a uniform distribution in the range [1 − δ/2,1 + δ/2]. The tablets nucleate simultaneously and grow to confluence before successive layers are nucleated. The resulting pattern of tablet boundaries in each layer is a multiplicatively weighted Voronoi construction [27–29]. Variations in the in-plane extent of each tablet within each individual layer are ignored, so that the volume of a given tablet is the product of its area and the layer separation. Because faster-growing tablets tend to take up a greater fraction of the area, and because a tablet will nucleate with the same velocity as the tablet below it with probability 1 − ǫ, when ǫ is small, there is a tendency for faster-growing tablets to take up more and more of the total area (and hence volume) as the growth proceeds. The ordering is not perfect, however, because some misoriented tablets are nucleated. Model of crystal orientation in nacre 5 3. Effective theory for the evolution of probability distribution of tablet orientations The theory in this section is formulated in terms of crystal orientations, as opposed to the formulation in terms of growth rates in the previous section. It is useful to considerthisformulationbecausetheX-PEEMexperimentsmeasurecrystalorientations as opposed to growth rates, so comparison of theoretical predictions to experimental results is facilitated. Motivated by the ordering of tablet c-axis orientations observed experimentally, our model posits that variability in the orientations of the c-axes of the tablets gives rise to variability in the tablet growth rates, with tablets with c-axes aligned perpendicular to the layer plane growing the fastest. This section presents and analyzes dynamical equations that govern the evolution of the probability distribution describing the number of tablets with different crystal orientations. These equations can be solved analytically in the limit that the fraction of misoriented tablets is very small, which is the parameter regime relevant to the experiments of Ref. [8]. We find that the geometry of the dynamical model leads to the result that the distribution of tablet orientations (or, equivalently, growth rates) has an extremely sharp peak whose width is exponentially small in ǫ, the fraction of misoriented tablets. The analytic theory presented here for the evolution of tablet orientations in nacre is closely related to mutation-selection models studied in the context of population biology [17–22]. The faster growth of tablets of a particular orientation is analogous to theprocessofselectioninpopulationbiology,wherespecieswithhigherfitnessreproduce faster than species with lower fitness. In the growth model examined in this paper, after all the tablets in a given layer grow to confluence, the next layer nucleates. When a tablet in the next layer nucleates, one of two things can happen. The first, more likely, possibility is that the nucleated tablet has the same orientation as the tablet directly below its nucleation site (analogous to inheritance, where the descendant has thesamefitnessastheparent), whilethesecondisthatthetablethasarandomlychosen orientation (analogous to mutation, where the fitness of the descendant differs from that of the ancestor by a random amount). Selection in the growth model is reflected in the tendency for a larger fraction of the area to be filled with tablets with higher growth rates, which increases the probability that a given randomly located nucleation point is located over a tablet with a higher growth rate. Theanalytictheorypresented hereissimilar toGinzburg-Landautheories[30]often usedinstatisticalphysics [31]todescribephasetransitionsincondensed mattersystems, in that the functional forms follow from symmetry and dimensionality considerations and may involve unknown coefficients. Even if the numerical values of these coefficients are not known, the analytic theory is very useful for obtaining insight into the interplay between the various parameters in the problem. However, additional insight is obtained if the Ginzburg-Landau parameters can be related to the parameters of the growth model described above in Sec. 2.1, which is done in Sec. 4 using a simple mean-field Model of crystal orientation in nacre 6 theory and also by comparing the predictions of the analytical model to the results of numerical simulations. We define φ (γ)dγ to be the fraction of tablets in layer ℓ that are misoriented from ℓ the layer normal by angles that are between γ and γ + dγ. We assume that tablets with γ = 0 have the largest rate of growth, so their share of the area in layer ℓ+1 will tend be greater than in layer ℓ, thus leading to a larger fraction of tablets with crystal axes oriented parallel to the layer normal; the function w(γ) governs this tendency. We define χ (γ)dγ to be the fraction of the area in layer ℓ after its growth is completed that ℓ is oriented in the range of angles between γ and γ+dγ from the layer normal, and that 1 χ (γ) = w(γ)φ (γ) , (1) ℓ ℓ N ℓ where, for each ℓ, the normalization factor N is determined by the normalization of the ℓ probability: π 1 = dγ χ (γ) (2) ℓ Z0 π ⇒ N = dγ w(γ)φ (γ) . (3) ℓ ℓ Z0 The maximum growth rate is at γ = 0, and one expects w(γ) to have a quadratic maximum at γ = 0. We will scale w(γ) so that w(0) = 1. We assume that w(γ) is a functiononlyofγ (inotherwords, thatitdependsonlyonthedegreeofmisorientationof the c-axis and not on the orientation of the a and b axes). One expects the dependence of w(γ) on γ near its maximum at γ = 0 to be quadratic, which, if the distribution is not too broad, can be approximated as a Gaussian, αγ2 w(γ) ∝ 1− ≈ e−αγ2/2 , (4) 2 (cid:18) (cid:19) where α is a numerical coefficient. The tendency for fast-growing tablets to take up an increasing fraction of the total tablet area is analogous to the effects of selection in population biology, where organisms with higher fitness tend to comprise an increasing fraction of the population. We then assume that most of the tablets that nucleate in layers above the first layer have the same crystal orientations as those of the tablets just below their nucleation sites, but that there is a small probability ǫ that a tablet nucleates with a value of γ that is chosen at random from a normalized probability distribution f(γ). The nucleation of misoriented tablets in the nacre growth model is analogous to the effects of mutation in a population genetics model. The combination of these growth and nucleation terms leads to a dynamical equation governing the behavior of φ (γ): ℓ φ (γ) = ǫf(γ)+(1−ǫ)χ (γ) ℓ ℓ 1 = ǫf(γ)+ (1−ǫ)φ (γ)w(γ) . (5) ℓ−1 N ℓ To complete the definition of the model, one must specify an appropriate form for f(γ),thefunctiondescribing thedistributionofanglesofmisorientedtabletorientations. Model of crystal orientation in nacre 7 We will show below that the behavior depends only on the properties of f(γ) as γ → 0. In Appendix A it is found that the generic behavior for f(γ) for small γ is for f(γ) to be proportionaltoγ asγ → 0. Thisresult isintuitively reasonableforthissystem geometry because the angular area between γ and γ +dγ is proportional to γ as γ → 0. At some points in the analytic treatment below, we will choose the specific, mathematically convenient form f(γ) = βγe−βγ2/2. This choice does not affect any of the results, because the behavior of f(γ) as γ → 0 determines the asymptotic behavior. We characterize the behavior of this model using the methods of Refs. [16–18]. First we note that when ǫ = 0, so that the co-orientation of tablets in successive layers is perfect, this model is easily solved for any initial distribution, φ (γ). Since increasing 1 ℓ by one multiplies the un-normalized φ (γ) by w(γ), it follows immediately that ℓ φ (γ) ∝ φ (γ)(w(γ))(ℓ−1) . (6) ℓ 1 (Note that it is sufficient to compute the un-normalized distribution since the normalization factor for any given ℓ can always be obtained via Eq. (2).) The long- time behavior for any w(γ) with a single quadratic maximum at γ = 0 depends only on the curvature in w(γ) near γ = 0. For the specific choice w(γ) = exp(−αγ2/2), one finds φ (γ) ∝ φ (γ)exp(−α(ℓ−1)γ2/2). The width of the distribution decreases as the ℓ 1 square root of the number of layers and becomes arbitrarily narrow as the number of layers tends to infinity. Now consider the effects of nonzero but small ǫ, so that a nonzero fraction of tablets nucleate that are misoriented. The intuitive picture of the process in this regime is that the distribution of angles gets narrower unless nucleation of misoriented tablets occurs. The misorientations prevent the peak from narrowing indefinitely, so after many layers φ (γ) approaches a stationary distribution that does not change as ℓ increases further. ℓ Since a fraction ǫ of tablets is misoriented at each layer, one might expect that the peak in the distribution narrows for ∼ (1/ǫ) layers, so that this naive argument leads to the expectation that the eventual width of the distribution should be proportional to ǫ1/2. However, it is shown below that the subtle interplay between the effects of mutation and selection results in a peak width that can be exponentially small in ǫ as ǫ → 0. The high degree of tablet orientation obtained using this mechanism could be a significant advantage in the biological context. 3.1. Steady state behavior of the model Firstwefindthefixedpointbehaviorforthismodel,inwhichtheprobabilitydistribution function φ (γ) approaches a limit φ∗(γ) that is independent of ℓ. We expect that this ℓ fixed point distribution is reached in the limit ℓ → ∞. ∗ The equation determining the fixed point probability distribution φ (γ) is 1 φ∗(γ) = ǫf(γ)+ (1−ǫ)φ∗(γ)w(γ) , (7) N∗ with π ∗ ∗ N = dγ φ (γ)w(γ) . (8) Z0 Model of crystal orientation in nacre 8 The solution to this equation is ǫf(γ) ∗ φ (γ) = , (9) 1− 1−ǫw(γ) N∗ with π ∞ 1−ǫ m ∗ N = dγǫf(γ)w(γ) w(γ) . (10) N∗ Z0 m=0(cid:20) (cid:21) X We now define 1−ǫ v = , (11) N∗ π I = dγ(w(γ))m+1f(γ) , (12) m Z0 and rewrite Eq. (10) as ∞ 1−ǫ = ǫ vmI . (13) m v m=0 X A useful explicit form for the area distributions can be obtained in the parameter regime in which ǫ, the fraction of misaligned tablets, is small. When ǫ ≪ 1, the right handsideofEq.(13)canbeoforderunityonlyifvmI decaysslowlyforlargearguments. m Because of its definition, v > 0, and it will be seen below that normalization of the probability implies that v ≤ 1. Therefore, when m ≫ 1, because w(γ) has a single maximum at γ = 0, the integrand in I from Eq. (12) is very sharply peaked near m γ = 0. Therefore, the integration interval can be extended to [0,∞], one can assume that w(γ) is a Gaussian, w(γ) = exp(−αγ2/2), and only the behavior of f(γ) for small γ is relevant. If the orientations of the misoriented tablets are chosen uniformly at random in three-dimensional space, then, as discussed in Appendix A, f(γ) is proportional to γ as γ → 0, and for small argument f(γ) = βγ (with β a constant of order unity), so that ∞ I = β dγ γ exp(−α(m+1)γ2/2) = β/(α(m+1)). Therefore, when ǫ ≪ 1, m 0 ∞ R 1−ǫ β vm = ǫ . (14) v α m+1 m=0 X (Eq. 14 makes it particularly clear that a solution is possible only if v ≤ 1.) Using the identity [32] ∞ xk = −ln(1−x) , (15) k k=1 X one finds 1−ǫα = −ln(1−v) , (16) ǫ β so that 1−ǫ α(1−ǫ) N∗ = ≈ (1−ǫ) 1+ǫ − . (17) v βǫ (cid:18) (cid:20) (cid:21)(cid:19) Model of crystal orientation in nacre 9 Therefore, φ∗(γ) can be written ǫf(γ) ∗ φ (γ) = . (18) 1− 1−exp −α(1−ǫ) exp[−αγ2/2] βǫ The width of the probabili(cid:16)ty distribhution φ∗(iγ) in γ, estim(cid:17)ated by finding the value of ∗ ∗ γ at which φ (γ) = φ (0)/2, is 2 γ = e−(1−ǫ)α/(2βǫ) , (19) 1/2 α r which is extremely small when ǫ is small. 3.2. Dynamics of the analytic model The dynamics of the model defined in Eq. (5) can of course be obtained numerically. Analytic insight can also be obtained, following Ref. [18], by writing Eq. (7) for the fixed point function as ∞ 1−ǫ m ∗ φ (γ) = ǫf(γ) w(γ) , (20) N∗ m=0(cid:18) (cid:19) X and comparing this expression to the solution to Eq. (5), which can be written ℓ−1 m−1 1−ǫ ℓ 1−ǫ φ (γ) = ǫf(γ) w(γ) + w(γ) φ (γ) . (21) ℓ 1 N N "m=1 n=1 (cid:18) n (cid:19) n=1(cid:18) n (cid:19) # X Y Y If one assumes that the term proportional to φ (γ) is negligible and that the 1 normalizations N can be approximated as being the same, N = N, independent n n of n, then one obtains an expression for φ (γ): ℓ ℓ−1 1−ǫ m φ (γ) ≈ ǫf(γ) w(γ) . (22) ℓ N m=1(cid:18) (cid:19) X This expression is the sum of contributions that can be interpreted as describing the contribution of the population that has undergone a given number of selection events since the last mutation, and Ref. [18] shows that it agrees well with numerical solutions of Eq. (5). 4. Relating the parameters in the effective theory of Sec. 3 to the parameters in the growth model of Subsec. 2.1. Sec. 3 presents an analysis of a theory in which one writes an equation for φ (γ), the ℓ fraction of the area in layer ℓ in which the tablet orientation angle is γ. The model analyzed there contains three parameters. The first, α, defined in Eq. 4, governs the degree to which tablets with small values of γ grow faster than tablets with larger values of γ. The second, ǫ, is the fraction of nucleation sites that have tablets with a random orientation instead of the same orientation as the tablet below. The third parameter, β, specifies the width of the distribution governing the distribution of angles of the Model of crystal orientation in nacre 10 misoriented tablets. In contrast, the original growth model described in Subsec. 2.1 has two parameters, ǫ, the fraction of misoriented tablets, and δ, which governs the range of in-plane growth speeds of the misoriented tablets. This section describes the relationships between the two sets of parameters. The value of ǫ is the same in the analytic theory for orientations as in the growth model (hence the use of the same symbol). The parameters β and δ are related in a straightforward fashion, as discussed below in subsection 4.2. Most of this section will focus on the third parameter in the effective theory, α, which we show depends on the geometry of the nucleation sites in the growth model in subsection 4.3. We will present a simple mean-field theory for estimating the value of α for uncorrelated and random nucleation site locations which yields the correct order of magnitude for the value in subsection 4.3.1. We compare the mean field predictions with the results of simulations of the growth model in subsection 4.3.2. 4.1. Analytic theory for model formulated in terms of tablet growth rates In the growth model that is simulated numerically, the tablet orientation angles are not considered explicitly because the calculation is formulated using tablet growth rates, and the simulation is performed by choosing an initial configuration of nucleation sites with a distribution of growth rates, and then allowing the tablets to grow from these nucleation sites until they reach confluence. We choose to define a new variable x that ranges between 0 and 1, and define P (x), the probability distribution describing the ℓ relative frequencies of tablets in layer ℓ with different values of x [16]. It is natural to interpret x as v/v , where v is the in-plane tablet growth speed and v is the tablet max max growth velocity at the orientation where this velocity is maximum. As discussed above and in Appendix A, when the model is formulated in terms of orientation angles, the “fitness” function specifying the changes in the fractions of the area covered by tablets with different angles of misorientation between successive layers is expected to have a quadratic maximum at γ = 0. Because the tablet growth velocity itself depends quadratically on γ near γ = 0 (again, because misorientations by γ and −γ are equivalent and so yield the same growth velocity), the “fitness” function that specifies the change in relative area of thedifferent values of x, w˜(x), depends linearly on x near x = 1. Recalling that the probability distributions in the model are normalized, so that the overall scale of w˜(x) is arbitrary, and that the behavior is dominated by the behavior near x = 1, two ways of parameterizing this dependence are to (1) fix the value of w˜(1) = 1 and specify the slope w˜′, or (2) to write w˜(x) = xξ and specify ξ. ′ These two forms are equivalent near x = 1, with w˜ (x) = ξ; we will choose to use the power-law form w˜(x) = xξ. To obtain additional insight into the relationship between the two formulations of themodel,wereformulatetheanalytictheoryintermsofthevariablex(thisformulation is very similar to that in Ref. [16]). Assuming that tablets with normalized growth rate x grow to have an area that is proportional to xξ for some ξ, and assuming that there