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

Shiva Rudraraju, Derek Moulton, Régis Chirat, Alain Goriely, Krishna Garikipati

<figure><img src="image_1.png" /></figure>

## To cite this version:

Shiva Rudraraju, Derek Moulton, Régis Chirat, Alain Goriely, Krishna Garikipati. A computational framework for the morpho-elastic development of molluskan shells by surface and volume growth. PLoS Computational Biology, 2019, 15 (7), pp.e1007213. 10.1371/journal.pcbi.1007213. hal-04843926

HAL Id: hal-04843926

https://hal.science/hal-04843926v1

Submitted on 18 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

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

Shiva Rudraraju<sup>1</sup>, Derek E. Moulton<sup>2</sup>, Régis Chirat<sup>3</sup>, Alain Goriely<sup>1,2</sup>,
Krishna Garikipati<sup>4*</sup>

1 Department of Mechanical Engineering, University of Wisconsin-Madison, Madison, Wisconsin, United States of America, 2 Mathematical Institute, University of Oxford, Oxford, United Kingdom, 3 UMR CNRS 5276 LGL-TPE, Université Lyon1, 69622 Villeurbanne Cedex, France, 4 Departments of Mechanical Engineering and Mathematics, Michigan Institute for Computational Discovery & Engineering, University of Michigan, Ann Arbor, Michigan, United States of America

* krishna@umich.edu

<figure><img src="image_3.png" /></figure>

## 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.

## Author Summary

Molluska are the second most diversified phylum of the animal kingdom, and their evolutionary success can be partly attributed to the hard shell that provides both protection and support to the soft body. The distinctive anatomical features of these hard shells are their rich pigmentation patterns and complex structural ornamentations. While the pigmentation patterns are primarily of biochemical origin, the ornamentations result from mechanical deformation of the mantle due to growth induced forces. This mechanical basis of “growth and form” has been previously investigated using simplified morpho-

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Rudraraju S, Moulton DE, Chirat R, Goriely A, Garikipati K (2019) A computational framework for the morpho-elastic development of molluskan shells by surface and volume growth. PLoS Comput Biol 15(7): e1007213. https://doi.org/10.1371/journal.pcbi.1007213

**Editor:** Qing Nie, University of California Irvine, UNITED STATES

**Received:** January 31, 2019

**Accepted:** June 27, 2019

**Published:** July 29, 2019

**Copyright:** © 2019 Rudraraju et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting Information files.

**Funding:** The authors received no specific funding for this work.

**Competing interests:** The authors have declared that no competing interests exist.
mechanical models, but restricted to reduced geometric representations. Here we propose a three-dimensional computational framework coupling morphology, incremental surface growth by accretion, and morpho-elastic volume growth, to enable an improved representation of the growth and structural parameters controlling the evolution of these ornamentations. We study the effects of growth rates, and of previously deposited shell geometries on the resulting modes of mantle deformation, and present a “phase diagram” of morphogenesis in molluskan shells. 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.

## 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–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–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. [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.

