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

TL;DR

This paper introduces a comprehensive 3D computational model that simulates mollusk shell development by coupling surface accretion and volume growth, providing insights into the morpho-elastic mechanisms shaping shell morphology.

Contribution

It presents a novel, general computational framework that integrates pre-existing shell shape, surface accretion, and volume growth to explain shell morphogenesis.

Findings

The model captures stepwise shell growth and deformation patterns.

Growth rates and initial geometries significantly influence shell morphology.

Connections to complex shell ornamentations are demonstrated.

Abstract

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.

Equations

$$\dot{\mbox{\boldmath$F$}}^{\text{g}}=\varepsilon_{2}\mbox{\boldmath$F$}^{\text{g}}\mbox{\boldmath$s$}_{2}\otimes\mbox{\boldmath$s$}_{2}.$$

$$\mbox{\boldmath$F$}^{\text{g}}_{t_{k+1}}=\mbox{\boldmath$F$}^{\text{g}}_{k}\exp[\varepsilon_{2_{k}}\Delta t\mbox{\boldmath$s$}_{2}\otimes\mbox{\boldmath$s$}_{2}].$$

$$\mbox{\boldmath$P$}=\partial\psi/\partial\mbox{\boldmath$F$}^{\text{e}}.$$

$$\text{Div}\,\mbox{\boldmath$P$}_{t_{k+1}}=0\;\text{in}\;\Omega^{\text{m}}_{0_{t_{k}}}.$$

$$y^{\prime \prime \prime \prime }(x)+(\gamma -1)y^{\prime \prime }(x)+k\gamma y(x)=0,$$

$$D{\nabla }^{4}w+N{w}_{xx}=0,$$

$$w=0,w_{xx}+\nu w_{yy}=0\text{on }x=0,A,$$

$$w=0,w_y=0\text{on }y=0,$$

$$w_{yy}+\nu w_{xx}=0,w_{yyy}+(2-\nu )w_{xxy}=0\text{on }y=B,$$

$$w(x,y)=\sin \left(\frac{m\pi x}{A}\right)f(y)$$

$$A^4{\lambda }^4-\frac{A^2{m}^2}{{\pi }^2}(\tilde{N}-2{\lambda }^2)+m^4{\pi }^4=0,$$

$${\lambda }_1^{\pm }=\pm i\frac{\sqrt{Am\sqrt{\tilde{N}}\pi +m^2{\pi }^2}}{A},$$

$${\lambda }_2^{\pm }=\pm \frac{\sqrt{Am\sqrt{\tilde{N}}\pi -m^2{\pi }^2}}{A},$$

$$f(y)=c_1\exp (i{\lambda }_1^{+}y)+c_2\exp (i{\lambda }_1^{-}y)+c_3\cos ({\lambda }_2^{+}y)+c_4\sin ({\lambda }_2^{+}y).$$

$$f(0)=f^{\prime }(0)=0,f^{\prime \prime }(B)-\nu \frac{m^2{\pi }^2}{A^2}f(B)=f^{\prime \prime \prime }(B)-(2+\nu )\frac{m^2{\pi }^2}{A^2}{f}^{\prime }(B)=0.$$

Full text

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

Shiva Rudraraju 1, Derek E. Moulton 2, Régis Chirat 3, Alain Goriely2, Krishna Garikipati 4*

1 Department of Mechanical Engineering, University of Wisconsin-Madison, Wisconsin, USA.

2 Mathematical Institute, University of Oxford, Oxford, UK.

3 Université Lyon1, CNRS UMR 5276 LGL-TPE, France.

4 Departments of Mechanical Engineering and Mathematics, University of Michigan, Ann Arbor, Michigan, USA. *Corresponding author. Email: krishna@umich.edu.

Abstract

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 [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 remain 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 from the 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 fully 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. [28, 21]). 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.

Mollusk shell growth mechanics

Molluskan shells grow via an accretive process occurring at the shell margin by an organ called the mantle, which is a thin elastic membrane lining the inner surface of the shell. Over each increment, the mantle extends slightly beyond the calcified shell edge, while adhering to the rigid shell. The mantle then secretes matrix proteins, which, through biomineralization and calcification harden into a new layer of shell.

