A computational framework for the morpho-elastic development of molluskan shells by surface and volume growth

Mollusk shells are an ideal model system for understanding the morpho-elastic basis of morphological evolution of invertebrates’ exoskeletons. During the formation of the shell, the mantle tissue secretes proteins and minerals that calcify to form a new incremental layer of the exoskeleton. Most of the existing literature on the morphology of mollusks is descriptive. The mathematical understanding of the underlying coupling between pre-existing shell morphology, de novo surface deposition and morpho-elastic volume growth is at a nascent stage, primarily limited to reduced geometric representations. Here, we propose a general, three-dimensional computational framework coupling pre-existing morphology, incremental surface growth by accretion, and morpho-elastic volume growth. We exercise this framework by applying it to explain the stepwise morphogenesis of seashells during growth: new material surfaces are laid down by accretive growth on the mantle whose form is determined by its morpho-elastic growth. Calcification of the newest surfaces extends the shell as well as creates a new scaffold that constrains the next growth step. We study the effects of surface and volumetric growth rates, and of previously deposited shell geometries on the resulting modes of mantle deformation, and therefore of the developing shell’s morphology. Connections are made to a range of complex shells ornamentations.


Introduction
With around 200,000 living species, molluska are the second most diversified phylum of the animal kingdom, including gastropods (snails, slugs), bivalves (mussels, oysters,...), cephalopods (squids, Nautilus,...) and five other classes [1] occupying a wide range of marine, freshwater, and terrestrial habitats.The huge morphological diversity among classes makes mollusks particularly interesting from an evolutionary perspective, notably with regard to questions related to the origin, evolution, and disparity of their body plan and their shell [2,3].The evolutionary success of mollusks, spanning over 540 million years, can be at least partly attributed to the shell that provides both protection and support to the soft body [4].Beyond their obvious aesthetic appeal, molluskan shells are an important research area in different fields.They have become exemplar model systems for studying the processes of biomineralization, a topic attracting a great deal of interest: from materials science to biomedical applications [5,6].Recent studies have begun to identify genes involved in these complex processes and to analyse how they are developmentally regulated (e.g.[7]), although the physical mechanisms underlying the morphogenesis of the shell ultrastructures remain poorly understood [8].Recent attention has also been given to the formation and differentiation of the shell-secreting mantle margin during development [9] and its morphological variations among classes [10,11].Detailed microscopic studies continue to provide important details about the structure and mutual relationships between the mantle, periostracum, and shell [12].However, despite their importance to many fields, the morphogenetic processes underlying the diversity of shapes remains elusive.This poor state of knowledge may lead to an incomplete, if not a distorted, view of the mechanisms underlying their morphological evolution.
Several interesting theories have addressed the formation of pigmentation patterns.However, these theoretical models invoking either reaction-diffusion chemical systems [13] or nervous activity in the mantle epithelial cells [14] cannot, by themselves, explain the emergence of three-dimensional forms that are subject to forces during the organism's development and life span.Indeed, while colour patterns on surfaces are primarily of biochemical origin, the formation of three-dimensional ornamentations such as ribs, tubercles, and spines is mostly a mechanical problem resulting from force generation on the mantle during growth, and its distortion in response to the force.Early studies have also considered the role of mechanics in the development of molluskan shells [15][16][17][18][19].More recently, some of the authors have developed a general framework of mollusk shell morphogenesis based on continuum theories of growth and mechanics [20,21].These models have been used to study the development and evolution of shell shape from a biophysical perspective [22][23][24][25][26].In particular, these morpho-mechanical models suggest that three-dimensional ornamentations, either parallel (i.e.commarginal ribs) or orthogonal (i.e.antimarginal ribs) to the growth lines do not require prefiguring patterns at the molecular level but may emerge de novo as a result of the dynamic balance of mechanical stresses intrinsic to the secreting system constrained in its growth by the calcified shell edge to which it adheres.
Following these simplified models, we present a full three dimensional numerical framework to study the accretive growth and nonlinear elastic deformations of the secreting mantle.As a first application, we study the effect of the geometry of the calcified shell edge, surface growth rate and morpho-elastic growth rate of the mantle on the resulting elastic deformation modes.We next study how these parameters may interact during shell development to generate diverse forms.Our main motivation for focusing on generic physical processes involved in development is that they may shape living beings in a predictive way and partly determine the spectrum of forms that have been and could have been generated during evolution.This outlook can be traced back 100 years to the pioneering work of D'Arcy Thompson, whose 1917 tome "On Growth and Form" [27] continues to inspire a growing community of researchers in various fields of theoretical, evolutionary and developmental biology (e.g.[21,28]).In this perspective, computational models of morphogenesis constitute an important tool to uncover the non-contingent rules that physical processes introduced in the development and evolution of forms.Molluskan shells grow via an accretive process occurring at the shell margin by a part of the anatomy called the mantle, which is a thin elastic membrane lining the inner surface of the shell.At each growth increment, the mantle extends slightly beyond the calcified shell edge, adheres to the rigid shell, and secretes matrix proteins, which, through biomineralization and calcification harden into a new layer of shell.