## 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
<figure><img src="image_5.png" /><figcaption>**Fig 1. Schematic of shell growth.** Growth process of a molluscan shell surface depicted through the steps of volume growth of mantle tissue, morphoelastic deformation, and shell surface growth via secretion and calcification. The calcified region of the shell is indicated in yellow.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g001

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, $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 over its pre-existing extent (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 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 $s_3$. We have $s_3 = s_1 \times s_2$, and more specifically, $s_i \cdot s_j = \delta_{ij}$, where the triad $\{s_1, s_2, s_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 [20, 21, 29–32]. 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, $\mathcal{S}_{\tau} \subset \mathbb{R}^3$. We regard $\mathcal{S}_{\tau}$ as a one-parameter family of surfaces, generated by $\tau \in [0, T]$, from a reference surface $\mathcal{S}_0 \subset \mathbb{R}^3$. The generating curve, $\Gamma_{\tau} \subset \partial\mathcal{S}_{\tau}$, is the leading edge of $\mathcal{S}_{\tau}$ and evolves along $s_1$ (see Fig 2). For points $X(\tau) \in \Gamma_{\tau}$ and $X(0) \in \Gamma_0$, where $\Gamma_0 \subset \partial\mathcal{S}_0$ is the initial generating curve at time $\tau = 0$, we have $X(\tau) = \chi_{\tau}(X(0))$. Surface growth occurs by extension of the mantle along the boundary curve $\Gamma_{\tau}$, which advances with the velocity $\dot{\chi} = v_1s_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
<figure><img src="image_6.png" /><figcaption>**Fig 2.** Local coordinate system on the surface patch. *s*₁ is the direction of surface growth due to mantle extension, *s*₂ is the direction of volume growth along the mantle margin, and *s*₃ = *s*₁ × *s*₂ is the normal to the local surface patch. The calcified shell edge after time τ₁ forms the generating curve Γτ₁ for the time step, [τ₁, τ₂]; the leading edge of the grown and deformed mantle strip then forms the generating curve, Γτ₂ for the next time step.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g002
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 $\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 $s_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, $X(\tau) \in \omega_{0\tau}$ if it was deposited at time $\tau$. The point-to-point map of material points $X(\tau)$ from the reference surface $\omega_{0\tau}$ to $x(t; \tau)$ on the current surface, $\omega_{t\tau}$, is $x(t; \tau) = \varphi(X(\tau), t) = X(\tau) + u(X(\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 $F(X(\tau), t) = \partial\varphi(X(\tau), t)/\partial X$. Morpho-elastic growth of the soft mantle tissue is modeled by the multiplicative decomposition of the deformation gradient $F(X(\tau), t) = F^e(X(\tau), t)F^g(X(\tau), t)$ into elastic and growth components, respectively [20, 21, 31, 33].

The intuitive idea is that with the mid-surface of the shell being represented by $\mathcal{S}_t$, the mantle’s “preferred” state is given by the growth tensor $F^g$ relative to $\mathcal{S}_t \times (-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 \ge \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 $s_2$, so that, accounting for the appropriate tangent spaces between which $F^g$ is imposed,

$$ \dot{F}^g = \varepsilon_2 F^g s_2 \otimes s_2. \qquad (1) $$

Here, $\varepsilon_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 $\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, $\mathcal{C}_{\tau} \in \varphi(\mathcal{S}_{\tau})$, is the interface between the rigid shell and material recently secreted by the mantle. The velocity of $\mathcal{C}_{\tau}$ is $v^e$, 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 \in [0, T]$ is discretized by instants $t_0, t_1, \dots, t_N$, into sub-
<figure><img src="image_7.png" /><figcaption>**Fig 3. The kinematics of growth.** The observed deformation gradient, *F*, is composed of an incompatible growth component, *F*<sup>g</sup>, and another incompatible, but elastic component, *F*<sup>e</sup>, which restores compatibility of *F*.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g003

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_i}$, advances by $v_1 \Delta t s_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: $C_{t_k} = \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 s_1$ from $\omega_{t_{t_k}}$ (expressed as a deformed configuration) to a new reference surface $\omega_{t_{0|t_{k+1}}}$. This allows
<figure><img src="image_8.png" /><figcaption>**Fig 4.** Mathematical model of molluscan shell growth over (t<sub>k</sub>, t<sub>k+1</sub>] through a sequence of surface growth, volume growth, secretion and calcification. Beginning with a calcified mantle current configuration Ω<sub>0t<sub>k+1</sub></sub><sup>m</sup>, growth over (t<sub>k</sub>, t<sub>k+1</sub>] is modelled via the following sequence of steps: (1) Surface growth—leading surface ω<sub>t<sub>k</sub></sub> is displaced due to mantle extension by v<sub>1</sub>Δt s<sub>1</sub> to a new reference surface ω<sub>0t<sub>k+1</sub></sub><sup>m</sup>, defining a strip of the mantle in its reference configuration Ω<sub>0t<sub>k+1</sub></sub><sup>m</sup>. (2) Volume (morphoelastic) growth— Ḟ<sup>g</sup> is imposed on the mantle strip Ω<sub>0t<sub>k+1</sub></sub><sup>m</sup> over Δt = t<sub>k+1</sub> - t<sub>k</sub>, driving its nonlinear deformation into the current configuration Ω<sub>t<sub>k+1</sub></sub><sup>m</sup>. (3) Shell growth occurs by secretion on the mantle strip Ω<sub>t<sub>k+1</sub></sub><sup>m</sup>, followed by (4) calcification of the mantle strip Ω<sub>t<sub>k+1</sub></sub><sup>m</sup> at t<sub>k+1</sub>. During the volume growth of the mantle strip from its reference configuration Ω<sub>0t<sub>k+1</sub></sub><sup>m</sup> to its current configuration Ω<sub>t<sub>k+1</sub></sub><sup>m</sup>, boundary conditions are applied on the trailing surface (highlighted in red) and the lateral surfaces (highlighted in yellow) of the mantle strip. The front surface of the mantle strip is highlighted in green.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g004

us to define a strip of the mantle in its reference configuration Ω<sub>0t<sub>k+1</sub></sub><sup>m</sup> bounded by the surface
ω<sub>t<sub>k</sub></sub> at its trailing edge and ω<sub>0t<sub>k+1</sub></sub> at its leading edge. See Fig 4.

**Volume growth.** Volume growth of the mantle is obtained by integrating Eq 1. We exploit the exponential map:

$$
F_{t_{k+1}}^{g} = F_{k}^{g} \exp[\varepsilon_{2k} \Delta t s_{2} \otimes s_{2}]. \quad (2)
$$

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_{0t_{k+1}}^m$, we have $F_k^g = 1$.

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

$$
\mathbf{P} = \partial\psi / \partial\mathbf{F}^e. \qquad (3)
$$
It is governed by the quasistatic stress equilibrium equation imposed at time $t_{k+1}$:

$$ \mathrm{Div} \, P_{t_{k+1}} = 0 \quad \text{in} \quad \Omega_{t_{k+1}}^m. \qquad (4) $$

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

**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_{t_{k+1}}^m$.

**Calcification.** The final step of the algorithm is propagation of the calcification front so that $\mathcal{C}_{t_{k+1}} = \varphi(\Gamma_{t_{k+1}})$. The mantle strip is calcified into its deformed configuration, $\Omega_{t_{k+1}}^m$.

**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_{t_{k+1}}^m$. Consequently, the mantle in its reference configuration $\Omega_{0_{t_{k+1}}}^m$ attaches to the rigid material surface $\omega_{t_{k+1}}$. An alternate model with possible biological relevance is to assume that $\Omega_{t_{k+1}}^m$ has not been calcified, but remains elastic. Then the attachment of the mantle reference configuration $\Omega_{0_{t_{k+1}}}^m$ at time $t_{k+1}$ and its morphoelastic volume growth over $[t_{k+1}, t_{k+2}]$ further deforms $\Omega_{t_{k+1}}^m$, also. Eq (4) is then to be solved over $\Omega_{t_{k+1}}^m \cup \Omega_{0_{t_{k+1}}}^m$. Variants on this idea also are admissible, including complete calcification of a proper subset of $\Omega_{t_{k+1}}^m$ by time $t_{k+1}$, so that $\mathcal{C}_{t_{k+1}}$ does not coincide with $\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 Fig 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_{k+1}}$ 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 Eq (2) subject to the constitutive law Eq (3) and the governing Eq (4). Following Remarks 1 and 2, secretion of shell material is a virtual step over the deformed mantle configuration $\Omega_{t_{k+1}}^m$. Calcification is imposed by turning $\Omega_{t_{k+1}}^m$ rigid for time $t \ge t_k$. 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 $u = 0$ on the trailing surface (boundary) $\omega_{t_{k+1}}$ of the mantle where it meets the rigid shell. The lateral surfaces (boundaries), $\partial\Omega_{0_{t_{k+1}}}^m \setminus \omega_{t_{k+1}} \setminus \omega_{0_{t_{k+1}}}$, are subject to Dirichlet conditions
<figure><img src="image_9.png" /><figcaption>**Fig 5.** Space-time discretization in the finite element framework. Evolution of the mollusk shell surface through surface growth and morphoelastic volume growth of mantle strips followed by their calcification. Also see S1 Movie for the time evolution of a representative molluskan shell surface through the accretive growth of 20 mantle strips.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g005

on either **u**, or its normal component, with traction-free Neumann conditions, **PN** = **0** on the remaining surfaces (boundaries) $\omega_{0_{ik+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:

1. The *active mantle width*, i.e. the amount of surface growth occurring in each time increment, given by $\delta s = |v_1 s_1 - v^c| \Delta t$.

2. The *volumetric growth increment* over each time increment: $\delta g$, which is related to the growth strain rate $\dot{F}^g = \varepsilon_2 F^g s_2 \otimes s_2$ by $\delta g = \varepsilon_2 \Delta t$;

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
<figure><img src="image_10.png" /><figcaption>**Fig 6.** Parameters controlling the morphology of shell ornamentations. The reference generating curve, Γ₀, with its curvature, κ along *s*₂; the active mantle width, *δs*; volume growth strain increment, *δg*.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g006

them on a mechanistic basis while making connections to ornamentations observed on mol-
luskan 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 *s*₂ on the reference curve, Γ₀. The effects will be characterized by the mode number (the number of crests and valleys) along the lip of active mantle in the *s*₂ 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, Γ₀, which is an arc of a circle with curvature κ = 0.01. Fig 7 shows the resulting morphologies when the active mantle width δs = |v₁ s₁ – vᶜ|Δt is varied by up to eightfold. We see that an increase in δs leads to a decrease in mode number, which can be understood as follows: increasing δs places the free edge of the active mantle strip, ω₀ₜₖ₊₁, further from the rigid boundary ωₜₖ (k = 1,2,...)
<figure><img src="image_11.png" /><figcaption>**Fig 7.** Effect of incremental active mantle widths on mantle deformation for fixed $\delta g$ and $\kappa$. Increasing the active mantle width $\delta s = |v_1 s_1 - v^*|\Delta t$ over $\Gamma_0$ leads to morphologies bearing a similarity with the ornamentations on *Clinocardium nuttallii* (upper inset) and *Tridacna squamosa* (lower inset). Dirichlet boundary conditions, $u = 0$, are applied on the trailing surface (boundary) and traction-free Neumann conditions, $PN = 0$, are applied on the remaining surfaces (boundaries). See Fig 4 for location of the trailing surface, front surface and the lateral surfaces. Also see S3 Movie for a morphology that is similar to the case $\delta s = \delta s^*$, and bears comparison to the ornamentation on members of the class *Bivalvia*.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g007

in the $s_1$-direction, thus decreasing its structural stiffness to bending. The excess material (along $s_2$) created by the increment in volume growth strain $\delta 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 $\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 $s_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 [22, 24]. In the linearized system, with the rod and foundation extending along the x-axis, the deformed rod has shape $y(x)$ satisfying [36]

$$y'''(\mathbf{X}) + (\gamma - 1)y''(\mathbf{X}) + k\gamma y(\mathbf{X}) = 0, \quad (5)$$
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 i n x)$, 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 $s_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 $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 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 $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 $\delta s = 1$.

In a time-continuous process, there would exist a critical time, $t_{cr}$, and corresponding amount of morphoelastic volume growth strain, $\delta g_{cr} = \int_0^{t_{cr}} \varepsilon_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 $\delta 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 $\delta g = \delta g^*$) initiates bifurcation into a mode with three discernible crests, $n = 3$. A growth strain increment to $\delta g = 2\delta 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 $\delta g = 3\delta g^*$, $4\delta g^*$ only increases the amplitude of the seventh mode, with no further bifurcations.
<figure><img src="image_12.png" /><figcaption>**Fig 8.** Scaling study of mode number *n* versus the incremental active mantle width *δs*. Shown is the dependence of the mode number of the deformed mantle strip on the incremental active mantle width, obtained from the finite element model (FE) and the buckling analysis of a plate.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g008

**Remark 3:** A careful examination of the bifurcated shape with mode number *n* = 3 for *δg* = *δg*\* reveals a deformation with amplitude that is highest at the midpoint of the arc of Γ₀, 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 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 *δ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,
<figure><img src="image_13.png" /><figcaption>**Fig 9.** Effect of volume growth strain increment. $\delta g = \varepsilon_2 \Delta t$ on the amplitude and mode numbers of the deformed mantle. The observed high mode number morphologies are similar to the ornamentations observed in bivalves like *Clinocardium nuttallii* (inset). Dirichlet boundary conditions, $\mathbf{u} = \mathbf{0}$, are applied on the trailing surface (boundary) and traction-free Neumann conditions, $\mathbf{PN} = \mathbf{0}$, are applied on the remaining surfaces (boundaries). See Fig 4 for location of the trailing surface, front surface and the lateral surfaces.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g009

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 $\mathbf{u} = \mathbf{0}$ on the trailing surface (boundary), $\omega_{0\nu_0}$, and the lateral surfaces (boundaries) of the mantle, which are perpendicular to $\mathbf{s}_2$.

We find that the deformed configuration of the mantle $\Omega_{\nu_0}^m$ is biased toward developing curvature of the same sign as the reference curvature, $\kappa$, along $\mathbf{s}_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_{\nu_0}^m$, 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^*$.
<figure><img src="image_14.png" /><figcaption>**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 generating curve. These reference curves bear similarity to the shape of the mantle surface (highlighted in red) found in (a) *Pterynotus phyllopterus*, (b) *Nucella freycineti*, (c) *Normal Bolinus brandaris* and (d) *Abnormal Bolinus brandaris* (insets). Dirichlet boundary conditions, *u* = 0, are applied on the trailing surface (boundary) and the lateral surfaces (boundaries) of the mantle, which are perpendicular to *s*2, and traction-free Neumann conditions, *PN* = 0, are applied on the front surface (boundary). The underlying spatial discretization (mesh) is also shown on the model geometries. Also see S2–S5 Movies for the evolution of mantle deformation and accretive growth over planar, arc, and closed circular geometries of the reference curves.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g010

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^*$, forming a reference configuration $\Omega_{t_{0}}^{m}$, which, under morphoelastic volume growth over ($t_0, t_1$] deforms into $\Omega_{t_{0}}^{m}$. Upon calcification, the curved mantle edge $\omega_{t_1}$ provides the reference curvature for subsequent growth. The second reference configuration $\Omega_{t_1}^{m}$ laid down by surface growth has only half the initial mantle width: $\delta s_1 = 0.5\delta s^*$. However, with the same morphoelastic volume growth $\delta g_1 = \delta g^*$ over ($t_1, t_2$], a mode of higher mode number ($n=3$) develops: secondary ornamentation is achieved on configuration $\Omega_{t_1}^{m}$. No further surface growth occurs, but the volume growth undergoes another increment of $\delta g_2 = 2\delta g^*$ 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_{t_2}^{m}$. In the interval ($t_1, t_3$], the calcification front is stationary at a distance of $\delta s^*$ along the nominal $s_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
<figure><img src="image_15.png" /><figcaption>**Fig 11.** Progression of curvature-compatible ornamentation with volume growth strain increments, $\delta g$. The deformed mantles show marginally preferred localization around points of high curvature and thereafter the amplitude increases with volume growth strain increments. Some of these shapes with different amplitudes can be observed in the shells of the species *Bolinus brandaris* (bottom row of inset images). Dirichlet boundary conditions, $u = 0$, are applied on the trailing surface (boundary) and traction-free Neumann conditions, $PN = 0$, are applied on the remaining surfaces (boundaries).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g011

resolving details beyond secondary modes becomes challenging due to the computational 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.
<figure><img src="image_16.png" /><figcaption>**Fig 12.** Influence of reference curvature singularities, and of smooth curvatures. 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. Dirichlet boundary conditions, *u* = 0, are applied on the trailing surface (boundary) and traction-free Neumann conditions, *PN* = 0, are applied on the remaining surfaces (boundaries).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g012

**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_{t_1}^m$, changes into negative Gaussian curvature after the final volumetric growth increment in $\Omega_{t_2}^m$. 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^*$ formed by surface growth over ($t_0, t_1$] experiences volume growth $\delta g = \delta g^*$, but the strip of length $\delta s = \frac{1}{2} \delta s^*$ from surface growth over ($t_1, t_2$] experiences a total volume growth of $\delta g = 3\delta g^*$. We are therefore prompted to consider that the rate of growth strain $\epsilon_2$ is an increasing function along the $s_1$ direction, i.e., $\partial\epsilon_2/\partial\xi_1 > 0$, where $\xi_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.
<figure><img src="image_17.png" /><figcaption>**Fig 13.** Hierarchical ornamentation arising from temporally varying surface growth, *δs* volume growth strains, *δg* *ε*₂. 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 computational model. The dotted white line indicates the location of the fixed calcification front between (*t*₁, *t*₃]. Dirichlet boundary conditions, *u* = 0, are applied on the trailing surface (boundary) and the lateral surfaces (boundaries) of the mantle, which are perpendicular to *s*₂, and traction-free Neumann conditions, *PN* = 0, are applied on the front surface (boundary). Inset image of *Hexaplex cichoreum* modified from source [37]. Original images licensed under the Creative Commons Attribution-Share Alike License.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g013

An example of directly imposing such spatial variation is shown in Fig 14, where the profile of ε₂(ξ₁) has low but positive curvature ∂²ε₂/∂ξ² > 0 for small ξ₁, changing smoothly to high curvature ∂²ε₂/∂ξ² ≫ 0 for larger ξ₁. The variation in the rate of growth strain generates a deformation with a finite component in the negative *s*₁ direction, i.e. the mantle “arches back” to accommodate the excess length, creating ornamentation with negative Gaussian curvature.
<figure><img src="image_18.png" /><figcaption>**Fig 14. Spatially varying volume growth.** We impose volume growth strain increments that vary along the $ξ_1$ direction that is the curvilinear coordinate defining $s_1$, 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 (c) *Hexaplex chicoreum*. Dirichlet boundary conditions, $u = 0$, are applied on the trailing surface (boundary) and the lateral surfaces (boundaries) of the mantle, which are perpendicular to $s_1$, and traction-free Neumann conditions, $PN = 0$, are applied on the front surface (boundary). Also see S6 Movie for the evolution of deformation leading to a morphology with strongly negative Gaussian curvature.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g014

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–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 $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 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. Fig 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 Fig 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.
<figure><img src="image_19.png" /><figcaption>Fig 15. Phase diagram representing the effect of the growth and geometric parameters—growth strain increment ($\delta_g$), active mantle width ($\delta_g$), and curvature of the reference curve ($\kappa$), on the morphology of shell ornamentations.</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1007213.g015

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.

## Supporting information

**S1 Text. Buckling analysis of a plate.**

(PDF)

**S1 Movie. Space-time discretization of the molluskan shell in the finite element framework.** Shown at the top is the evolution of the geometry (finite element mesh) for 20 surface growth increments (addition of the mantle in 20 incremental strips, each being four elements wide) of a representative molluskan shell in its reference configuration, without the morphoelastic volume growth increments. As a result the reference configuration remains a flat plate. The accompanying computation at the bottom shows the 20 growth increments with the complete model (surface growth, volume growth and calcification of 20 mantle strips in sequence). Beginning with a flat plate geometry, each surface growth increment is followed by its morphoelastic volume growth increment and calcification. Calcification is the final stage of the sequence for each such mantle strip of four elements. Therefore, at any instant, it is only the mantle strip at the leading edge that undergoes morphoelastic volume growth. The mantle strips that grew before it have already undergone calcification. This is the case for S1 Movie-S3 and S5 Movies.

(MP4)

**S2 Movie. Influence of the geometry of reference curves on antimarginal ornamentation: Planar geometry.** Shown is the evolution of 10 mantle strips starting from a flat plate geometry of the reference curve. The contour colors indicate the normalized displacement magnitude. As in S1 Movie, at any instant, it is only the mantle strip at the leading edge that undergoes morphoelastic volume growth. The mantle strips that grew before it have already undergone calcification.

(MP4)

**S3 Movie. Influence of the geometry of reference curves on antimarginal ornamentation: Arc geometry.** Shown is the evolution of 4 mantle strips starting from an arc geometry of the reference curve. The contour colors indicate the normalized displacement magnitude. As in
S1 and S2 Movies, at any instant, it is only the mantle strip at the leading edge that undergoes morphoelastic volume growth. The mantle strips that grew before it have already undergone calcification.

(MP4)

**S4 Movie. Influence of the geometry of reference curves on antimarginal ornamentation: An arc with positive and negative curvature.** Shown is the evolution of a single mantle strip starting from a reference curve that is an arc with curvature of changing signs. Note that there is no surface growth in this movie. The trailing edge is the only calcified part of the shell. The snap-through events seen during the deformation of the mantle strip are the elastic bifurcations (buckling modes) triggered by growth over this geometry. There are several bifurcations occurring in rapid succession due to volume growth within a single growth increment, in this simulation. Because of the stiff, nonlinear response of the shell undergoing elastic bifurcations, the single increment of volume growth is numerically applied as 20 load steps, and the evolution of the deformed geometry after each load step is shown in the corresponding evolution of the 20 frames shown in this movie. The contour colors indicate the normalized displacement magnitude.

(MP4)

**S5 Movie. Influence of the geometry of reference curves on antimarginal ornamentation: Closed circular geometry.** Shown is the evolution of 3 mantle strips starting from a circular geometry of the reference curve, and is representative of growth over a closed shell geometry. The contour colors indicate the normalized displacement magnitude. As in S1–S3 Movies, at any instant, it is only the mantle strip at the leading edge that undergoes morphoelastic volume growth. The mantle strips that grew before it have already undergone calcification.

(MP4)

**S6 Movie. Evolution of deformation leading to a backward arching morphology.** The contour colors indicate the normalized displacement magnitude. Note that there is no surface growth in this movie. The trailing edge is the only calcified part of the shell.

(MP4)

## Author Contributions

**Conceptualization:** Shiva Rudraraju, Derek E. Moulton, Alain Goriely, Krishna Garikipati.

**Formal analysis:** Derek E. Moulton, Krishna Garikipati.

**Investigation:** Derek E. Moulton, Alain Goriely.

**Methodology:** Shiva Rudraraju.

**Software:** Shiva Rudraraju.

**Validation:** Régis Chirat, Krishna Garikipati.

**Writing – original draft:** Shiva Rudraraju, Derek E. Moulton, Alain Goriely, Krishna Garikipati.

**Writing – review & editing:** Shiva Rudraraju, Derek E. Moulton, Régis Chirat, Alain Goriely, Krishna Garikipati.

## References

1. Ponder W. F. and Linderg D. R. Molluscan evolution and phylogeny. In *Phylogeny and Evolution of the Mollusca*, pages 1–17. University of California Press, Berkeley, 2008.
2. Wanninger A. *Evolutionary Developmental Biology of Invertebrates 2: Lophotrochozoa (Spiralia)*. Springer, 2015.

3. Sigwart J. D. Zoology: molluscs all beneath the sun, one shell, two shells, more, or none. *Current Biology*, 27(14):R708–R710, 2017. https://doi.org/10.1016/j.cub.2017.05.075 PMID: 28743018

4. Runnegar B. Early evolution of the mollusca: the fossil record. *Origin and evolutionary radiation of the Mollusca*, 1996.

5. Berland S., Delattre O., Borzeix S., Catonné Y., and Lopez E. Nacre/bone interface changes in durable nacre endosseous implants in sheep. *Biomaterials*, 26(15):2767–2773, 2005. https://doi.org/10.1016/j.biomaterials.2004.07.019 PMID: 15585281

6. Espinosa H. D., Juster A. L., Latourte F. J., Loh O. Y., Gregoire D., and Zavattieri P. D. Tablet-level origin of toughening in abalone shells and translation to synthetic composite materials. *Nature communications*, 2:173, 2011. https://doi.org/10.1038/ncomms1172 PMID: 21285951

7. Jackson D. J., Wörheide G., and Degnan B. M. Dynamic expression of ancient and novel molluscan shell genes during ecological transitions. *BMC evolutionary biology*, 7(1):160, 2007. https://doi.org/10.1186/1471-2148-7-160 PMID: 17845714

8. Zlotnikov I. and Schoeppler V. Thermodynamic aspects of molluscan shell ultrastructural morphogenesis. *Advanced Functional Materials*, 27(28):1700506, 2017. https://doi.org/10.1002/adfm.201700506

9. Audino J. A., Marian J. E. A. R., Wanninger A., and Lopes S. G. B. C. Mantle margin morphogenesis in nodipecten nodosus (mollusca: Bivalvia): new insights into the development and the roles of bivalve palial folds. *BMC developmental biology*, 15(1):22, 2015. https://doi.org/10.1186/s12861-015-0074-9 PMID: 26017922

10. McDougall C., Green K., Jackson D. J., and Degnan B. M. Ultrastructure of the mantle of the gastropod halotiis asinina and mechanisms of shell regionalization. *Cells Tissues Organs*, 194(2-4):103–107, 2011. https://doi.org/10.1159/000324213 PMID: 21525717

11. Westermann B., Schmidtberg H., and Beuerlein K. Functional morphology of the mantle of nautilus pompilius (mollusca, cephalopoda). *Journal of morphology*, 264(3):277–285, 2005. https://doi.org/10.1002/jmor.10321 PMID: 15803486

12. Checa A. G., Salas C., Harper E. M., and de Dios Bueno-Perez J. Early stage biomineralization in the periostracum of the 'living fossil' bivalve neotrigonia. *PLoS One*, 9(2):e90033, 2014. https://doi.org/10.1371/journal.pone.0090033

13. Meinhardt H. and Meinhardt H. *Models of biological pattern formation*, volume 6. Academic Press London, 1982.

14. Boettiger A., Ermentrout B., and Oster G. The neural origins of shell structure and pattern in aquatic mollusks. *Proceedings of the National Academy of Sciences*, pages pnas–0810311106, 2009. https://doi.org/10.1073/pnas.0810311106

15. Checa A. A model for the morphogenesis of ribs in ammonites inferred from associated microsculptures. *Palaeontology*, 37(4):863–888, 1994.

16. Checa A. G. and Crampton J. S. Mechanics of sculpture formation in magadiceramus? rangatira rangatira (inoceramidae, bivalvia) from the upper cretaceous of new Zealand. *Lethaia*, 35(4):279–290, 2002.

17. Checa A. G. and Jiménez-Jiménez A. P. Rib fabrication in ostreoidea and plicatuloidea (bivalvia, pteriomorphia) and its evolutionary significance. *Zoomorphology*, 122(3):145–159, 2003. https://doi.org/10.1007/s00435-003-0080-5

18. Hammer O. A theory for the formation of commarginal ribs in mollusc shells by regulative oscillation. *Journal of Molluscan Studies*, 66(3):383–392, 2000. https://doi.org/10.1093/mollus/66.3.383

19. Morita R. Finite element analysis of a double membrane tube (dms-tube) and its implication for gastropod shell morphology. *Journal of Morphology*, 207(1):81–92, 1991. https://doi.org/10.1002/jmor.1052070110 PMID: 29865519

20. Ambrosi D., Ateshian G. A., Arruda E. M., Cowin S. C., Dumais J., Goriely A., Holzapfel G. A., Humphrey J. D., Kemkemer R., Kuhl E., Olberding J. E., Taber L. A., and Garikipati K. Perspectives on biological growth and remodeling. *J. Mech. Phys. Solids*, 59:863–883, 2011. https://doi.org/10.1016/j.jmps.2010.12.011 PMID: 21532929

21. Goriely A. *The mathematics and mechanics of biological growth*, volume 45. Springer, 2017.

22. Chirat R., Moulton D. E., and Goriely A. Mechanical basis of morphogenesis and convergent evolution of spiny seashells. *Proceedings of the National Academy of Sciences*, 110(15):6015–6020, 2013. https://doi.org/10.1073/pnas.1220443110

23. Erlich A., Moulton D. E., Goriely A., and Chirat R. Morphomechanics and developmental constraints in the evolution of ammonites shell form. *Journal of Experimental Zoology Part B: Molecular and Developmental Evolution*, 326(7):437–450, 2016. https://doi.org/10.1002/jez.b.22716
24. Erlich A., Howell R., Goriely A., Chirat R., and Moulton D. E. Mechanical feedback in seashell growth and form. *The ANZIAM Journal*, pages 1–26, 2018.

25. Moulton D. E., Goriely A., and Chirat R. Mechanical growth and morphogenesis of seashells. *Journal of theoretical biology*, 311:69–79, 2012. https://doi.org/10.1016/j.jtbi.2012.07.009 PMID: 22820493

26. Moulton D.E., Goriely A., and Chirat R. The morpho-mechanical basis of ammonite form. *Journal of theoretical biology*, 364:220–230, 2015. https://doi.org/10.1016/j.jtbi.2014.09.021 PMID: 25258003

27. Thompson D'A. W. *On Growth and Form*. Cambridge University Press, 1942.

28. T. E. Saunders and P. W. Ingham. Mod special issue celebrating 100 years since "on growth and form" by d'arcy wentworth thompson., 2017.

29. Rodriguez E. K., Hoger A., and McCulloch A. Stress-dependent finite growth in soft elastic tissue. *J. Bio-mech.*, 27:455–467, 1994. https://doi.org/10.1016/0021-9290(94)90021-3 PMID: 8188726

30. Ben Amar M. and Goriely A. Growth and instability in elastic tissues. *J. Mech. Phys. Solids*, 53:2284–2319, 2005. https://doi.org/10.1016/j.jmps.2005.04.008

31. Garikipati K. The kinematics of biological growth. *Appl. Mech. Rev.*, 62:030801, 2009. https://doi.org/10.1115/1.3090829

32. Goriely A. and Ben Amar M. Differential growth and instability in elastic shells. *Phys. Rev. Lett.*, 94 (19):198103, 2005. https://doi.org/10.1103/PhysRevLett.94.198103 PMID: 16090217

33. Garikipati K., Arruda E. M., Grosh K., Narayanan H., and Calve S. A continuum treatment of growth in biological tissue: the coupling of mass transport and mechanics. *J. Mech. Phys. Solids*, 52(7):1595–1625, 2004. https://doi.org/10.1016/j.jmps.2004.01.004

34. Alzetta G., Arndt D., Bangerth W., Boddu V., Brands B., Davydov D., Gassmoeller R., Heister T., Heltai L., Kormann K., Kronbichler M., Maier M., Pelteret J.-P., Turcksin B., and Wells D. The *deal.II* library, version 9.0. *Journal of Numerical Mathematics*, 26(4):173–183, 2018. https://doi.org/10.1515/jnma-2018-0054

35. S. Rudraraju and K. Garikipati. shellMorph. https://github.com/mechanoChem/shellMorph, 2019.

36. Moulton D. E, Lessinnes T., and Goriely A. Morphoelastic rods Part 1: A single growing elastic rod. *Journal of the Mechanics and Physics of Solids*, 61(2):398–427, 2012. https://doi.org/10.1016/j.jmps.2012.09.017

37. Didier Descouens. Hexaplex cichoreum. https://commons.wikimedia.org/wiki/File:Hexaplexcichoreum_Profils.jpg, 2009 (accessed December 30, 2018).

38. Meinhardt H. The algorithmic beauty of sea shells. *Cahiers de biologie marine*, 36(3):239–240, 1995.

39. Goriely A., Geers M. G. D., Holzapfel G. A., Jayamohan J., Jérusalem A., Sivaloganathan S., Squier W., van Dommelen J. A. W., Waters S. L., and Kuhl E. Mechanics of the brain: perspectives, challenges, and opportunities. *Biomech. Model. Mechanobiol.*, 14:931, 2015. https://doi.org/10.1007/s10237-015-0662-4 PMID: 25716305

40. Narayanan H., Verner S. N., Mills K. L., Kemkemer R., and Garikipati K. In silico estimates of the free energy rates in growing tumor spheroids. *Journal of Physics: Condensed Matter*, 22(19):194122, 2010. https://doi.org/10.1088/0953-8984/22/19/194122 PMID: 21386444

41. Rudraraju S., Mills K. L., Kemkemer R., and Garikipati K. Multiphysics modeling of reactions, mass transport and mechanics of tumor growth. In *Computer Models in Biomechanics*, pages 293–303. Springer, 2013.

42. Mills K. L., Kemkemer R., Rudraraju S., and Garikipati K. Elastic free energy drives the shape of prevascular solid tumors. *PLOS ONE*, 9(7):1–7, 07 2014. https://doi.org/10.1371/journal.pone.0103245

43. Bayly P. V., Taber L. A., and Kroenke C. D. Mechanical forces in cerebral cortical folding: a review of measurements and models. *J. Mech. Behav. Biomed. Mater.*, 29:568–581, 2014. https://doi.org/10.1016/j.jmbbm.2013.02.018 PMID: 23566768