Within this process is an interesting mechanical interaction, due to the fact that the form taken by the mantle along the growing shell front is fixed in the calcified edge, while the form of the calcified edge partly determines the shape of the mantle at the next growth increment [25]. As the mantle may have grown since the last shell secretion, its margin may be longer than the shell edge, and hence attachment to the shell may induce deformation of the mantle tissue that is then recorded and fixed in the shell shape upon secretion and calcification. We introduce our foundational notions of two distinct modes of growth. A process that creates new surface where none existed before is labelled as surface growth. In contrast, if growth takes place by deposition where material already exists, i.e., over pre-existing volume, followed by elastic relaxation, the effect is to locally increase the material volume without adding new surface. For this reason we label it as volume growth. 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 while also growing along the shell margin (volume growth), and (2) adhering to the rigid shell, creating an elastic deformation (morphoelasticity); in this deformed configuration, (3) new shell material is secreted (surface growth of the shell), and thus (4) 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 the distinction between the growing shell and the growing mantle, as well as the additional process 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 moving boundary condition. To simplify the description, in this paper we formulate a mathematical description of the process that treats the mantle and shell edge as a single elastic object undergoing surface growth, volume growth, and calcification. 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 process shown in Fig 1. In step (1) of the figure, extension of the mantle edge beyond the shell edge is modelled as surface growth. Volume growth is assumed to occur only in the direction parallel to the shell edge, and produces an excess of length of mantle relative to the calcified shell edge. Upon attachment, the calcified shell edge acts as a boundary condition constraining one edge of the mantle margin in step (2). Elastic deformation is determined via mechanical equilibrium within the framework of finite elasticity (and computed via the finite element method). This is the morphoelastic volume growth step that causes out-of-plane deformation of the mantle. Subsequent secretion of new shell material then occurs over the extended and deformed mantle in step (3) and the calcification front advances. Configurations of both the mantle and the shell are updated in step (4), 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.

Methods

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 (creation of new surface) by mantle extension along a unit vector, \mbox{\boldmaths}_{1}, which is tangential to the shell surface and perpendicular to the mantle margin, i.e. \mbox{\boldmaths}_{1} denotes the general direction of shell growth; (b) growth in size of the mantle over its pre-existing extent (therefore “volume” growth) manifesting in its expansion along another unit vector \mbox{\boldmaths}_{2}, which is tangential to the shell surface and in the direction of the mantle margin and shell edge (since the mantle takes on the shape of the shell surface near its leading edge, the mantle margin tracks the shell edge); and (c) formation of crests and valleys nominally perpendicular to the undeformed shell surface, and along the unit vector \mbox{\boldmaths}_{3}. We have \mbox{\boldmaths}_{3}=\mbox{\boldmaths}_{1}\times\mbox{\boldmaths}_{2}, and more specifically, \mbox{\boldmaths}_{i}\cdot\mbox{\boldmaths}_{j}=\delta_{ij}, where the triad \{\mbox{\boldmaths}_{1},\mbox{\boldmaths}_{2},\mbox{\boldmaths}_{3}\} changes along the curved shell surface (Fig 2). The third mechanism above arises from elastic bifurcations from a smooth surface, and post-bifurcation deformation driven by “excess” mantle growth relative to the previous shell increment. As explained in the Mollusk shell growth mechanics subsection of the Introduction, 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 [29, 20, 30, 31, 32, 21]. Although surface growth due to mantle extension perpendicular to the shell edge occurs only following volume growth along the shell edge, we have described the steps (a-c) in the order of surface growth, volume growth and morphoelastic deformation. This is for mathematical purposes only. In our model, these mechanisms are not separated in time. 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, $\mathscr{S}_{\tau}\subset\mathbb{R}^{3}$. We regard $\mathscr{S}_{\tau}$ as a one-parameter family of surfaces, generated by $\tau\in[0,T]$, from a reference surface $\mathscr{S}_{0}\subset\mathbb{R}^{3}$. The generating curve, $\Gamma_{\tau}\subset\partial\mathscr{S}_{\tau}$, is the leading edge of $\mathscr{S}_{\tau}$ and evolves along \mbox{\boldmaths}_{1} (see Fig 2). For points \mbox{\boldmathX}(\tau)\in\Gamma_{t} and \mbox{\boldmathX}(0)\in\Gamma_{0}, where $\Gamma_{0}\subset\partial\mathscr{S}_{0}$ is the initial generating curve at time $\tau=0$, we have \mbox{\boldmathX}(\tau)=\mbox{\boldmath\chi}_{\tau}(\mbox{\boldmathX}(0)). Surface growth occurs by extension of the mantle along the boundary curve $\Gamma_{\tau}$, which advances with the velocity \dot{\mbox{\boldmath\chi}}=v_{1}\mbox{\boldmaths}_{1}. In our computational studies, we will consider spatial and temporal variations in $v_{1}$. In principle, $v_{1}$ depends on space and time through quantities such as the density, stress, and chemical fields, among others, but we neglect such details in this preliminary communication. We approximate the shell and the extended mantle as maintaining a constant thickness along \mbox{\boldmaths}_{3} throughout the growth process.

Volume growth of the mantle