Mollusk shell growth mechanics
Within this process is an interesting mechanical interaction, due to the fact that the mantle is itself a part of the living mollusk, and thus undergoing growth separate from the growth of the shell.As the mantle may have grown since the last shell secretion, the mantle margin may be longer than the shell edge, and hence attachment to the shell may induce deformation of the mantle tissue that is then "cemented" in the shell shape upon secretion and calcification.From a mechanical point of view, shell growth may thus be summarized by the steps illustrated schematically in Fig. 1: (1) the mantle extends beyond the shell edge (surface growth of the mantle) while also growing along the shell margin (volumetric growth), followed by (2) attachment to the rigid shell, creating (3) an elastic deformation (morphoelasticity); in this deformed configuration, (4) new shell material is secreted, which can be viewed as surface growth of the shell, and thus (5) a new layer of shell appears in the shape of the deformed mantle, which undergoes biomineralization and calcification, and the process repeats.
The same basic process occurs in all shell-building mollusks, and yet produces a hugely diverse output of shell shapes and ornamentations.A general goal is to produce a mathematical and computational framework to explore this diversity: in particular how mechanical properties of the mantle, growth rates, and geometry conspire to produce the beautiful and varied outcomes observed in this phylum.However, a complete mathematical description is inherently challenging, as it links complex shell geometry (helicospiral, e.g.), elements of both surface and volume growth, nonlinear elastic deformations, and calcification.Previous work by some of the authors has approached this problem in a setting of one-dimensional elasticity, treating the interaction between the mantle margin and shell as a rod on an evolving elastic foundation.Here, our objective is to develop an algorithmic approach and computational tools to model the problem in a setting of three-dimensional nonlinear morphoelasticity.
However, the mathematical details of such a computation are quite complex; indeed the combination of surface and volume growth is itself a significant challenge in biomechanics, with evolving reference configurations and multiple growth tensors.Here, we have the added complexity of shell calcification.Hence, a proper description involves the delicate treatment of surfaces evolving due to combined mechanisms of growth, mechanical forces, and a calcification front that plays the role of a boundary condition.The mathematical details are provided in the following subsection (which may be skipped by the reader whose interest lies primarily in the outcome of applying the growth models).In short, our algorithmic approach is to execute the steps shown in Fig 1 .In step (1), surface growth happens by extension of the mantle, while its volumetric growth of the mantle is assumed to occur only in the direction parallel to the shell edge.The calcified portion of the shell then acts as a boundary condition constraining one edge of the mantle margin in step (2).Elastic deformation in step (3) is determined via mechanical equilibrium within the framework of finite elasticity (and computed via finite elements).Subsequent secretion, i.e. surface growth of the shell, then occurs over the extended and deformed mantle in step (4).Configurations of both the mantle and the shell are updated in step (5), and the process is repeated.As discussed further below, there are variations possible on exactly how this process is implemented in a computational setting.Our objective in this paper is not to exhaustively explore all possibilities, but rather to demonstrate the general validity of the algorithm and examine some basic properties of the shell patterns that emerge as output.