We next consider volume growth of the mantle by expansion along \mbox{\boldmaths}_{2}, which is also the tangent to $\Gamma_{\tau}$ (both surface and volume growth of the mantle as described are a consequence of growth in size of the mollusk. It is for purposes of mathematical modelling that we have distinguished the process into surface and volume growth.). Due 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 \mbox{\boldmaths}_{2}. We adopt the framework of finite strain elasticity, with one important variation on the traditional theme: The reference configuration of a material point is determined by its deposition time. A family of reference surfaces (2-manifolds in $\mathbb{R}^{3}$) is defined: $\omega_{0_{\tau}}=\Gamma_{\tau}\times(-h/2,h/2)\subset\mathbb{R}^{3}$, parameterized by the time of deposition, $\tau$. A material point lies on a reference surface, \mbox{\boldmathX}(\tau)\in\omega_{0_{\tau}} if it was deposited at time $\tau$. The point-to-point map of material points \mbox{\boldmathX}(\tau) from the reference surface $\omega_{0_{\tau}}$ to \mbox{\boldmathx}(t;\tau) on the current surface, $\omega_{t_{\tau}}$, is \mbox{\boldmathx}(t;\tau)=\mbox{\boldmath\varphi}(\mbox{\boldmathX}(\tau),t)=\mbox{\boldmathX}(\tau)+\mbox{\boldmathu}(\mbox{\boldmathX}(\tau),t), where $u$ is the displacement field. Note that $\omega_{t_{\tau}}\in\mathbb{R}^{3}$ also is a 2-manifold. The primary strain quantity is the deformation gradient, defined as \mbox{\boldmathF}(\mbox{\boldmathX}(\tau),t)=\partial\mbox{\boldmath\varphi}(\mbox{\boldmathX}(\tau),t)/\partial\mbox{\boldmathX}. Morpho-elastic growth of the soft mantle tissue is modeled by the multiplicative decomposition of the deformation gradient \mbox{\boldmathF}(\mbox{\boldmathX}(\tau),t)=\mbox{\boldmathF}^{\text{e}}(\mbox{\boldmathX}(\tau),t)\mbox{\boldmathF}^{\text{g}}(\mbox{\boldmathX}(\tau),t) into elastic and growth components, respectively  [33, 31, 20, 21].

The intuitive idea is that with the mid-surface of the shell being represented by $\mathscr{S}_{t}$, the mantle’s “preferred” state is given by the growth tensor \mbox{\boldmathF}^{\text{g}} relative to $\mathscr{S}_{t}\times(-h/2,h/2)$. Because of its attachment to the rigid shell the mantle cannot attain \mbox{\boldmathF}^{\text{g}}, but only $F$, with \mbox{\boldmathF}^{\text{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\geq\tau$ 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 \mbox{\boldmaths}_{2}, so that, accounting for the appropriate tangent spaces between which \mbox{\boldmathF}^{\text{g}} is imposed,

[TABLE]

Here, $\varepsilon_{2}$ is the rate of the growth strain along \mbox{\boldmaths}_{2}. As with the surface growth velocity, we will consider spatial and temporal variations in $\varepsilon_{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 outer surface.

The calcification front

While the mantle can be treated as an elastic solid, the calcified shell itself is rigid. An advancing calcification front, \mathscr{C}_{\tau}\in\mbox{\boldmath\varphi}(\mathscr{S}_{\tau}), is the interface between the rigid shell and material recently secreted by the mantle. The velocity of $\mathscr{C}_{\tau}$ is \mbox{\boldmathv}^{\text{c}}, which lies in the plane defined by \{\mbox{\boldmaths}_{1},\mbox{\boldmaths}_{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\in[0,T]$ is discretized by instants $t_{0},t_{1},\dots,t_{N}$, into sub-intervals $[t_{0},t_{1}],\dots,[t_{N-1},t_{N}]$, where $t_{0}=0$ and $t_{N}=T$. For simplicity, we also consider deposition times $\tau=t_{0},t_{1},\dots$. In a time step $\Delta t=t_{k+1}-t_{k}$, the leading surface of material secreted by the mantle, $\omega_{t_{t}}$, advances by v_{1}\Delta t\mbox{\boldmaths}_{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 volume 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 volume growth, (3) shell growth by secretion on the mantle’s current configuration, followed by (4) calcification. However, the discrete model operates with time steps $\Delta 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 above four processes are implemented in parallel over $[t_{k},t_{k+1}]$. We note that the time discretization reflects a time-discontinuous process of growth and calcification.

Surface growth

We assume that at time $t_{k}$, the shell has been fully calcified: \mathscr{C}_{t_{k}}=\mbox{\boldmath\varphi}(\Gamma_{t_{k}}). In the time interval $[t_{k},t_{k+1}]$, the leading surface is displaced due to mantle extension by v_{1}\Delta t\mbox{\boldmaths}_{1} from $\omega_{t_{t_{k}}}$ (expressed as a deformed configuration) to a new reference surface $\omega_{0_{t_{k+1}}}$. This allows us to define a strip of the mantle in its reference configuration $\Omega^{\text{m}}_{0_{t_{k}}}$ bounded by the surface $\omega_{t_{t_{k}}}$ at its trailing edge and $\omega_{0_{t_{k+1}}}$ at its leading edge. See Fig. 4.

Volume growth

Volume growth of the mantle is obtained by integrating Eqn (1). We exploit the exponential map:

[TABLE]

Since each increment of volume growth over a time step $\Delta t=t_{k+1}-t_{k}$ occurs relative to a new reference configuration, e.g., $\Omega^{\text{m}}_{0_{t_{k}}}$, we have \mbox{\boldmathF}^{\text{g}}_{k}=\mbox{\boldmath1}.

The actual deformation gradient achieved is \mbox{\boldmathF}_{k+1}, with the elastic component \mbox{\boldmathF}^{\text{e}}_{t_{k+1}}=\mbox{\boldmathF}_{t_{k+1}}\mbox{\boldmathF}_{t_{k+1}}^{\text{g}^{-1}} being governed by nonlinear elasticity. With a strain energy density function \psi(\mbox{\boldmathF}^{\text{e}}) that satisfies frame invariance (so that \psi(\mbox{\boldmathF}^{\text{e}})=\tilde{\psi}({\mbox{\boldmathF}^{\text{e}}}^{\text{T}}\mbox{\boldmathF}^{\text{e}})), the first Piola-Kirchhoff stress tensor is

[TABLE]

It is governed by the quasistatic stress equilibrium equation imposed at time $t_{k+1}$:

[TABLE]

In our computations we apply Dirichlet boundary conditions \mbox{\boldmathu}=\mbox{\boldmath0} on the trailing surface (boundary) $\omega_{t_{t_{k}}}$ of the mantle where it meets the rigid shell. A combination of fixed Dirichlet, \mbox{\boldmathu}=\mbox{\boldmath0}, and traction-free Neumann conditions, \mbox{\boldmathP}\mbox{\boldmathN}\big{|}_{t_{k+1}}=\mbox{\boldmath0} are applied on the remaining surfaces (boundaries) $\partial\Omega^{\text{m}}_{0_{t_{k+1}}}\backslash\omega_{t_{t_{k}}}$. This defines the morphoelastic growth problem for mapping the mantle strip from its reference configuration $\Omega^{\text{m}}_{0_{t_{k}}}$ to its deformed configuration $\Omega^{\text{m}}_{t_{t_{k}}}$.

Secretion

Following morphoelastic volume 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 $\Omega^{\text{m}}_{t_{t_{k}}}$.

Calcification

The final step of the algorithm is propagation of the calcification front so that \mathscr{C}_{t_{k+1}}=\mbox{\boldmath\varphi}(\Gamma_{t_{k+1}}). The mantle strip is calcified into its deformed configuration, $\Omega^{\text{m}}_{t_{t_{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 volume growth occur. The actual formation of new shell material by secretion over the mantle, and the calcification of this material, follow once the current mantle configuration has been defined by surface and volume growth.

Remark 2: The steps presented above impose full calcification of the secreted material in deformed configuration $\Omega^{\text{m}}_{t_{t_{k}}}$. Consequently, the mantle in its reference configuration $\Omega^{\text{m}}_{0_{t_{k+1}}}$ attaches to the rigid material surface $\omega_{t_{t_{k+1}}}$. An alternate model with possible biological relevance is to assume that $\Omega^{\text{m}}_{{t}_{t_{k}}}$ has not been calcified, but remains elastic. Then the attachment of the mantle reference configuration $\Omega^{\text{m}}_{0_{t_{k+1}}}$ at time $t_{k+1}$ and its morphoelastic volume growth over $[t_{k+1},t_{k+2}]$ further deforms $\Omega^{\text{m}}_{t_{t_{k}}}$, also. Equation (4) is then to be solved over $\Omega^{\text{m}}_{t_{t_{k}}}\cup\Omega^{\text{m}}_{0_{t_{k+1}}}$. Variants on this idea also are admissible, including complete calcification of a proper subset of $\Omega^{\text{m}}_{t_{t_{k}}}$ by time $t_{k+1}$, so that $\mathscr{C}_{t_{k+1}}$ does not coincide with \mbox{\boldmath\varphi}(\Gamma_{0_{t_{k+1}}}).

Implementation

The above formulation for surface and volume growth has been implemented in a finite element framework. An in-house C++ code based on the deal.II [34] open source finite element library is used to implement the model of surface and volume growth depicted in Figure  5. Key highlights of this computational implementation are the ability to handle growing meshes (to model the growth of the reference configuration in a Lagrangian setting) and related dynamic updates to the solution data structures (global vectors and matrices). Simulations presented in this work use hexahedral elements with a linear/quadratic Lagrange basis, and one to four layers of elements through the thickness. The code base is publicly available as a GitHub repository [35].

The attachment of the mantle to $\omega_{t_{t_{k}}}$ and its extension up to $\omega_{0_{t_{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 $\Omega^{\text{m}}_{t_{t_{k}}}$. Calcification is imposed by turning $\Omega^{\text{m}}_{t_{t_{k}}}$ rigid for time $t\geq 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 volume growth and the reference generating curve, $\Gamma_{0}$. In each case, the boundary condition is \mbox{\boldmathu}=\mbox{\boldmath0} on the trailing surface (boundary) $\omega_{t_{t_{k}}}$ of the mantle where it meets the rigid shell. The lateral surfaces (boundaries), $\partial\Omega^{\text{m}}_{0_{t_{k}}}\backslash\omega_{t_{t_{k}}}\backslash\omega_{0_{t_{k+1}}}$, are subject to Dirichlet conditions on either $u$, or its normal component, with traction-free Neumann conditions, \mbox{\boldmathP}\mbox{\boldmathN}=\mbox{\boldmath0} on the remaining surfaces (boundaries) $\omega_{0_{t_{k+1}}}$. See Fig. 5 for an illustration of the evolving mantle and calcified shell configurations.

Results

In the framework constructed above, there are three key parameters determining the ornamentation pattern that develops as the shell grows:

The active mantle width, i.e. the amount of surface growth occurring in each time increment, given by \delta s=|v_{1}\mbox{\boldmaths}_{1}-\mbox{\boldmathv}^{\text{c}}|\Delta t. 2. 2.

The volumetric growth increment over each time increment: $\delta g$, which is related to the growth strain rate \dot{\mbox{\boldmathF}}^{\text{g}}=\varepsilon_{2}\mbox{\boldmathF}^{\text{g}}\mbox{\boldmaths}_{2}\otimes\mbox{\boldmaths}_{2} by $\delta g=\varepsilon_{2}\Delta t$; 3. 3.

The initial curvature, $\kappa$, of $\Gamma_{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 the distribution of curvature along \mbox{\boldmaths}_{2} on the reference curve, $\Gamma_{0}$. The effects will be characterized by the mode number (the number of crests and valleys) along the lip of active mantle in the \mbox{\boldmaths}_{2} direction, the amplitudes and the locations of crests and valleys. 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.

Effect of parameters on morphology

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 generating curve, $\Gamma_{0}$, which is an arc of a circle with curvature $\kappa=0.01$. Fig. 7 shows the resulting morphologies when the active mantle width \delta s=|v_{1}\mbox{\boldmaths}_{1}-\mbox{\boldmathv}^{\text{c}}|\Delta t is varied by up to eightfold. We see that an increase in $\delta s$ leads to a decrease in mode number, which can be understood as follows: increasing $\delta s$ places the free edge of the active mantle strip, $\omega_{0_{t_{k+1}}}$, further from the rigid boundary $\omega_{t_{t_{k}}}$ ($k=1,2,\dots$) in the \mbox{\boldmaths}_{1}-direction, thus decreasing its structural stiffness to bending. The excess material (along \mbox{\boldmaths}_{2}) created by the increment in volume growth strain $\delta g$ can therefore be accommodated by an increased deformation in the \mbox{\boldmaths}_{3} direction, without paying a large strain energy penalty. As a result, each increment in $\delta s$ induces fewer wave crests/valleys, with progressively larger wavelength and amplitude. Fig. 7 also presents a comparison with the ornamentations on Clinocardium nuttallii and Tridacna squamosa, the active mantle width being much larger in the second species (giant clam), and both species differing in amplitude and wavelength of the antimarginal ribs in a manner that is consistent with the model predictions.

Here, it is instructive to compare to an analogous reduced order model: a growing 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 \mbox{\boldmaths}_{3} direction in our framework). This system has been formulated in detail by some of the authors [36] and forms the basis of previous mechanical descriptions of shell morphogenesis [24, 22]. In the linearized system, with the rod and foundation extending along the $x$-axis, the deformed rod has shape $y(x)$ satisfying [36]

[TABLE]

Here $\gamma>1$ describes the axial growth, analogous to $\delta 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\sim\exp(2\pi inx)$, the preferred bifurcation mode corresponds to the smallest value of $\gamma^{*}>1$ for which (5) has a solution; this is found to be $\gamma^{*}=1+2k+2(k+k^{2})^{1/2}$, from which we obtain that the mode number at buckling satisfies $n=\sqrt{\gamma^{*}-1}/(2\sqrt{2}\pi)$. From this we can extract the scaling relationship $n\sim\sqrt{k}$ for large $k$.

Based on the intuitive argument above, we would posit an inverse relationship between $k$ and $\delta s$, e.g. $k\sim\delta s^{\alpha}$ with $\alpha<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 $\delta 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, following crests or troughs, and ignoring the dependence of amplitude and wavelength on the coordinate in the \mbox{\boldmaths}_{2}-direction. To validate the computational model, we also include the critical mode as calculated from a buckling analysis of a plate (see S1 Text). From the slope of this log-log plot, we get $n\sim\delta s^{-1}$, and so $k\sim\delta s^{-2}$. In principle, one could use this map to more systematically parameterize the foundation in the 1D morphoelastic rod framework. Computing the morphology of the shell edge as a 1D (geometrically nonlinear) structure has the advantage of decreased computation time, though with potential inaccuracies due to loss of detail. A systematic means of determining $k$ provides a very useful step in alleviating this, though it remains an interesting question how far into the post-buckling regime the relation $k\sim\delta s^{-2}$ holds.

Morphoelastic volume 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 $\delta g$ along \mbox{\boldmaths}_{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 \mbox{\boldmaths}_{3} direction.

Fig. 9 illustrates the effect of varying the incremental growth strain $\delta g$ by up to a factor of $4\times$ while holding the incremental active mantle width fixed at $\delta s=1.0$. The reference generating curve $\Gamma_{0}$ has curvature of $\kappa=0.01$ in the middle of its extent in the \mbox{\boldmaths}_{2} direction, decreasing toward zero near the ends. Note that this induces a length scale: the radius of curvature $r=100$, relative to which $\delta s=1$.

In a time-continuous process, there would exist a critical time, $t_{\text{cr}}$, and corresponding amount of morphoelastic volume growth strain, $\delta g_{\text{cr}}=\int_{0}^{t_{\text{cr}}}\varepsilon_{2}(t)\mathrm{d}t$ 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 $\delta g_{\text{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 $\delta g=\delta g^{\ast}$) initiates bifurcation into a mode with three discernible crests, $n=3$. A growth strain increment to $\delta g=2\delta g^{\ast}$ 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 $\delta g=3\delta g^{\ast},4\delta g^{\ast}$ only increases the amplitude of the seventh mode, with no further bifurcations.

Remark 3: A careful examination of the bifurcated shape with mode number $n=3$ for $\delta g=\delta g^{\ast}$ reveals a deformation with amplitude that is highest at the midpoint of the arc of $\Gamma_{0}$, where the curvature $\kappa=0.01$, decreases in the two immediate neighbor crests, and decays by more than an order of magnitude toward the lateral edges, where $\kappa\to\infty$. Inclusion of all crests regardless of amplitude would raise the inferred mode 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 $\delta g=\delta g^{\ast}$, 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 $\delta g=2\delta g^{\ast}$ the strain energy settles into the seventh mode, whose prominence is magnified. In contrast, for $\kappa=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 $\kappa\to 0$ and uniform, where the amplitude remains uniform. The localization of mode shapes is a consequence of curvature, which we examine in greater detail in the next subsection on Curvature.

Curvature

Taking a cue from the localization of modes in high $\kappa$ regions of the reference generating curve, $\Gamma_{0}$, we next consider this effect in greater detail. We consider three geometries for $\Gamma_{0}$: a line with curvature $\kappa\to\infty$, a curve with $\kappa$ having almost uniform sign, and a second curve with $\kappa$ changing sign along the arc. The result appears in Fig 10. For all three cases in the figure, the Dirichlet boundary conditions are \mbox{\boldmathu}=\mbox{\boldmath0} on the trailing surface (boundary), $\omega_{0_{t_{0}}}$, and the lateral surfaces (boundaries) of the mantle, which are perpendicular to \mbox{\boldmaths}_{2}.

We find that the deformed configuration of the mantle $\Omega^{\text{m}}_{t_{t_{0}}}$ is biased toward developing curvature of the same sign as the reference curvature, $\kappa$, along \mbox{\boldmaths}_{2}. With $\Gamma_{0}$ a straight line in Fig 10a, the mantle deforms upward and downward equally, i.e. curvatures of both signs are seen in $\Omega^{\text{m}}_{t_{t_{0}}}$, and the crests/troughs have the same amplitude. The (mostly) single-signed $\kappa$ of Fig 10b promotes like-signed crests and suppresses oppositely signed ones. This is also apparent in Fig 10d (and to a lesser degree in Fig 10c), which has regions where $\kappa$ 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 $\delta g=\delta g^{\ast}$.

This result has an intriguing relevance for mollusk shell ornamentation: the reference shape of the shell on which ornamentation appears, modelled here by $\Gamma_{0}$, is generally convex and surrounds the mollusk body. 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 complex developmental process to ensure that the ornamentation is built in the outward direction. The results here suggest, rather, that the geometry and growth mechanisms naturally conspire to bias the pattern outward.

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, and 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 Fig 10(c), having varying curvature, and impose volume growth strain increments ranging from $\delta g=0.0$ to $\delta 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. The pronounced localization into spines has been explained by some of the authors of this communication via the added mechanism of spatially varying material properties [36]. Fig 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 also replicates the trend of favoring compatible crests.

Hierarchical 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 presented in the earlier sections 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 could explain the secondary and tertiary crests and valleys that are visible along the shell edge in species such as Hexaplex cichoreum. We recognize these as the potential remnants of increasing mode number during shell growth. Noting that the decrease in mode number with increasing active mantle width, as shown earlier, implies an increase in mode number with decreasing active mantle width, and recalling the magnification of higher modes with increasing volume growth strain rates, we consider the following protocol of surface and volume growth rates in Fig 13: initial surface growth lays down a mantle of width $\delta s_{0}=\delta s^{\ast}$, forming a reference configuration $\Omega^{\text{m}}_{0_{t_{0}}}$, which, under morphoelastic volume growth over $(t_{0},t_{1}]$ deforms into $\Omega^{\text{m}}_{t_{t_{0}}}$. Upon calcification, the curved mantle edge $\omega_{t_{t_{1}}}$ provides the reference curvature for subsequent growth. The second reference configuration $\Omega^{\text{m}}_{0_{t_{1}}}$ laid down by surface growth has only half the initial mantle width: $\delta s_{1}=0.5\delta s^{\ast}$. However, with the same morphoelastic volume growth $\delta g_{1}=\delta g^{\ast}$ over $(t_{1},t_{2}]$, a mode of higher mode number ($n=3$) develops: secondary ornamentation is achieved on configuration $\Omega^{\text{m}}_{t_{t_{1}}}$. No further surface growth occurs, but the volume growth undergoes another increment of $\delta g_{2}=2\delta g^{\ast}$ over $(t_{2},t_{3}]$. Another bifurcation into a higher mode, $n=4$, is seen and further detail of secondary ornamentation is visible on the configuration $\Omega^{\text{m}}_{t_{t_{2}}}$. In the interval $(t1,t_{3}]$, the calcification front is stationary at a distance of $\delta s^{\ast}$ along the nominal \mbox{\boldmaths}_{1} direction, as indicated by the dotted white line in Fig 13.

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 bifurcation 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 computational complexity associated with ensuring curvature continuity in the \mbox{\boldmaths}_{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 Gaussian 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 valleys form with positive Gaussian curvature. One exception is Fig. 10b. This case is explained by the strain energy due to high curvature, $\kappa$ of the reference curve, $\Gamma_{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 $\Omega^{\text{m}}_{t_{t_{1}}}$, changes into negative Gaussian curvature after the final volumetric growth increment in $\Omega^{\text{m}}_{t_{t_{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 $\delta s^{\ast}$ formed by surface growth over $(t_{0},t_{1}]$ experiences volume growth $\delta g=\delta g^{\ast}$, but the strip of length $\delta s=\frac{1}{2}\delta s^{\ast}$ from surface growth over $(t_{1},t_{2}]$ experiences a total volume growth of $\delta g=3\delta g^{\ast}$. We are therefore prompted to consider that the rate of growth strain $\varepsilon_{2}$ is an increasing function along the \mbox{\boldmaths}_{1} direction, i.e., $\partial\varepsilon_{2}/\partial\xi_{1}>0$, where $\xi_{1}$ is the curvilinear coordinate defining \mbox{\boldmaths}_{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 $\varepsilon_{2}(\xi_{1})$ has low but positive curvature $\partial^{2}\varepsilon_{2}/\partial\xi_{1}^{2}>0$ for small $\xi_{1}$, changing smoothly to high curvature $\partial^{2}\varepsilon_{2}/\partial\xi_{1}^{2}\gg 0$ for larger $\xi_{1}$. The variation in the rate of growth strain generates a deformation with a finite component in the negative \mbox{\boldmaths}_{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 of bivalves, cephalopods and gastropods 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(b). Here we have another instance of a mechanical basis for a feature of ornamentation in mollusk shells 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

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 [38] 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 [33, 39]. 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 [40, 41, 42], aspects of cardiovascular systems, and the folding of soft, layered structures during morphogenesis [43], 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, volume 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 \mbox{\boldmaths}_{1} and \mbox{\boldmathv}^{\text{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 listed above for immersed boundary type methods could be approached by nonuniform element sizes in the advancing surface, and stepped fronts approximating changing directions on average.

Most crucially, our work has identified the prominent role played by geometry in controlling mantle deformation under the driver of morphoelastic growth. The active mantle width, which is a direct outcome of the surface growth rate, has a very visible influence on the mode number, mode shape and, as we have demonstrated, on the appearance of hierarchical ornamentation. The curvature of reference surfaces is the other parameter by which geometry acts directly to control the locations of crests and valleys. The evolving reference configurations, as each generation of active mantle is enslaved to the reference curvature of the previous generation, present a pathway for coupling of morphoelastic volume growth with surface growth and curvature. Figure 15 is a “phase diagram” illustrating salient aspects of our studies over this parameter space.

We have not attempted to compile computational demonstrations that match molluskan morphologies with high fidelity. While, in our opinion, most of our computations match well with features of actual molluskan morphologies, others such as in Figures 10c and 10d are less satisfying, especially in representation of spiny outgrowths (see the following paragraph in this regard). It is also true, however, that more complexity could be introduced to the model. Contact mechanics is one such addition, which is on our critical path, but must await a future communication. Another is the spatial variation of material properties, which we have already addressed before (see following paragraph). While a detailed tuning over such effects may add insight to mechanisms, several trends are visible in what we have explored here. We have experimented with some variations on our basic themes; variations that typically are not describable with analytic forms, but are ubiquitous in nature, thus making them obvious candidates for computational exploration. This is the rational for investigations of temporal and spatial variations in growth.

We note that previous work by some of the authors [36] already has highlighted the role of a spatial variation in material properties in shaping the sharp spines seen on the shells of bivalves, cephalopods and gastropods. The present communication adds to this emerging picture of the influence of physics, by shifting the focus to geometry, which, driven by morphoelastic growth, acts through the mechanism of surface growth and a parametric dependence on curvature.

Of particular interest is the further coupling of this framework with pattern formation by a range of reaction-transport equations. This would make accessible well-studied developmental milestones such as the patterning and morphogenesis of limbs and digits. Other problems in morphogenesis, such as the formation of skeletons are also within reach of our computational framework. Finally, we note that we have modelled antimarginal ornamentation events via bursts of growth. It is more likely that the growth rate does not change significantly at the location of ornamentations, but rather that the thickness of the mantle decreases, creating an increase in the length of the active mantle strip and a corresponding decrease in stiffness without requiring a strong increase in material. This effect would be interesting to incorporate in future studies, as the decrease in stiffness would further amplify the amplitudes of transverse deformation, and also may be non-uniform along the length of the mantle strip.

S1 Text. Buckling analysis of a plate

We compute here the relationship between buckling mode and active mantle width via the buckling of a plate. We consider a plate with zero reference curvature of length $A$ in the $x$ direction and width $B$ in the $y$ direction (the $x$ and $y$ directions here correspond to the $s_{2}$ and $s_{1}$ directions in the main text, respectively, so we are primarily interested in the case $A\gg B$). The governing equation for the transverse deformation $w(x,y)$ of the plate is [2, 1]

[TABLE]

where $D$ is the bending modulus and $N$ is a compressive force due to growth in the $x$ direction (defined as positive here). For boundary conditions, we take the plate to be clamped on one long edge, free on the other long edge, and simply supported on the two short edges. These conditions read

[TABLE]

where $\nu$ is the Poisson ratio. The system (S1), (S2) has solution

[TABLE]

where $m$ is the buckling mode. Taking $f=\exp(i\lambda y)$ yields the characteristic equation

[TABLE]

where $\tilde{N}=N/D$. This has roots

[TABLE]

a nontrivial solution only existing if $\tilde{N}>m^{2}\pi^{2}/A^{2}$. The function $f(y)$ thus takes the form

[TABLE]

The boundary conditions in the $y$-direction translate to

[TABLE]

Imposing these conditions yields an eigenvalue problem for the critical compression $\tilde{N}=\tilde{N}^{*}$, and the critical buckling mode is determined by finding the integer value of $m$ at which $\tilde{N}^{*}$ is minimized. The points in Fig 8 of the main text were produced by fixing $A=20$, $\nu=0.3$, varying $B$ between 0.5 and 2.5, and computing the critical mode at each width.