Mathematical framework
We detail the mathematical framework that is the foundation for the eventual computational treatment of the growth processes outlined above.This requires descriptions of surface and volume growth, elasticity and calcification.In our model, de novo shell material is configured by a combination of mechanisms among those introduced above: (a) surface growth by mantle extension along a unit vector, s 1 , which is tangential to the shell surface and perpendicular to the mantle margin, i.e. s 1 denotes the general direction of shell growth; (b) growth in size of the mantle (therefore, "volume" growth) manifesting in its expansion along another unit vector s 2 , which is tangential to the shell surface and in the direction of the mantle margin;1 and (c) formation of crests and troughs nominally perpendicular to the undeformed shell surface along the unit vector s 3 .To be precise, we have where the triad {s 1 , s 2 , s 3 } changes along the curved shell surface (Fig 2).The third mechanism arises from elastic bifurcations from a smooth surface, and post-bifurcation deformation driven by "excess" mantle growth.This is the morphoelastic mechanism, which is susceptible to a continuum mechanical treatment.It is key to development of the elaborate, antimarginal decorations of molluskan shells [20,21,[29][30][31][32].At the outset we remark that the need for precision in describing the array of configurations and mechanisms leads to complexity of notation.

Surface growth of the mantle
The mid-surface of the shell is represented by the surface, S τ ⊂ R 3 .We regard S τ as a one-parameter family of surfaces, generated by τ ∈ [0, T ], from a reference surface S 0 ⊂ R 3 .The generating curve, Γ τ ⊂ ∂S τ , is the leading edge of S τ and evolves perpendicular s 1 (see Fig 2).For points X(τ ) ∈ Γ t and X(0) ∈ Γ 0 , where Γ 0 ⊂ ∂S 0 is the initial generating curve at time τ = 0, we have X(τ ) = χ τ (X(0)).Surface growth occurs by extension of the mantle along the boundary curve Γ τ , which advances with the velocity χ = v 1 s 1 .In our computational studies, we Local coordinate system on the surface patch.s 1 is the direction of surface growth due to mantle extension, s 2 is the direction of volume growth along the mantle margin, and s 3 = s 1 × s 2 is the normal to the local surface patch.The calcified shell edge after time τ 1 forms the generating curve Γ τ1 for the τ th 1 step; the leading edge of the grown and deformed mantle strip then forms the generating curve, Γ τ2 for time step (τ 1 , τ 2 ).
will consider spatial and temporal variations in v 1 . 2We approximate the shell and the extended mantle as maintaining a constant thickness along s 3 throughout the growth process.

Volume growth of the mantle
We next consider volume growth of the mantle by expansion along s 2 , which is also the tangent to Γ τ . 3Due to this mechanism of growth the arc length of the fully relaxed mantle increases over time.Our treatment is focused on the kinematic manifestation of possibly inhomogeneous volume growth along s 2 .We adopt the framework of finite strain elasticity, with one important variation on the traditional theme: We regard surface growth as a process of deposition of material.The reference configuration of a material point is then determined by its deposition time.A family of reference surfaces (2-manifolds in R 3 ) is defined: 3 , parameterized by the time of deposition, τ .A material point lies on a reference surface, X(τ ) ∈ ω 0τ if it was deposited at time τ .The point-to-point map of material points X(τ ) from the reference surface ω 0τ to x(t; τ ) on the current surface, ω tτ , is x(t; τ ) = ϕ(X(τ ), t) = X(τ ) + u(X(τ ), t), where u is the displacement field.Note that ω tτ ∈ R 3 also is a 2-manifold.The primary strain quantity is the the deformation gradient, defined as F (X(τ ), t) = ∂ϕ(X(τ ), t)/∂X.Morpho-elastic growth of the soft mantle tissue is modeled by the multiplicative decomposition of the deformation gradient F (X(τ ), t) = F e (X(τ ), t)F g (X(τ ), t) into elastic and growth components, respectively [21,31].
The intuitive idea is that with the mid-surface of the shell being represented by S t , the mantle's "preferred" state is given by the growth tensor F g relative to S t × (−h/2, h/2).Because of its attachment to the rigid shell the mantle cannot attain F g , but only F , with F e being the elastic incompatibility.This multiplicative decomposition of the kinematics is the framework of morphoelasticity.It depends on the time of deposition of material points, and therefore on evolving reference configurations, and is depicted in Fig 3 .In the above parametrizations t ≥ τ is understood.In what follows, we will suppress functional and parameteric dependencies wherever there is no danger of confusion.
As expressed above, our key kinematic assumption on volume growth of the mantle is that it occurs only along s 2 , so that, accounting for the appropriate tangent spaces between which F g acts, ( Here, ε 2 is the rate of the growth strain along s 2 .As with the surface growth velocity, we will consider spatial and temporal variations in ε 2 .

Secretion of shell material
The scaffold for de novo deposition of shell material is the mantle that has undergone an increment of surface and volume growth.New shell material is secreted on the mantle's dorsal surface.

The calcification front
While the mantle can be treated as an elastic solid, the calcified shell itself is rigid.An advancing calcification front, C τ ∈ ϕ(S τ ), is the interface between the rigid shell and material recently secreted by the mantle.At time τ the calcified shell is in the configuration denoted by Ω s t . 4The velocity of C τ is v c , which lies in the plane defined by {s 1 , s 2 }.

Algorithmic formulation and implementation
The first step toward an algorithmic implementation is a discretization of the continuous processes of surface growth, morphoelastic volume growth, and evolution of the calcification front.The time interval of interest t ∈ [0, T ] is discretized by instants t 0 , t 1 , . . ., t N , into sub-intervals [t 0 , t 1 ], . . ., [t N −1 , t N ], where t 0 = 0 and t N = T .For simplicity, we also consider deposition times τ = t 0 , t 1 , . . . .In a time step ∆t = t k+1 − t k , the leading surface of material secreted by the mantle, ω tt , advances by δss 1 = v 1 ∆ts 1 .Fig 4 depicts the relevant geometry (generating curve, mantle front surface in reference and deformed configurations parameterized by deposition time) and the growth processes driving the mantle's shape by surface growth and morphoelastic volumetric growth.Secretion of shell material and calcification are implied, but not shown The preceding mathematical model is continuous in time.It describes the biological processes in the sequence of (1) the mantle's surface growth (extension), (2) morphoelastic volumetric growth, (3) shell growth by secretion on the mantle's current configuration, followed by (4) calcification.However, the discrete model operates with time steps ∆t = t k+1 −t k .While the time-continuous setting led to complex notation for evolving configurations, the time-discrete setting allows some simplifications in this regard.The four stages are implemented in a time-parallel setting, meaning that no order is implied over [t k , t k+1 ].We note that the time discretization reflects a time-discontinuous, e.g. a diurnal process of growth and calcification.

Surface growth
At time t k , the shell has been fully calcified: C t k = ϕ(Γ t k ).In the time interval [t k , t k+1 ], the leading surface ω tt k (in its deformed configuration) is displaced due to extension of the mantle by v 1 ∆ts 1 to a reference surface ω 0t k+1 .This allows us to define a strip of the mantle in its reference configuration Ω m 0t k bounded by the surface ω tt k at its trailing edge and ω 0t k+1 at its leading edge.See Fig. 4 2

.2.2 Volume growth
Volume growth of the mantle is obtained by integrating Eqn (1).We exploit the exponential map: where The actual deformation gradient achieved is F k+1 , with the elastic component F e k+1 = F k+1 F g −1 k+1 being governed by nonlinear elasticity.With a strain energy density function ψ(F e ) that satisfies frame invariance (so that ψ(F e ) = ψ(F eT F e )), the first Piola-Kirchhoff stress tensor is It is governed by the quasistatic stress equilibrium equation imposed at time t k+1 : In our computations we apply Dirichlet boundary conditions u = 0 on the trailing surface ω 0t k of the mantle where it meets the rigid shell.A combination of fixed Dirichlet, u = 0, and traction-free Neumann conditions, P N k+1 = 0 are applied on the remaining boundaries ∂Ω m 0t k+1 \ω tt k .This defines the morphoelastic growth problem for mapping the mantle strip from its reference configuration Ω m 0t k to its deformed configuration Ω m tt k .

Secretion
Following morphoelastic volumetric growth of the mantle in the algorithmic step described above, a virtual step occurs, in which new material is secreted over the deformed configuration of the mantle Ω m tt k .

Calcification
The final step of the algorithm is propagation of the calcification front so that C t k+1 = ϕ(Γ t k+1 ).The mantle strip is calcified into its deformed configuration, Ω m tt k .Remark 1: The above algorithm is a manifestation of our observation that it is over the mantle that both surface growth and morpho-elastic volumetric growth occur.The actual formation of new shell material by secretion over the mantle and calcification follow once the current mantle configuration has been defined by surface and volumetric growth., also.Equation ( 4) is then to be solved over of . Variants on this idea also are admissible, including complete calcification of a subset of Ω m tt k by time t k+1 , so that C t k+1 does not coincide with ϕ(Γ 0t k+1 ).

Implementation
The formulation has been implemented in a finite element framework using standard hexahedral elements.Two to five layers of elements are used through the thickness.The attachment of the mantle to ω tt k and its extension up to ω 0t k+1 is implemented by extending the finite element mesh by one or more rows of elements as shown in Fig 5 .This is followed by the growth law in Equation ( 2) subject to the constitutive law Equation ( 3) and the governing equations (4).Following Remarks 1 and 2, secretion of shell material is a virtual step over the deformed mantle configuration Ω m

Results
In the framework constructed above, there are three key parameters determining the ornamentation pattern that develops through shell growth: 1.The active mantle width, i.e. the amount of surface growth occurring in each time increment, given by δs = |v 1 s 1 − v c |∆t.
2. The volume growth over each time increment: δg, which is related to the growth rate 3. The initial curvature, κ, of Γ 0 .
These three governing parameters are illustrated in Fig. 6.Our objective is to explore the morphological space of patterns that results from variations in these parameters and explain them on a mechanistic basis while making connections to ornamentations observed on molluskan shells.
In the first instance, we study the morphologies that result from varying each parameter in isolation.In each case, we use our computational framework to impose either (1) a single, finite increment of surface growth manifested in a specified active mantle width, or (2) an increment of morphoelastic volume growth, or (3) observe the influence of distribution of curvature along s 2 on the reference curve, Γ 0 .The effects will be characterized by the wave number (the number of crests and troughs) along the lip of active mantle in the s 2 direction, the amplitudes and the locations of crests and troughs.In a growing molluskan shell, these effects are coupled, and potentially dynamically changing through development.We therefore proceed next to analyze the coupled effects of variations in the above three parameters, as well as of spatially and temporally varying surface and volume growth.The reference generating curve, Γ 0 , with its curvature, κ along s 2 ; the active mantle width, δs; volumetric growth increment, δg.

Variation in surface growth via the active mantle width
The effect of varying surface growth by the active mantle width is studied for fixed volume growth rate and reference curvature.With our computational formulation, we solve for the resulting shape after a single surface growth increment of active mantle.We start from a reference shell edge, Γ 0 , that is an arc of a circle with curvature κ = 0.01.Fig. 7 shows the resulting morphologies when varying the active mantle width δs = |v 1 s 1 − v c |∆t by up to eightfold.We see that an increase in δs leads to a decrease in wave number, which can be understood as follows: Increasing δs places the free edge of the active mantle strip, ω 0t k+1 , further from the rigid boundary ω tt k (k = 1, 2, . . . ) in the s 1 -direction, thus decreasing its structural stiffness to bending.The excess material (along s 2 ) created by the volume growth increment δg can therefore be accommodated by an increased deformation in the s 3 direction, without paying a large strain energy penalty.As a result, each increment in δs induces fewer wave crests/troughs, each with a larger wavelength and amplitude.Fig. 7 also presents a comparison with the ornamentations on Clinocardium nuttallii and Timbellus phyllopterus.Here, it is instructive to compare to an analogous reduced order model: a one-dimensional rod on an elastic foundation.In this model, an elastic rod that has an excess of length due to axial growth is connected elastically to a rigid support: a curve representing the calcified shell edge.The support provides an external force that resists displacement of the rod away from the foundation (i.e.displacement in the s 3 direction in our framework).This system has been formulated in detail by some of the authors [34] and forms the basis of previous mechanical descriptions of shell morphogenesis [22,24].In the linearized system, with the rod and foundation extending along the x-axis, the deformed rod has shape y(x) satisfying [34] y (x) + (γ − 1)y (x) + kγy(x) = 0, Here γ > 1 describes the axial growth, analogous to δg in the computational framework.The parameter k is proportional to the stiffness of the elastic foundation, and therefore models the stiffness of the active mantle strip to deflections of the mantle margin.Considering for simplicity an infinite rod and seeking a solution of the form y ∼ exp(2πinx), the preferred bifurcation mode corresponds to the smallest value of γ * > 1 for which (5) has a solution; this is found to be γ * = 1 + 2k + 2(k + k 2 ) 1/2 , from which we obtain that the wave number at buckling satisfies n = √ γ * − 1/(2 √ 2π).From this we can extract the scaling relationship n ∼ √ k for large k.Based on the intuitive argument above, we would posit an inverse relationship between k and δs, e.g.k ∼ δs α with α < 0, i.e. an increase in strip width acts to decrease the effective foundation stiffness, leading to a decrease in bifurcation mode.To further explore this relationship, we extract the dependence of mode number n on δs from the simulations presented in Fig. 7, and plot the comparison in Fig. 8.Because of the highly nonlinear, post-bifurcation states of deformation, n was defined as the number of waves, ignoring the dependence of amplitude and wavelength on the coordinate in the s 2 -direction.From the slope of this log-log plot, we get n ∼ δs −1 , and so k ∼ δs −2 .

Growth rate
In our model, the morphoelastic volume growth of the soft mantle is the origin of the pattern of bifurcation: This mechanism generates an excess of length in the active mantle strip relative to the rigid shell edge, whose distribution is given by the variation of δg along s 2 .A compressive stress is thus induced in the mantle.As is well-understood within the theory of finite strain elasticity, when this in-surface stress exceeds a critical threshold, a bifurcation occurs and the mantle relaxes into a lower energy configuration by deflecting out of the mid-surface in the local s 3 direction.Fig. 9 illustrates the effect of varying the growth increment δg by up to a factor of 4× while holding the incremental active mantle width fixed at δs = 1.0.The reference curve Γ 0 has curvature of κ = 0.01 in the middle of its extent in the s 2 direction, decreasing toward zero near the ends.Note that this induces a length scale: the radius of curvature r = 100, relative to which δs = 1.
In a time-continuous process, there would exist a critical time, t cr , and corresponding amount of morphoelastic volumetric growth, g cr = tcr 0 ε 2 (t)dt for which the compressive stress crosses the critical threshold and the corresponding bifurcation mode appears.Modes with greater numbers of waves have increased energy.Then, as volume growth continues beyond g cr the compressive mantle stress and amplitude of the first observed bifurcation mode both increase, until the stress exceeds a second critical threshold corresponding to another bifurcation and a higher mode appears.This effect is demonstrated in Fig 9, where initially the compressive stress (corresponding to δg = δg * ) initiates bifurcation into a mode with three discernible crests, n = 3.A growth increment to δg = 2δg * gives a compressive mantle stress that exceeds a higher critical threshold, and is accommodated by an increase in the mode number to n = 7. Increasing the growth to δg = 3δg * , 4δg * only increases the amplitude of the seventh mode, with no further bifurcations.Remark 3: A careful examination of the bifurcated shape with wave number n = 3 for δg = δg * reveals a deformation with amplitude that is highest at the midpoint of the arc of Γ 0 , where the curvature κ = 0.01, decreases in the two immediate neighbor crests, and decays by more than an order of magnitude toward the lateral edges, where κ → ∞.Inclusion of all crests regardless of amplitude would raise the inferred wave number to n = 7 even for this first volume growth increment.We understand this as a cascade in which the lowest mode to appear (n = 3) is localized to the higher curvature region.At δg = δg * , there is a super-position, with the n = 5 and n = 7 modes also present, but at lower amplitude.The coincidence of crests for modes n = 3, 5, 7 in the high curvature region leads to the higher amplitudes there.At the very next increment to δg = 2δg * the strain energy settles into the seventh mode, whose prominence is magnified.In contrast, for κ = 0.01, but uniform as in Fig. 7, a single mode is observed, whose amplitude is uniform provided it does not merge into the lateral boundary.Also consider Fig. 10a with κ → ∞ and uniform, where the amplitude remains uniform.The localization of mode shapes is a consequence of curvature, which we examine in greater detail in Section 3.1.3.

Curvature
Taking a cue from the localization of modes in high κ regions of the reference curve, Γ 0 , we next consider this effect in greater detail.We consider three geometries for Γ 0 : a line with curvature κ → ∞, a curve with κ having almost uniform sign, and a second curve with κ changing sign along the arc.The result appears in Fig 10 .For all three cases in the figure, the Dirichlet boundary conditions are u = 0 on the trailing surface, ω 0t 0 , and the lateral surfaces of the mantle, which are perpendicular to s 2 .
We find that the deformed configuration of the mantle Ω m tt 0 is biased toward developing curvature of the same sign as the reference curvature, κ, along s 2 .With Γ 0 a straight line in Fig 10a, the mantle deforms upward and downward equally, i.e. curvatures of both signs are see in Ω m tt 0 , and the crests/troughs have the same amplitudes.The (mostly) single-signed κ of Fig 10b promotes like-signed crests and suppresses oppositely signed ones.This is also apparent in Fig 10c, which has regions where κ takes on positive and negative signs.Thus, we see the influence of geometry in inducing compliance to mantle deformation by forming crests that are compatible, and resistance to forming crests that are incompatible with the reference curvature, respectively.This pattern of deformation is consistent with the greater amplitude of the central, compatible crest in Fig. 9 with δg = δg * .This result has an intriguing relevance for mollusk shell ornamentation: the reference shape of the shell on which ornamentation appears, modelled here by Γ 0 , is generally convex and surrounds the mollusk living inside.It would be disadvantageous for the ornamentation to appear inward, as this would intrude on the mollusk's living space.A natural question then is whether the mollusk must execute a developmental plan to ensure that the ornamentation is built in the outward direction.The results here suggest, however, that the geometry and growth mechanisms naturally conspire to bias the pattern outward.The influence of the geometry of reference curves on antimarginal ornamentation.For fixed active mantle width, the amplitudes of crests in the deformed configurations are magnified if they are compatible, and attenuated if they are incompatible, respectively, with the reference curvature.These reference curves bear similarity to the shape of the mantle edge found in (a) Clinocardium, (b) Pterynotus phyllopterus and (c) Bolinus brandaris.The underlying spatial discretization (mesh) is also shown on the model geometries.

More complex patterning
The influences of the surface and volume growth rates, and of the geometry via reference curvature, have been established.We now consider the combination of these effects, their temporal and spatial variations, in two mechanisms that lead to more complex ornamentations.

Progression of ornamentation with volume growth
As demonstrated above, geometry exerts its influence by magnifying the amplitudes of crests in mantle deformation that are compatible, and attenuating those crests that are incompatible, respectively, with reference curvature.It is of interest to study the progression of these crests with continued volume growth.With this aim, we consider a shell edge with the geometry of

Hiearchical ornamentation by temporal variation of growth rates
As a second approach to complex patterning we study the effect of multiple generations of surface and volume growth.
Since the examples in Section 3.1 have already considered a range of either surface or volume growth rates varying individually, we now consider the effect of combining these growth rates while also allowing them to vary in time.Our aim in this section is to identify a mechanism that explains the secondary and even tertiary crests and troughs that are visible along the shell edge in species such as Hexaplex cichoreum.We recognize these as the remnants of increasing wave number during shell growth.Noting that the decrease in wave number with increasing active mantle width in Section 3.1.1also implies an increase with decreasing active mantle width, and recalling the magnification of higher modes with increasing volume growth rates in Section 3.1.2,we consider the following protocol of surface The resulting shell morphology thus has a hierarchical structure to its ornamentation, with higher modes appearing over configurations that initially have lower modes.Such features are present in the ornamentations of a number of muricid species including Hexaplex cichoreum and Hexapelx duplex.Indeed, as shown in the H. cichoreum shell in Fig 13, a tertiary buckling mode also appears, with some shells even showing quaternary modes in a fractal-like structure.Such morphological features are within reach of our model in principle, although resolving details beyond secondary modes becomes challenging due to the complexity associated with ensuring curvature continuity in the s 1 direction with accumulation of high curvature crests.The combined modulation of surface growth rate (active mantle width), volume growth rate and curvature presents a simple mechanical basis for the morphogenesis of hierarchical ornamentation in seashells, which has not previously been described.

Ornamentation with negative Gaussiann curvature due to spatial variation in growth rate
An examination of the mantle deformation in Fig. 7 and Figs 9-13 shows that the majority of crests and troughs form with positive Gaussian curvature.One exception is Fig. 10b.This case is explained by the strain energy due to high curvature, κ of the reference curve, Γ 0 , being relieved by development of negative Gaussian curvature in the active mantle strip.The other interesting case is in Fig. 13, where the positive Gaussian curvature after the first two growth increments, i.e., over (t 0 , t 1 ] and (t 1 , t 2 ] up to mantle configuration Ω m tt 1 , changes into negative Gaussian curvature after the final volumetric growth increment in Ω m tt 2 .Taking a cue from the temporal variation in volume growth over (t 2 , t 3 ], we recognize that it also imposes a spatial variation: The mantle strip of width δs * formed by surface growth over (t 0 , t 1 ] experiences volume growth δg = δg * , but the strip of length δs = 1 2 δs * from surface growth over (t 1 , t 2 ] experiences a total volume growth of δg = 3δg * .We are therefore prompted to consider that the growth rate ε 2 is an increasing function along the s 1 direction, i.e., ∂ε 2 /∂ξ 1 > 0, where ξ 1 is the curvilinear coordinate defining s 1 .In this instance, within a single growth increment there is greater excess of length at the leading edge of the active mantle strip compared to the trailing edge.The elastic mantle attains a locally energy minimizing configuration by adopting negative Gaussian curvature of the deformed mantle. An example of directly imposing such spatial variation is shown in Fig 14, where the profile of ε 2 (ξ 1 ) has low but positive curvature ∂ 2 ε 2 /∂ξ 2 1 > 0 for small ξ 1 , changing smoothly to high curvature ∂ 2 ε 2 /∂ξ 2 1 0 for larger ξ 1 .The variation in growth rate generates a deformation with a finite component in the negative s 1 direction, i.e. the mantle "arches back" to accommodate the excess length, creating ornamentation with negative Gaussian curvature.
While, as suggested by our computations and demonstrated in Fig. 7 and Figs 9-13, antimarginal ornamentation in shells is often with Gaussian curvature that is positive or appears to vanish, there are a number of species that display such negative Gaussian curvature.In Fig 14(c) we show a top view of Hexaplex chicoreum, which displays a strongly backward arching ornamentation, similar to the mantle deformation in Fig. 14.Here we have another instance of a mechanical basis for a feature of ornamentation in seashells for which no mechanistic explanation has previously been advanced, and that can be reproduced in our computational framework for coupled surface and volume growth.

Discussion and concluding remarks
Mechanics has been recognized as a framework for explaining biological growth and form since at least the appearance of D'Arcy Thompson's work [27].However, a large part of the literature on morphological aspects of growth since the 1970s, such as that assembled by Meinhardt [41] and others, has focused on applying analytic or semi-analytic generating curves to the forms of shells, horns and antlers.The coupling of three-dimensional form to material forces and displacements, one aspect of which is morphoelasticity, has remained a more difficult problem.The difficulty stems from the complexity attained by the coupled equations, especially where nonlinear elasticity appears, and has been very well laid out in [20] and [21].Consequently, it is only with the marriage of mathematics and numerical methods that general, three-dimensional, initial and boundary value problems have been solved [42,43].The literature on computational treatments of biological growth also has, in our eyes, suffered a limitation: Problems addressable by the model of inhomogeneous, volume growth, i.e., morphoelasticity, have formed the mainstay of this body of work.Effective as this treatment has been in explaining tumor growth [44,45], aspects of cardiovascular systems, and the folding of soft, layered structures during morphogenesis [46], it cannot be elegantly extended to accretive, surface growth.For such problems, the morphoelastic treatment is restricted to representing advancing fronts by a thickening surface layer.Under its effect, the front's motion is an emergent phenomenon that is controlled by local, pointwise, volumetric growth.Neither the elaborate, generated surfaces, nor their accompanying elastic fields can be represented by such an application of volume growth with its foundations in local volume changes, rather than de novo deposition of material.
Against this backdrop, we have presented, to the best of our knowledge, the first combined computational framework for accretive, surface growth and local, morphoelastic, volume growth.The mathematical basis for this framework in terms of generating surfaces, evolving reference configurations and moving fronts has been crucial because it has provided a rigorous foundation on which to elucidate the discretized, space-time formulation, as well as the finite element framework.The discretized space-time is a faithful reflection of the coupled processes of accretive surface growth and morphoelastic volume growth.There are alternatives to the finite element framework, however.Variants of level-set, phase field and immersed boundary methods would allow propagation of surface growth and the calcification front by fractions of an element width.Arbitrary changes in the propagation directions s 1 and v c could also be easily represented.We do not, however, see that this restriction to propagation by integral element widths presents a fundamental limitation in the shell morphologies and ornamentations that can be represented by the approach presented here.The advantages

Fig 2 .
Fig 2. Local coordinate system on the surface patch.s 1 is the direction of surface growth due to mantle extension, s 2 is the direction of volume growth along the mantle margin, and s 3 = s 1 × s 2 is the normal to the local surface patch.The calcified shell edge after time τ 1 forms the generating curve Γ τ1 for the τ th 1 step; the leading edge of the grown and deformed mantle strip then forms the generating curve, Γ τ2 for time step (τ 1 , τ 2 ).

Fig 3 .
Fig 3.The kinematics of growth.The observed deformation gradient, F , is composed of an incompatible growth component, F g , and another incompatible, but elastic component, F e , which restores compatibility of F .

Fig 4 .
Fig 4. Surface and volume growth.Extension of the mantle by v 1 ∆t in each time increment, ∆t, and its morphoelastic volumetric growth.Shell secretion on the mantle and calcification are implied in this figure.

Remark 2 :
The steps presented above impose full calcification of the secreted material in deformed configuration Ω m tt k .Consequently, the mantle in its reference configuration Ω m 0t k+1 attaches to the rigid material surface ω tt k+1 .An alternate model with possible biological relevance is to assume that Ω m tt k has not been fully calcified, leaving an elastic deformed configuration Ω m tt k .Then the attachment of the mantle reference configuration Ω m 0t k+1 at time t k+1 and its morphoelastic volume growth over [t k+1 , t k+2 ] further deforms Ω m tt k tt k .Calcification is imposed by turning Ω m tt k rigid for time t ≥ t k+1 .A number of examples are considered of mollusk shells displaying the shapes and marginal ornamentation that best demonstrate the interplay between surface growth, morphoelastic volumetric growth and calcification.In each case, the boundary condition u = 0 on the trailing surface ω tt k of the mantle where it meets the rigid shell.The "lateral" boundaries, ∂Ω m 0t k \ω 0t k \ω 0t k+1 , are subject to Dirichlet conditions on either u, or its normal component, with traction-free Neumann conditions, P N = 0 on the remaining boundaries ω 0t k+1 .See Fig. 5 for an illustration of the evolving mantle and calcified shell configurations.

Fig 5 .
Fig 5. Space-time discretization in the finite element framework: Evolution of the seashell surface through surface growth and morphoelastic volumetric growth of mantle strips followed by their calcification.

Fig 7 .
Fig 7. Effect of incremental active mantle widths on mantle deformation for fixed δg and κ: Increasing the active mantle width δs = |v 1 s 1 − v c |∆t over Γ 0 leads to morphologies bearing a similarity with the ornamentations on Clinocardium nuttallii.The Dirichlet boundary conditions applied are u = 0 on the trailing surface (ω tt k ).Inset image of Clinocardium nuttallii modified from source [33].Original image licensed under Creative Commons Attribution-NonCommercial 4.0 International License.

Fig 8 .
Fig 8. Scaling study of wave number n versus the incremental active mantle width δs.

Fig 9 .
Fig 9. Effect of volume growth increment δg = ε 2 ∆t on the amplitude and wave numbers of the deformed mantle.The observed high mode number morphologies are similar to the ornamentations observed in Bivalvia class of mollusks like Clinocardium nuttallii.The Dirichlet boundary conditions applied are u = 0 on the trailing surface (ω 0t 0 ).Inset image of Clinocardium nuttallii modified from source [33].Original image licensed under Creative Commons Attribution-NonCommercial 4.0 International License.

Fig 10 .
Fig 10.The influence of the geometry of reference curves on antimarginal ornamentation.For fixed active mantle width, the amplitudes of crests in the deformed configurations are magnified if they are compatible, and attenuated if they are incompatible, respectively, with the reference curvature.These reference curves bear similarity to the shape of the mantle edge found in (a) Clinocardium, (b) Pterynotus phyllopterus and (c) Bolinus brandaris.The underlying spatial discretization (mesh) is also shown on the model geometries.
Fig 10(c), having varying curvature, and impose volumetric growth ranging from δg = 0.0 to δg = 0.56.Several of these mantle deformations are shown in Fig 11.The trend observed in Fig 10-of magnification and attenuation of mantle deformation that is respectively compatible and incompatible with the reference curvature-continues.Favored crests display progressive magnification of amplitudes with continued growth.Also shown are Bolinus brandaris shells with progressively increasing amplitudes of crests corresponding to the mantle deformations in our computations. 5Fig 12 examines the smoothness of geometry.Mild reference curvature singularities leave virtually no visible trace on mantle deformation following volume growth.However, strong reference curvature singularities promote compatible crests, and remain visible as mild singularities in the deformed mantle lip.Smoothly varying reference curvature replicates the trend of favoring compatible crests.

Fig 11 .
Fig 11.Progression of curvature-compatible ornamentation with volume growth δg.The deformed mantles show marginally preferred localization around points of high curvature and thereafter the amplitude increases with growth.Some of these shapes with different amplitudes can be observed in the shells of the species Bolinus brandaris.Inset images of Bolinus brandaris modified from source[35][36][37].Original images licensed under Creative Commons Attribution-Share Alike License.

Fig 12 .
Fig 12.The influence of mild versus strong reference curvature singularities, and of smooth curvatures.

Fig 13 .
Fig 13.Hierarchical ornamentation arising from temporally increasing volume growth rates ε 2 .In the Hexaplex cichoreum image shown in the inset, three levels of ornamentation hierarchy are shown: primary (red) as a low mode bifurcation from a flat surface, secondary ornamentation as a second mode bifurcation (magenta) and tertiary ornamentation mode as a third mode (blue).The corresponding first, second and third modes are traced on the mantle edge of the numerical model.Inset images of Hexaplex cichoreum and Hexaplex princeps modified from source[38,39].Original images licensed under Creative Commons Attribution-Share Alike License.

Fig 14 .
Fig 14.Spatially varying volume growth We impose volume growth that varies along s 1 direction, with an increasing gradient toward the leading edge, as shown in (a).The result appears in (b), displaying large, negative Gaussian curvature, mimicking the strongly backward arching morphology observed in a number of shell species, for example as seen also in Hexaplex chicoreum.Inset image of Hexaplex chicoreum modified from source [40].Original image licensed under Creative Commons Attribution-Share Alike License.