# MECHANICAL FEEDBACK IN SEASHELL GROWTH AND FORM

Alexander Erlich, R. Howell, A. Goriely, R. Chirat, D. Moulton

► To cite this version:

Alexander Erlich, R. Howell, A. Goriely, R. Chirat, D. Moulton. MECHANICAL FEEDBACK IN SEASHELL GROWTH AND FORM. The ANZIAM Journal, 2018, 59 (4), pp.581-606. 10.1017/S1446181118000019. hal-04800368

HAL Id: hal-04800368

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

Submitted on 8 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mechanical feedback in seashell growth and form

A. Erlich<sup>a,f</sup>, R. Howell<sup>b,c,f</sup>, A. Goriely<sup>d</sup>, R. Chirat<sup>e</sup>, D.E. Moulton<sup>d,*</sup>

<sup>a</sup>School of Mathematics, University of Manchester

<sup>b</sup>Randall Division of Cell and Molecular Biophysics, King's College London

<sup>c</sup>Francis Crick Institute, London.

<sup>d</sup>Mathematical Institute, University of Oxford, Oxford, UK

<sup>e</sup>Université Lyon 1, CNRS UMR 5276 LGL-TPE, France

<sup>f</sup>Co-first author

## Abstract

Mollusc seashells grow through the local deposition and calcification of material at the shell opening by a soft and thin organ called the mantle. Through this process, a huge variety of shell structures are formed. Previous models have shown that these structural patterns can largely be understood by examining the mechanical interaction between the deformable mantle and the rigid shell aperture to which it adheres. In this paper we extend this modelling framework in two distinct directions. For one, we incorporate a mechanical feedback in the growth of the mollusc. Second, we develop an initial framework to couple the two primary and orthogonal modes of pattern formation in shells, which are termed antimarginal and commarginal ornamentation. In both cases we examine the change in shell morphology that occurs due to the different mechanical influences and evaluate the hypotheses in light of the fossil record.

Keywords: morphogenesis, morphomechanics, mathematical model, mollusc

## 1. Introduction

Seashells provide an intriguing case study for the role of mechanics in pattern formation. Shells are formed incrementally at the shell opening, termed the aperture, by a part of the mollusc's anatomy called the mantle. The mantle is a soft and thin organ that adheres to the shell opening and secretes a new layer of shell material, which then calcifies and hardens [1]. Previous models [2, 3, 4, 5] have shown that various structural patterns found on mollusc shells can be largely understood as a result of the mechanical interactions between the soft mantle and the rigid shell to which it adheres and which itself creates. Structural patterns, or ornamentations, in shells are usually classified as being commarginal or antimarginal. Essentially, the difference is in the direction of

<footnote>*Corresponding author</footnote>

<footnote>Email address: moulton@maths.ox.ac.uk (D.E. Moulton)</footnote>
<figure><img src="image_2.png" /></figure>

Figure 1: Examples of a) commarginal ornamentation, where the pattern forms parallel to the growth lines, and b) antimarginal ornamentation, with the pattern forming orthogonal to the growth lines. The shells are a) *Promicroceras planicosta* and b) *Tridacna squamosa*. The direction of growth lines are highlighted in red. Scale bars are 10 mm.

the pattern: in the former, the pattern forms parallel to the growth lines (i.e. in the direction of growth), such as in the ribs of ammonite shells (Fig 1a); for antimarginal ornamentation, the pattern forms orthogonal to the growth lines, as seen for instance in spines and ridges (Fig 1b).

The basic premise of these previous models of the mechanical shell forming process is to view the mantle margin as a growing elastic material with small ratio of width and depth to length, such that the mantle is subject to forces due to its constrained geometry and its adherence to the shell aperture. At each increment of growth, mechanical equilibrium (the system evolves quasistatically) determines the shape of the mantle, and this shape is then ‘sealed’ in the next iteration of the shell through calcification. The system is similar to an elastic rod on a foundation [6], in which the rod continually increases in length through development, and the foundation evolves following the equilibrium shape of the rod.

These models are able to generate a range of realistic shell morphologies as the geometric and growth parameters are varied and provide a simple and natural physical basis for many of the patterns seen in mollusc shells, both in extant species and in the fossil record. However, there are two particular simplifying assumptions present in these models that we address here.

1. No mechanical feedback in the growth process was considered. That is, the growth of the mollusc was a fixed input parameter to the system, completely decoupled from the mechanical forces and stress experienced through development.

2. Commarginal and antimarginal patterning are completely separated. That is, each model has restricted the class of deformations of the model so that only one type of pattern is permitted.
These assumptions enabled for analytical progress and investigation of an isolated phenomenon. In moving towards a more realistic description of shell growth, in this paper we extend the basic modelling framework with these assumptions in mind.

We first consider the extension of the modelling framework to include mechanical feedback on the growth rate of the mollusc. Over the past decades, mechanobiology has emerged as a discipline dedicated to the study of the effects of mechanical forces as epigenetic regulators of developmental processes. A growing body of studies have shown that mechanical forces generated within cells and tissues provide epigenetic regulatory signals that are as crucial as chemical factors for control of developmental processes. Thus, developmental control is now viewed as a mechanochemical process involving reciprocal interactions between both mechanical and chemical cues, mechanotransduction describing the processes that transduce mechanical stimuli into biochemical signals. These mechanical cues alter cellular signalling and switch cells between different fates (e.g., growth, differentiation, motility, apoptosis) by modulating cell shape and activating specific mechanotransduction pathways (e.g. [7]). More specifically, it has been shown that the geometry of the entire tissue and its resulting distribution of mechanical stress feeds back to regulate spatial patterns of cell proliferation [8], the positions of increased growth rate corresponding to sites under mechanical tension, the tissue geometry appearing not only as a consequence but also an active regulator of tissue growth [9]. An interesting aspect of mollusc shell growth is that the mechanical stresses in the secreting mantle emerge from its interactions with the hard shell which the mantle itself is secreting. A basic question is then: does the mechanical stress generated within the mantle feed back to alter the mollusc's growth rate during development?

We examine this question in Section 3. We again take as a starting point the model for ammonite rib formation derived in [4], and incorporate a mechanical feedback law into the previous modelling framework: the growth rate is reduced when in compression and increased when in tension. Mathematically, the full system becomes intractable analytically, however in the limit of small feedback we find the presence of a multiscale pattern, for which we can generate explicit asymptotic solutions using the method of multiple scales.

Regarding point 2, as stated antimarginal and commarginal ornamentation represent orthogonal directions of pattern formation. In practice, both patterns may occur, depending on the specific geometry, growth rate and mechanical properties, and the pattern that is ultimately seen may be dominated by one or the other or may be a combination of the two. However, fully incorporating both patterns involves a much broader potential morphology and is thus computationally quite challenging. As a first step towards this, here we consider a one-directional coupling of pattern selection by considering the possibility of antimarginal ornamentation forming during a growth process which produces a commarginal pattern. In particular, we combine the mechanical ideas of [4] and [3] – the intuitive idea underlying the formation of ammonite ribs [4] is that the closed ring of mantle tissue continuously oscillates between states of compression and tension, pushing the shell outward while in compression and pulling
the shell inward in tension, thereby producing an oscillatory shell opening. On the other hand, the physical principle we have suggested for the formation of spines [3] is that a growth spurt (rapid expansion of the mantle) creates enough compression in the mantle tissue to instigate a mechanical instability. The buckled shape of the mantle tissue is calcified and thus forms a template, effectively amplifying the pattern as growth continues.

In some sense, these patterns are distinct: the commarginal ribs of the ammonite occur as a continuous oscillation, while spine formation is a process occurring at discrete growth spurts. Nevertheless, the physical principle of the buckling that initiates the spines can certainly be present even without a growth spurt, and it is this aspect that we explore here. Our objective in Section 4 is to investigate the coexistence of the two patterns within the simplified geometry of a growing ring. We first formulate a model for the circumferential buckling of a growing elastic ring attached to a rigid foundation. This model extends previous work on buckling of morphoelastic rods [6], though as we demonstrate, the extension to a circular geometry is non-trivial. We then apply the buckling criterion to the case of a growing ammonite exhibiting commarginal ornamentation, in which case buckling in the circumferential direction during the compression phase of the oscillation creates an antimarginal patterning.

While shell growth forms an intriguing system to study for mechanical pattern formation, it is inherently difficult to provide experimental validation. For instance, in the case of commarginal ribs, ammonites have been extinct for 65 million years, and so one can only examine the (extremely extensive) fossil record. For each of the model extensions, our approach to evaluating model hypotheses is to compare the morphologies predicted from the models with the fossil record to determine whether there is sufficient evidence for the proposed effects. Moreover, our theoretical approach can potentially motivate new experimental studies on shell morphogenesis in living molluscs, notably in gastropods that may display commarginal or antimarginal ornamentations.

## 2. Model setup

The model considered here follows the same principles of our previous works [4, 5]: the shell producing mantle edge is treated as a growing elastic ring, with mechanical forces imposed by the interaction with the current shell edge and not-yet-calcified generative zone. The model is shown schematically in Fig. 2. We work within a fixed circular geometry, such that the shell radius is $r(z)$, where $z$ denotes the axial direction, i.e. the growth axis, and the shell makes angle $\phi(z)$ with the $z$ axis. In this formulation, the shell is akin to a growing cone, with shell coiling subsequently added as a secondary effect.

The stress free radius of the mantle edge (i.e. the radius that would be observed if the mantle were dissected from the shell) is denoted by $R(z)$. Intuitively, the basic idea is that a discrepancy between this radius and the radius of the calcified shell $r(z)$ will induce stress in the mantle and alter its size and thus the form of the developing shell. The region connecting the mantle edge to the already calcified portion of the shell is termed the generative zone. This
<figure><img src="image_3.png" /><figcaption>Figure 2: Setup for modelling commarginal ribbing. The shell edge and mantle edge (in its stress-free state) are treated as circles with radius $r(z)$ and $R(z)$, respectively, where $z$ denotes distance along the growth direction. The angle the shell opening makes with the growth axis is denoted $\phi(z)$.</figcaption></figure>

region is pliable, as it is yet to be calcified, but also forms a continuum with the rigid calcified portion of the shell, and thus acts like a torque spring, resisting a change in angle. That is, in the absence of any elastic forces in the mantle, the generative zone would maintain the same angle as the previously calcified portion of the shell edge to which it is attached, and work must be done to change the orientation.

The actual orientation of the generative zone and radius of the newly secreted shell are thus determined by the balance of the elastic force in the mantle and a generative zone force, which takes the form

$$f_{GZ} = k_{GZ} \frac{d\phi}{dz}, \qquad (1)$$

where $k_{GZ}$ is a stiffness parameter characterising the strength of the resistance to change orientation. In the circular geometry, all forces act in the radial direction and a radial force balance gives [6]

$$\frac{n}{R} + f_{GZ} = 0. \qquad (2)$$

where

$$n = E_s \left( \frac{r(z)}{R(z)} - 1 \right) \qquad (3)$$
is the resultant tangential force in the mantle, with $E_s$ the axial stiffness (Young's modulus times cross sectional area). We thus have the following differential equation for the angle $\phi(z)$:

$$ \frac{d\phi}{dz} = \frac{k}{R(z)} \left( 1 - \frac{r(z)}{R(z)} \right), \qquad (4) $$

where $k = E_s/k_{GZ}$ characterises the relative stiffness of mantle to generative zone. We also have the geometrical relation

$$ \frac{dr}{dz} = \tan \phi(z). \qquad (5) $$

In previous work, the reference radius $R(z)$ was taken as a prescribed function. As stated, a primary goal in this work is to incorporate the notion of mechanical feedback on growth, and thus we allow the reference radius to vary based on the internal mantle stress felt by the mantle. Here we consider the simplest form of a linear feedback:

$$ \frac{dR}{dz} = g + Mn. \qquad (6) $$

Here $g$ describes the basal expansion rate of the mantle. We take $g$ to be constant, denoting a base linear expansion of the mantle tissue through development, an assumption consistent with fossil measurements [4]. The second term in (6) captures mechanical feedback based on the internal stress $n$: if the mantle is in compression ($n < 0$), it reduces its expansion rate, thus relieving the degree of stress, while tension ($n > 0$) leads to an increased expansion rate. The constant parameter $M$ determines the strength of this feedback mechanism. Growth towards an optimal homeostatic stress $n^*$ may be considered by replacing $n$ with $n - n^*$ in (6). Also, while we restrict attention to linear feedback here, a non-linear feedback law based on Rayleigh dissipation was explored in the thesis work [10].

The system is then completed with initial conditions

$$ R(0) = R_0, \quad r(0) = R_0 + \delta R, \quad \phi(0) = \phi_0. \qquad (7) $$

### 2.1. Base case

The basic system described by equations (4)-(7) above (in the case of no mechanical feedback, i.e. $M = 0$) has been studied originally in [4]. Before proceeding, we briefly summarize the results of the base case.

In the case of a small initial displacement ($\epsilon = \delta R/R_0 \ll 1$), an asymptotic expansion may be performed [4], yielding the solution

$$ r(z) = R_0 + gz + \delta R A \left( 1 + \frac{g}{R_0} z \right)^{\frac{1}{2}} \cos \left[ \omega \log \left( 1 + \frac{g}{R_0} z \right) + \varphi \right] \qquad (8) $$

where

$$ A = \frac{\sqrt{1 + 4\omega^2}}{2\omega}, \quad \omega = \frac{\sqrt{4k(1 + g^2) - g^2}}{2g}, \quad \varphi = \arctan \left( \frac{1}{2\omega} \right). \qquad (9) $$
The first two terms in $r(z)$ defined in (8) coincide with the linearly expanding reference radius $R(z)$. The general behaviour then is that the shell oscillates about the reference radius, with increasing wavelength through development and a ribbing amplitude that grows like $\sqrt{z}$.

### 3. Mechanical feedback

In this section we extend the base case to consider the influence of the mechanical feedback term on the reference state $R(z)$. To simplify the analysis we choose $n^* = 0$ for the homeostatic stress. Equations (4)-(6) can then be written

$$ 
\begin{aligned}
r'(z) &= \tan \phi(z) \\
\phi'(z) &= \frac{k}{R(z)} \left( 1 - \frac{r(z)}{R(z)} \right) \\
R'(z) &= g + m \left( \frac{r(z)}{R(z)} - 1 \right)
\end{aligned} 
\quad (10) 
$$

where $m = ME_s$ and prime denotes differentiation.

#### 3.1. Displacement and non-dimensionalisation

To make analytical progress, it is convenient to work in terms of the displacement

$$ u(z) = r(z) - R(z). \quad (11) $$

The system (10) transforms to

$$ u'(z) = \tan \phi - g - m \frac{u(z)}{R(z)}, \quad \phi'(z) = -k \frac{u(z)}{R(z)^2}, \quad R'(z) = g + m \frac{u(z)}{R(z)}, \quad (12) $$

along with initial conditions

$$ u(0) = \delta R, \quad R(0) = R_0, \quad \phi(0) = \phi_0. $$

Next, we nondimensionalise:

$$ \tilde{z} = \ln \left( 1 + \frac{g}{R_0} z \right), \quad \tilde{u} = \frac{u}{\delta R}, \quad \tilde{R} = \frac{R}{R_0}. \quad (13) $$

The benefit of the definition of $\tilde{z}$ is it transforms the problem of integrating $u$ to a second order ODE with constant coefficients, as will become apparent below. In the case of a fixed linear expansion the ribbing pattern (8) oscillates with $\cos(\omega \log(1 + gR_0^{-1}z))$. With (13)$_1$, the oscillation simplifies to $\cos(\omega \tilde{z})$. Note also that the relationship between $z$ and $\tilde{z}$ is one-to-one, with $z = R_0 g^{-1} (e^{\tilde{z}} - 1)$ the unique inverse.

Inserting the above scalings into (12) gives
$$ 
\begin{aligned}
\epsilon g e^{-\tilde{z}} \tilde{u}'(\tilde{z}) &= \tan \phi(\tilde{z}) - g \left(1 + \epsilon \nu \frac{\tilde{u}}{\tilde{R}}\right), \\
g e^{-\tilde{z}} \phi'(\tilde{z}) &= -k \epsilon \frac{\tilde{u}}{\tilde{R}^2}, \\
\tilde{e}^{-\tilde{z}} R'(\tilde{z}) &= 1 + \epsilon \nu \frac{\tilde{u}}{\tilde{R}}.
\end{aligned} 
\quad (14) $$

The initial conditions scale to

$$ \tilde{u}(0) = 1, \quad \tilde{R}(0) = 1, \quad \phi(0) = \phi_0. \quad (15) $$

Here we have defined the dimensionless parameters

$$ \epsilon = \frac{\delta R}{R_0}, \quad \nu = \frac{m}{g}. $$

The parameter $\epsilon$ characterises the size of the initial displacement relative to radius, and is typically small, while $\nu$ relates the strength of the mechanical feedback compared to basal growth. To make progress, we now make the assumption that $\epsilon \ll 1$, and seek a solution as an asymptotic expansion of the form

$$ 
\begin{aligned}
\tilde{u} &\sim \tilde{u}_0 + \epsilon \tilde{u}_1 + \dots, \\
\tilde{R} &\sim \tilde{R}_0 + \epsilon \tilde{R}_1 + \dots, \\
\tilde{\phi} &\sim \tilde{\phi}_0 + \epsilon \tilde{\phi}_1 + \dots
\end{aligned} 
\quad (16) $$

At $O(1)$ we obtain the solution

$$ \tilde{\phi}_0 = \arctan g, \quad \tilde{R}_0 = e^{\tilde{z}}. $$

At $O(\epsilon)$ we have the system

$$ 
\begin{aligned}
g \tilde{u}'_0 + g \nu \tilde{u}_0 &= e^{\tilde{z}} (1 + g^2) \phi_1 \\
g \phi'_1 &= -k e^{-\tilde{z}} \tilde{u}_0 \\
\tilde{R}'_1 &= \nu u_0.
\end{aligned} 
\quad (17) $$

The first two equations can be written as a single constant-coefficient second order equation for $\tilde{u}_0$, whose general solution is

$$ \tilde{u}_0(\tilde{z}) = \exp\left(\frac{1-\nu}{2}\tilde{z}\right) (K_1 \cos(\omega_\nu \tilde{z}) + K_2 \sin(\omega_\nu \tilde{z})) \quad (18) $$

where

$$ \omega_\nu = \frac{\sqrt{4k(1+g^2) - g^2(1+\nu)^2}}{2g}. \quad (19) $$

Taking as initial condition

$$ \phi(0) = \phi_0 = \arctan g \Rightarrow \tilde{u}'_0(0) = -\nu, $$
we can determine $K_1 = 1$, $K_2 = -\frac{\nu+1}{2\omega_\nu}$. In terms of the original variables (13), and replacing the displacement (18) by the current radius $r = u + R$, we obtain

$$r(z) = R_0 + gz + \delta R A_\nu \left(1 + \frac{g}{R_0} z\right)^{\frac{1-\nu}{2}} \cos \left[\omega_\nu \log \left(1 + \frac{g}{R_0} z\right) + \varphi_\nu\right] \quad (20)$$

where

$$A_\nu = \frac{\sqrt{(\nu+1)^2 + 4\omega_\nu^2}}{2\omega_\nu}, \quad \varphi_\nu = \arctan \left(\frac{\nu+1}{2\omega_\nu}\right). \quad (21)$$

Setting $\nu = 0$ in (20) reproduces the solution without feedback, (8), i.e. $A_\nu$, $\varphi_\nu$, and $\omega_\nu$ tend to $A$, $\varphi$, and $\omega$, respectively, as $\nu \to 0$. In particular, noting that

$$\exp(\tilde{z}/2) \sim z^{1/2},$$

we see the ribbing amplitude increases as $\sqrt{z}$ when $\nu = 0$. However, for any $\nu > 0$, this amplitude growth will be decreased. Moreover, for $\nu > 1$, the amplitude of the oscillations switches from growth to decay. In terms of dimensional variables, for a given expansion rate $g$, since $\nu = ME_s/g$, the linearised model predicts a decaying ribbing amplitude if the feedback parameter $M$ satisfies

$$M > gE_s^{-1}.$$

Figure 3 shows the oscillation pattern of $u(z)$ according to (20). For the three chosen values of $\nu$, the amplitude of the oscillations is growing, constant and decaying, respectively. We also see a good match between the numerical solution of (12) and the asymptotic result (20).

## 3.2. Morphological trends in fossil record

A primary observation from ammonites' shells is that rapidly expanding shells tend to become smoother during development while the slowly expanding shells tend to remain more strongly ornamented at later stages of the development. This trend is generally supported by the basic morphomechanical model without feedback as described in Section 2 and analysed in [4, 5]. As shown above, the primary effect of the mechanical feedback is to change the amplitude of the ribbing pattern, in particular in the regime of moderate feedback a significant decay in amplitude is predicted.

It is difficult to conclude whether such effects are present in growing shells, in particular in ammonite shells, for which we can only examine the fossil record. Some clues that suggest the possible presence of a feedback mechanism can be observed in the species *Sonninia Propinquans*. The three shells shown in the top row of Figure 4 all show similar geometric characteristics, i.e. the rate of shell expansion and coiling are similar. Moreover, each of these shells is strongly ribbed in the juvenile stages, and each shell becomes smoother by the adult stage, but to a varying degree. In terms of the classic Raup morphospace [11], these shells have geometrical properties (relatively high expansion rate and eccentricity of the cross section) such that the mechanical model would
<figure><img src="image_4.png" /><figcaption>Figure 3: Oscillation pattern of $u(z)$ according to (20) for increasing values of $\nu$ (0.1, 1 and 5). As the mechanical feedback is increased, the amplitude of the oscillations switches from growing to constant to decaying. Also plotted is the numerical solution of the full system (12) (blue solid curve). Other parameters: $g = 0.067$, $k = 10$, $R(0) = 10$mm, $\delta R = 0.05$mm (which implies $\epsilon = 0.005$).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Photographs of *Sonninia Propinquans* and modelling results. The coiling parameters for the surface II.a were determined by fitting curves to photograph I.a. The other simulated shells, II.b and II.c, differ only in the mechanical feedback parameter *v*. Other parameters: *g*<sub>a</sub> = 0.0588, *τ* = 1.83, *k* = 850, *ε* = 0.03, *â*<sub>0</sub> = 16.32mm (see Appendix A).</figcaption></figure>

predict them to be weakly ornamented in the adult stage (see Fig. 5 of [4]). Following the protocol to extract the relevant parameters for the ribbing model (see Appendix B) and simulating a 3D shell in the case of zero feedback (ν = 0) produces a weakly ornamented shell qualitatively similar to the most ribbed specimen, I.a. Incorporating mechanical feedback, on the other hand, generates a much smoother adult shell (all other parameters held constant), and is thus able to account for the different degrees of ribbing observed<fnref n="1" />. While this is by no means proof, and other explanations may well account for the varying smoothness, such observations in the fossil record do at least lend credence to the idea. A similar mechanical feedback on growth rate could indeed account for the extensive intraspecific variability in the ribbing pattern in some species displaying low variation in shell geometry (e.g. [12]).

Morphospace

Allowing for the presence of mechanical feedback, the shell morphology of
ammonites predicted by the mechanical model above can largely be classified
by two distinct parameters: the expansion rate *g* and the degree of mechanical
feedback as described by the dimensionless parameter *v*. To these, following the
work in [5] we can add a third dimension: in [5], the circular model was adapted

<footnote n="1">¹To produce Fig 4 we have incorporated the mechanical feedback in the elliptical framework outlined in [5]. Details are provided in Appendix A. The same trend would be observed within a circular cross-sectional geometry, but an ellipse produces the best visual fit to the given shells.</footnote>
to an elliptical geometry, where we investigated the effect of the eccentricity of
the cross-section on the ribbing pattern. A primary observation from this anal-
ysis was that the greater the eccentricity, the smoother the shell appeared, i.e.
increased non-uniformity of cross-sectional curvature correlated with a damping
out of the ribbing pattern.

Incorporating eccentricity, the output of the ribbing model can thus be char-
acterised by two purely geometric parameters – the expansion rate and cross-
sectional eccentricity – and one mechanical feedback parameter. In Fig. 5, we
present a sample of the range of shell shapes that exist in this 3-dimensional
morphospace. Each parameter has a strong and distinct impact on the rib-
bing pattern produced<fnref n="2" />. Such a formulation enables evolutionary trends to be
considered in a new, mechanical perspective. Indeed, this perspective suggests
that combinatorial variations of these parameters modulating the oscillatory be-
haviour of the shell secreting system may underlie the exceptional morphological
diversity of ammonite shells.

4. Commarginal and antimarginal link

In this section we explore the possibility of antimarginal pattern formation
induced through the oscillations of a commarginal pattern. That is, within the
system of the rib producing mechanical oscillations (a commarginal ornamenta-
tion pattern), we ask whether transverse buckling could occur in the compression
phase of the oscillation, which would form an antimarginal pattern. For this,
we first derive a buckling criterion for a growing planar extensible elastic ring
attached to a foundation.

4.1. *Derivation of mechanical equations*

We model the generative region as an extensible ring with initial radius, $A$, attached to a fixed, circular, Winkler foundation (representing the calcified shell) also of radius $A$. We then consider this ring growing so that it has a stress free radius of $\gamma A$. As the ring grows, it will eventually buckle into a non-circular shape and we wish to determine the buckling threshold, $\gamma_b$, at which buckling first occurs. To this end we first adapt the computations of [6, 13] for use in a planar circular geometry.

4.2. *Geometry*

Following the general notion of the decomposition of the deformation gra-
dient in morphoelasticity [14] and applied to elastic rods as in the framework
of morphorods [6], there are 3 frames we must consider for this problem: an
initial, stress free circular frame of radius $A$ and parameterised by arc length,
$S_0$; a hypothetical grown but stress free frame which remains circular but with
radius $\gamma A$ and parameterised by arc length $S$; and the current, stressed state,

<footnote n="2">²The shell simulating process is outlined in detail in Appendix B.3, including a web interface for simulating shell form and producing 3D printed shells of output.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 5: 3-dimensional morphospace. The triplet of values of expansion rate, reference eccentricity, and mechanical feedback are given for each surface in the form $(g, \hat{e}, \nu)$. An increase in the expansion rate $g$ increases involution and decreases ornamentation; an increase in eccentricity $\hat{e}$ increases compression and decreases ornamentation; increased mechanical feedback $\nu$ decreases ornamentation. Details regarding the eccentricity and the elliptical model are found in [5] and briefly summarised in Appendix A. Other parameters: $\tau = 1$, $k = 400$, $\epsilon = 0.09$, $a_0 = 16.32\text{mm}$ (see Appendix B).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Schematic diagram showing the three frames involved in the problem of a buckling ring. a) The ungrown (attachment) radius, b) the unstressed grown radius, c) the current, buckled state.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Schematic diagram showing the orientation of frames in the buckled state and indicating the measurement of the angles $\theta$ and $\phi$. The shell radius is drawn in red with the polar coordinate vectors drawn in black and the local coordinate vectors drawn in blue.</figcaption></figure>

parameterised by its arc length, $s$, as indicated in Fig 6. The objective is to determine the smallest value of $\gamma > 1$ for which the current state admits a non-circular, buckled solution, and to determine the corresponding mode number for this buckled state.

As we are working in an initially circular geometry, a polar coordinate system, {$e_r, e_\theta, e_z$} is an obvious choice, although note that the majority of the analysis requires only planar considerations. We also make use of a local, right handed, frame defined by the current (buckled) configuration, {$d_1, d_2, d_3$}, where $d_1$ is associated with the inward pointing normal, $d_2$ with $e_z$ and $d_3$ with the tangent. Then designating $\theta$ as the polar angle and defining $\phi$ to be the positive angle between $e_\theta$ and $d_3$ (see figure 7), we express the shape in the current configuration as

$$ \mathbf{r} = f(\theta)e_r. $$
The polar and local coordinate systems are related by

$$ \mathbf{d}_3 = \sin(\phi)\mathbf{e}_r + \cos(\phi)\mathbf{e}_\theta \qquad (22) $$

$$ \mathbf{d}_1 = -\cos(\phi)\mathbf{e}_r + \sin(\phi)\mathbf{e}_\theta. \qquad (23) $$

In addition to the growth factor, $\gamma$, we account for a stretch factor $\alpha$, as the ring is extensible; this is given by

$$ \alpha = \frac{ds}{dS} = \frac{ds}{d\theta} \frac{d\theta}{dS_0} \frac{dS_0}{dS}. $$

From the definition of the stress-free grown configuration

$$ \frac{dS_0}{dS} = \frac{1}{\gamma}. $$

As $s$ is an arclength parameterisation, properties of polar curves imply

$$ \frac{ds}{d\theta} = \sqrt{f^2 + f'^2}, $$

where $'$ refers to the derivative with respect to $\theta$. Hence, the stretch factor can be written

$$ \alpha = \frac{1}{\gamma} \frac{d\theta}{dS_0} \sqrt{f^2 + f'^2}. \qquad (24) $$

By the properties of the local basis

$$ \frac{d\mathbf{r}}{ds} = \mathbf{d}_3 \\ \Rightarrow \frac{d\mathbf{r}}{dS_0} = \alpha\gamma\mathbf{d}_3 = \frac{d\theta}{dS_0}\sqrt{f^2 + f'^2}(\sin(\phi)\mathbf{e}_r + \cos(\phi)\mathbf{e}_\theta). $$

Considering this derivative in the polar coordinate system gives

$$ \frac{d\mathbf{r}}{dS_0} = \frac{d\theta}{dS_0}(f'\mathbf{e}_r + f\mathbf{e}_\theta), $$

and enforcing equality here gives

$$ f' = \sqrt{f^2 + f'^2} \sin(\phi), \qquad (25) $$

$$ f = \sqrt{f^2 + f'^2} \cos(\phi). \qquad (26) $$

The strain curvature vector, $\mathbf{u}$, is defined by

$$ \frac{d\mathbf{d}_i}{dS} = \mathbf{u} \wedge \mathbf{d}_i. $$

Noting that

$$ \frac{d\mathbf{d}_3}{dS} = \frac{d\theta}{dS_0} \frac{dS_0}{dS} \frac{d\mathbf{d}_3}{d\theta} = \frac{1}{\gamma} \frac{d\theta}{dS_0} (1 - \phi') \mathbf{d}_1, $$

and similarly

$$ \frac{d\mathbf{d}_1}{dS} = -\frac{1}{\gamma} \frac{d\theta}{dS_0} (1 - \phi') \mathbf{d}_3, $$

we have

$$ \mathbf{u} = \frac{1}{\gamma} \frac{d\theta}{dS_0} (1 - \phi') \mathbf{d}_2. \qquad (27) $$
### 4.3. Mechanics

The buckling criterion is determined as a solution of mechanical equilibrium. First, we define the resultant force and moments in the polar coordinate system

$$ \mathbf{n}(\theta) = n_r \mathbf{e}_r + n_{\theta} \mathbf{e}_{\theta} $$

$$ \mathbf{m} = m \mathbf{e}_z. $$

Then, assuming the ring is at equilibrium, conservation of linear and angular momenta lead to

$$ \frac{d\mathbf{n}}{dS_0} + \mathbf{F} = 0 $$

$$ \frac{d\mathbf{m}}{dS_0} + \frac{d\mathbf{r}}{dS_0} \wedge \mathbf{n} = 0. $$

The external force term $\mathbf{F}$ derives from the attachment of the ring to the foundation. Assuming a Hookean relation, this takes the form

$$ \mathbf{F} = -\alpha\gamma k_f (f - A) \mathbf{e}_r, $$

where $k_f$ is a parameter describing the strength of attachment to the foundation. Note that the force is scaled by $\alpha\gamma$ to obtain an appropriate force per unit arc length $S_0$. The balance of linear momentum reduces, component-wise, to

$$ n'_{r} - n_{\theta} = \sqrt{f^2 + f'^2} k_f (f - R), \qquad (28) $$

$$ n'_{\theta} + n_{r} = 0. \qquad (29) $$

The balance of angular momentum, which has only one non-zero component in the $\mathbf{e}_z$ direction, reads

$$ m' + \sqrt{f^2 + f'^2} (\sin(\phi)n_{\theta} - \cos(\phi)n_{r}) = 0. \qquad (30) $$

To close the system we also specify the constitutive relations

$$ \mathbf{n} \cdot \mathbf{d}_3 = E_s (\alpha - 1) \qquad (31) $$

$$ m = E_b \mathbf{u} = E_b \frac{1}{\gamma} \frac{d\theta}{dS_0} (1 - \phi'), \qquad (32) $$

where $E_s$ and $E_b$ are the stretching and bending stiffness respectively.

### 4.4. Buckling criterion

The equations (25-32), along with the assumption of periodicity govern the shape of the buckled ring. To investigate the nature of solutions, and find the values of $\gamma$ that can induce buckling, we conduct a linear analysis. Noting that $\phi \equiv 0$ in the circular state, we make a small angle approximation:

$$ \phi = \epsilon\phi_1 + \dots, \qquad (33) $$
where $0 < \epsilon \ll 1$. As mentioned above, all variables except $\theta$ itself are defined in terms of $\theta$. To simplify the problem, we would like all terms to have the same dependent variable so we would like a characterisation of $d\theta/dS_0$ in terms of $\theta$. Noting that

$$ \theta(0) = 0; \quad \theta(2\pi A) = 2\pi, $$

and assuming $f$ is sufficiently smooth, a bijective mapping between $\theta$ and $S_0$ exists. Furthermore, for an approximately circular shape

$$ \frac{dS_0}{d\theta} = A + \epsilon\nu_1(\theta) + \dots, $$

where we have introduced an unknown function $\nu_1 = \mathcal{O}(1)$. Hence

$$ \frac{d\theta}{dS_0} = \frac{1}{\frac{dS_0}{d\theta}} = (A + \epsilon\nu_1 + \mathcal{O}(\epsilon^2))^{-1} \sim \frac{1}{A} - \frac{1}{A^2}\epsilon\nu_1 + \dots \quad (34) $$

In addition to these expansions we expand the remaining variables:

$$ \begin{aligned} f &= f_0 + \epsilon f_1 + \dots \\ n_{\theta} &= n_{\theta_0} + \epsilon n_{\theta_1} + \dots \\ n_r &= n_{r_0} + \epsilon n_{r_1} + \dots \\ m &= m_0 + \epsilon m_1 + \dots \end{aligned} $$

Considering $\mathcal{O}(1)$ terms of the equations, we find

$$ n_{\theta_0} = -k_f f_0 (f_0 - A), \quad (35) $$

$$ n_{r_0} = 0, \quad (36) $$

$$ m = \frac{E_b}{A\gamma}. \quad (37) $$

To determine $f_0$ we combine (35) with the constitutive relation (31), which at leading order reads

$$ n_{\theta_0} = E_s \left( \frac{f_0}{\gamma A} - 1 \right). \quad (38) $$

Eliminating $n_{\theta_0}$ from the system then requires $f_0$ to satisfy the following quadratic equation

$$ k_f f_0^2 + \left( \frac{E_s}{\gamma A} - A k_f \right) f_0 - E_s = 0. \quad (39) $$

Requiring a positive solution then uniquely determines $f_0$ as

$$ f_0 = \frac{-\left(\frac{E_s}{\gamma A} - A k_f\right) + \sqrt{\left(\frac{E_s}{\gamma A} - A k_f\right)^2 + 4 k_f E_s}}{2 k_f} \quad (40) $$

Note that in this linear stability analysis, the $\mathcal{O}(1)$ solution describes the circular shape prior to buckling, hence at leading order every quantity is constant. In
particular, $f_0$ gives the radius and $n_{\theta_0}$ the compressive stress for any growth $\gamma > 1$; when evaluated at the buckling threshold $\gamma_b$ (to be found), these provide expressions for the critical radius and stress.

At $\mathcal{O}(\epsilon)$ the equations (25-32) become

$$
\begin{align*}
f_1' &= f_0\phi_1 \\
n'_{r_1} - n_{\theta_1} &= k_f(2f_0 - A)f_1 \\
n_{\theta_1} + n_{r_1} &= 0 \\
m_1' + f_0n_{\theta_0}\phi_1 - f_0n_{r_1} &= 0 \\
n_{\theta_1} &= \frac{E_s}{\gamma} \left( \frac{f_1}{A} - \frac{f_0\nu_1}{A^2} \right) \\
m_1 &= \frac{E_b}{\gamma} \left( -\frac{\phi_1'}{A} - \frac{\nu_1}{A^2} \right)
\end{align*}
$$

Eliminating $\nu_1$ leaves a system of 5 linear ODEs in the variables ($f_1, \phi_1, n_{r_1}, n_{\theta_1}, m_1$). Defining

$$
\mathbf{X} = \begin{pmatrix} f_1 \\ \phi_1 \\ n_{r_1} \\ n_{\theta_1} \\ m_1 \end{pmatrix}, \qquad (41)
$$

this can be phrased in the form

$$
\mathbf{X}'(\theta) = M\mathbf{X}, \qquad (42)
$$

where

$$
M = \begin{pmatrix}
0 & f_0 & 0 & 0 & 0 \\
-1/f_0 & 0 & 0 & \gamma A / E_s f_0 & -\gamma A / E_b \\
k_f(2f_0 - A) & 0 & 0 & 1 & 0 \\
0 & 0 & -1 & 0 & 0 \\
0 & -f_0 n_{\theta_0} & f_0 & 0 & 0
\end{pmatrix} . \quad (43)
$$

Due to the required periodicity, we can look for solutions as a Fourier mode
**X** ~ e^iNθ, where N is an integer. Inserting into (42) gives

$$
(M - iNI)\mathbf{X} = \mathbf{0}.
$$

For non-trivial solutions, we thus require that

$$
\det(M - iNI) = 0, \qquad (44)
$$

which will yield a condition on the values of $\gamma$ required to achieve buckling mode $N$. This will allow us to determine the buckling threshold $\gamma_b$, which is the smallest value of $\gamma$ for any $N$; this is the point at which physical buckling will be observed, and the corresponding $N$ gives the mode of buckling.
Computing the determinant gives the condition

$$ \frac{Ni}{E_b E_s} (E_b E_s - 2E_b E_s N^2 + E_b E_s N^4 - E_b A^2 \gamma k_f + 2E_b A f_0 \gamma k_f - E_s A f_0 \gamma n_{\theta_0} + 2E_s A f_0^3 \gamma k_f - E_s A^2 f_0^2 \gamma k_f + E_s N^2 A f_0 \gamma n_{\theta_0}) = 0. \quad (45) $$

The factor of $Ni$ at the front shows that the equation reduces to a quartic in $N$ with real valued coefficients. Noting that $f_0$ and $n_{\theta_0}$ are given in terms of $\gamma$ by (38) and (40), (45) is of the form

$$ F(\gamma, N; E_s, E_b, A, k_f) = 0. \quad (46) $$

The parameters $E_s$, $E_b$, and $k_f$ are fixed material parameters, and $A$ is the fixed reference radius. For an elastica with rectangular cross section of height $h$, width $w$, and Young's modulus $E$, we have

$$ \begin{aligned} E_s &= Ewh, \\ E_b &= \frac{1}{12} E w^3 h. \end{aligned} $$

Similar to [3, 4], we take as base values $E = 1$MPa, $w = h = 0.5$mm. The shell radius $A$ changes throughout the formation of the shell. As a base value, $A = 1$cm will be used as an order of magnitude approximation. The parameter $k_f$ characterises the stiffness of the attachment of the generative region to the hardened shell but in practice it is not possible to obtain accurate measurements. Thus, we have investigated a range of values of $k_f$ to see when realistic behaviour is observed.

For fixed values of $E_s$, $E_b$, $k_f$, and $A$, buckling is determined by finding the first root of (46) for $\gamma > 1$, and then finding the minimum value over all modes $N$. The non-linear dependence of $f_0$ on $\gamma$ means that the roots cannot be computed analytically; we have thus implemented a numerical root finding and utilised a simple algorithm for minimisation over modes.

**Effect of foundation stiffness.** In Fig 8 we plot the buckling mode (a) and critical growth (b) as functions of foundation stiffness $k_f$. We see that increasing $k_f$ increases $N^*$ monotonically, while the critical growth $\gamma_b$ does not change monotonically with $k_f$, but rather has a minimum value and diverges as $k_f \to 0$.

To understand these features, it is worthwhile to compare and contrast with the case of buckling of a straight growing elastica, i.e. a rod growing between two walls, as considered in [6, 13]. In both geometries, buckling results from a competition between compressive energy, bending energy, and spring energy in the foundation. For a straight rod with no foundation, prior to buckling growth does not induce any deformation, the rod remains flat with increasing compressive stress. Buckling then occurs because at some point it is more beneficial to add bending energy in order to reduce the compressive energy. For a long straight rod, this tradeoff occurs with very little growth, as such a rod can endure almost no compression (note that in classic Euler buckling, the critical
<figure><img src="image_9.png" /><figcaption>Figure 8: The buckling mode $N_b$ (a) and critical growth $\gamma_b$ (b) as functions of foundation stiffness $k_f$. Here $A = 1\text{cm}$ is fixed.</figcaption></figure>

force scales like the inverse square of length [15]). Adding a foundation effec-
tively stabilises the rod: larger compression may be endured before buckling,
and when the rod does buckle this occurs at a higher mode, with more bending
energy spent to balance the spring energy.

In the ring geometry, there is a fundamental difference: without any attach-
ment to the foundation, the growing ring does deform, expanding uniformly
while developing zero compressive stress, and hence never buckles. This can
be seen by setting $k_f = 0$ in (35), (39), which yields $f_0 = \gamma A$, $n_{\theta_0} = 0$. This
explains the divergence in $\gamma_b$ as $k_f \to 0$ in Fig 8 (b). However, as in the flat
rod case, the foundation force (eventually) has a stabilising effect, so that for
large enough $k_f$ the critical growth again increases with increasing $k_f$. The
result is a minimum value in $\gamma_b(k_f)$. This qualitative difference also underlies
the increased computational complexity in the ring geometry, since the grown
but prebuckled shape is a nontrivial function of $\gamma$ (via (40)).

*Effect of radius.* Buckling mode and critical growth are plotted as functions of reference radius $A$ in Fig 9. Here, the behaviour is the same as in the straight rod geometry: the buckling mode increases and critical growth decreases, both monotonically, with increasing radius $A$. Note that $\gamma_b$ asymptotes to a value greater than 1 as $A \to \infty$. This limit is equivalent to the buckling growth obtained for an infinite straight rod, as computed in [6]. For discussion below, the buckling threshold $\gamma_b$ decays approximately as $1/A^2$.

4.5. *Link to commarginal patterning*

Having established a buckling criterion, we now consider the possibility of
antimarginal patterning forming during compression phases of commarginal os-
cillations. Here, careful consideration must be made to relate the two frame-
works. The commarginal ribbing model involves the evolution of two radii: the
growing reference radius $R(z)$ and the actual shell radius, $r(z)$, which oscil-
lates around the reference radius. The antimarginal model similarly involves
<figure><img src="image_10.png" /><figcaption>Figure 9: The buckling mode $N_b$ (a) and critical growth $\gamma_b$ (b) as functions of reference radius $R$. Here $k_f = 1$ kPa is fixed.</figcaption></figure>

two radii: the radius of the foundation, $A$, and the radius of the unbuckled, grown ring, $\gamma A$. The hardened shell plays the role of the foundation, with the still-calcifying generative zone applying the elastic adhesion, while the reference radius may be interpreted as the grown stress-free state. Thus we identify

$$ \gamma A = R(z), \quad A = r(z), $$

from which we can define the growth

$$ \gamma = \frac{R(z)}{r(z)}. $$

This expression makes it clear that the valleys of the commarginal ribs (or the local minima in the oscillatory shape, for which $r < R$) correspond to the regions of compression with $\gamma > 1$ and hence the location where buckling can potentially occur.

In Fig 10 we plot the shell radius $r(z)$, following equation (8), oscillating about the reference radius $R(z)$, which simply expands linearly without mechanical feedback in this case. The shaded yellow region denotes the buckling region, i.e. for each value of the reference radius (blue curve), if the mantle radius dips into the yellow region then the compression is sufficient to initiate circumferential buckling. In Fig 10(a), with $k_f = 1000$Pa, the shell radius enters the buckling region in each compression phase, and hence the model predicts that antimarginal buckling would occur in each valley of the commarginal pattern. Since these regions do not correspond to the sort of growth spurt that produces large scale ridges and ribs, we would only expect to see the formation of wrinkles in these regions. Indeed, there is such evidence in the fossil record, as for example illustrated in Fig. 11 and in many non-closely-related ammonites genera [16]. Whether the compression is sufficient to trigger buckling, however, depends on the geometrical and material parameters. For instance, Fig 10(b) plots the same buckling region for a softer foundation ($k_f = 40$Pa), and at
<figure><img src="image_11.png" /><figcaption>Figure 10: The radius of a growing shell exhibiting commariginal patterning. The buckling region is indicated as the shaded yellow region, determined by calculating the buckling threshold $\gamma_b$ for the given shell radius (red curve) and reference radius (blue line) at each point. In (a), $k_f = 1000$, and the radius enters the buckling region in each compression phase, i.e. the valleys of the ribbing pattern. In (b) $k_f = 40$ and buckling is triggered only at larger $z$. In (c), $k_f = 5$, and the softer foundation requires a greater deviation in radius for buckling than is produced during the ribbing pattern.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Photograph showing the presence of wrinkles (indicated by the arrows) in valleys between commarginal ribs in Quenstedtoceras lamberti (Callovian, Russia). These wrinkles are formed during the shell secretion process and are subsequently recorded in the calcified shell.</figcaption></figure>

this value the buckling threshold is only reached at larger $z$; in Fig 10(c) the foundation is even softer and buckling is never triggered.

A final observation involves the potential for and scaling of the pattern with size. As the shell develops and the radius generally increases, competing effects come into play. On the one hand, larger values of $A = r(z)$ lead to a decrease in $\gamma_b$, as seen in Fig 9(b), so that buckling will occur for smaller values of $\gamma$ as the radius increases. However, since $\gamma = R(z)/r(z)$, the extremal values of $\gamma$ achieved during the compression phase depend on the relative amplitude of the commarginal oscillations (amplitude divided by radius), which tend to decrease through development [4]. We can gain some insight by considering the case of a linear expansion with no mechanical feedback. From Equation (8), the maximum value of $\gamma = R/r$ produced by the oscillations occurs when the cosine equals negative one in $r(z)$, at which point we have

$$ \gamma_{\max} = \frac{R_0 + gz}{R_0 + gz - \epsilon R_0^{1/2} A (R_0 + gz)^{1/2}} \sim 1 + \frac{\epsilon R_0^{1/2} A(g)}{R(z)^{1/2}} + O(\epsilon^2), \quad (47) $$

where $A = A(g)$ depends on the expansion rate $g$ through (9). Meanwhile, as seen in Fig 9, the buckling threshold has the form

$$ \gamma_b = \gamma_1 + \frac{b}{R(z)^\alpha} \quad (48) $$

where $\gamma_1 > 1$ and $b > 0$ depend on material parameters, and the decay rate$^3 \alpha \approx 2$. Comparing (47) and (48), we see that $\gamma_b \gg \gamma_{\max}$ as $R \to 0$, while $\gamma_b > \gamma_{\max}$

<footnote>$^3$We have computed $\alpha$ through a least squares fit for several choices of stiffness parameters and found $\alpha \approx 2$ in all cases. The variance is irrelevant in the present discussion, since the arguments so long as $\alpha > 1/2$.</footnote>
in the limit $R \to \infty$. This suggests that if antimarginal buckling is to occur, it will not be present at the juvenile stage, since the buckling threshold is too high for small radii; it may occur for intermediate radii, and then would ultimately disappear again should the radius become large enough. This will all depend on the specific details of any given shell, i.e. the growth and material parameters. Also, the thickness of the mantle will increase through development, which may lead to a nonlinearly change in the stiffness parameters $E_b$ and $E_s$ in the buckling model, which would impact the realised values of $\gamma$. It is thus unlikely that there is a generic rule that universally holds. Moreover, the empirical evidence to test such predictions are slim, mainly because specimens must be exceptionally well preserved to display such wrinkles. Nevertheless, the wrinkles in the species *Quenstedtoceras lamberti*, as displayed in Fig 11, are indeed only present in the later stages, and we do not know of any specimens showing these wrinkles only at early stages. It is also generally observed in ammonites that antimarginal spines appear only during latter stages of development, which is also observed in many gastropods. As noted in [5], in some ammonite genera (e.g. *Aspidoceras* and *Collignoniceras*), tubercles or clavi (oscillations) during juvenile stages may also transform into spines at a later stage of development. This seems to suggest that in these genera, the buckling threshold is more likely to be met at larger radii.

## 5. Conclusions

In this paper we have extended the modelling framework for commarginal and antimarginal patterning in seashells in two distinct directions. First, we have incorporated a mechanical feedback law into the growth of the mollusc, following the simple principle that compressive forces in the mantle could slow the growth of the mollusc, while tensile forces could increase the growth rate. By analysing the behaviour of a linear feedback law in a circular shell margin geometry, we found that mechanical feedback has the general effect of diminishing the amplitude of the oscillations. Intuitively, this is not surprising: with feedback the mantle is varying its growth rate in an attempt to perfectly coincide with the shell margin. The result is a smoother shell than would otherwise be observed. Mechanical feedback also provides another important dimension in generating a “mechanical morphospace”, from which the great diversity of trends in shell patterning might be understood in terms of a small number of physical, mechanical parameters directly underlying the growth process. Here, we have only briefly examined the fossil record, and have found some evidence to corroborate the presence of mechanical feedback. A more thorough examination, both of the fossil record and possible mechanotransduction pathways in extant species, combined with comparison with model prediction, should be undertaken before concluding whether or not such mechanical feedback is or was actually present in shell building molluscs.

In the second half of this paper, we have taken a first modelling step towards combining the two orthogonal patterns appearing in shells: commarginal and antimarginal. This required the development of an antimarginal buckling
criterion within a ringed geometry, and connecting this to the commarginal
oscillation framework. We have found that the compression phase of the oscil-
lation pattern can, for realistic parameters, instigate a circumferential buckling,
but only for sufficiently stiff adherence between shell and mantle and/or for suf-
ficiently high ribbing amplitude. This framework may explain the presence of
wrinkling appearing in the valleys of certain ammonite species. While we have
speculated about the possible connection to larger scale antimarginal patterns
such as ridges or spines, it is important to note that such patterns do not fall
within the assumptions underlying the model we have developed. In particular,
we have only devised a buckling threshold via a linear stability analysis – the
evolution of the pattern beyond the initial buckling would require solving the
full nonlinear equations. Moreover, we have taken the commarginal pattern as
being fixed, while in reality if an antimarginal wrinkling were to occur, then the
ring geometry of the commarginal model would cease to be valid. From these
simple observations it should be clear that a full coupling of these patterns is
a complex biomechanical problem, and likely requires a sophisticated computa-
tional framework to explore fully. With such a framework in place, along with
coupling different modes of pattern, mechanical feedback could also be incor-
porated. This has the interesting feature that in a post-buckled, non-circular
geometry, the stress distribution will be non-uniform and hence the subsequent
growth will also be non-uniform. Whether there are regimes in which feedback
then serves to amplify the pattern remains to be seen. Such questions form
useful directions to pursue, not just in terms of generating a more realistic un-
derstanding of seashell morphology and the evolutionary trends of molluscs,
but as a generic framework for mechanical pattern formation in growing tissues.

**Appendix A. Elliptic geometry**

The effect of the eccentricity of the elliptic geometry on the ribbing pattern
is explored in detail in [5]. For the present purpose, and to understand the
effect of elliptic cross-sections in Figures 4 and 5, we require three functions,
two of which are independent: the semi-minor axis $b(z)$, the semi-major axis
$a(z)$ and the eccentricity $e(z)$ of the aperture ellipse (i.e. a measure of how
much it deviates from being circular). In the context of elliptic geometry, we
denote reference quantities with an overhat, i.e. $\hat{b}, \hat{a}, \hat{e}$. The eccentricity is
defined as $e = \sqrt{1 - \frac{b^2}{a^2}}$. We suppose that the mantle grows isometrically,
that is, without changing its reference eccentricity ($e = \hat{e} = \text{constant}$). An
implication of the isometric growth assumption is that quantities related to the
major axis (such as its initial reference length $\hat{a}_0$, grown reference length $\hat{a}$,
initial displacement $\delta a$ and expansion rate $g_a$) are functionally related to their
minor axis counterparts ($\hat{b}_0, \hat{b}, \delta b$ and $g_b$, respectively) via the fixed reference
eccentricity $\hat{e}$. This is important to note because similarly to (4)-(7), a problem
in isometrically growing elliptic geometry can be stated by providing the relevant
quantities related to, say, the major axis in addition to the eccentricity value $\hat{e}$.
The 3D surface plots shown in Figures 4 and 5 are indeed generated from major
<figure><img src="image_13.png" /><figcaption>Figure B.12: The logarithmic spirals (B.1) representing the shell centreline $\rho_m$, the dorsal side $\rho_d$ and the ventral side $\rho_v$. In the case of perfect contact between whorls, the tightness of coiling parameter is $\tau = 1$. In the case of overlap, $\tau > 1$. Figure adapted from [4].</figcaption></figure>

axis parameters $g_a$, $\hat{a}_0$, $\varepsilon = \delta a / \hat{a}_0$ as well as the eccentricity value $\hat{e}$ as given in the captions. An explicit form of the surface is detailed in Appendix B.

## Appendix B. Coiling and 3D seashell surface

In this section we relate the ribbing oscillations to the coiled 3D shells presented in Figures 3, 4 and 5. We split this task into two parts: First, we introduce coiling via a logarithmic spiral shape that forms the basis for the shell surface (i.e. the shape in the absence of ribbing) and relate the coiling parameters to the expansion rate $g_a$ of the major axis of the ribbing model. Then, we give a parametric representation of the fully ribbed and coiled seashell surface that underlies the figures in the main text.

### Appendix B.1. Description of shell coiling

A logarithmic spiral is a polar curve with radius $\rho$ and angle $\phi$ by $\rho(\phi) = c_1 e^{c_2 \phi}$, characterized by two coiling parameters $c_1, c_2$. In Cartesian coordinates $\{x, y\}$, the curve takes the form $x(\phi) = \rho(\phi) \cos \phi$, $y(\phi) = \rho(\phi) \sin \phi$. The (unribbed) ammonites' shell surface is modeled by three logarithmic spirals: The centreline spiral $\rho_m$, the dorsal spiral $\rho_d$, and the ventral spiral $\rho_v$. In polar coordinates, they are

$$ \rho_m(\phi) = c_1 e^{c_2 \phi}, \quad \rho_d(\phi) = (c_1 - \hat{a}) e^{c_2 \phi} \quad \text{and} \quad \rho_v(\phi) = (c_1 + \hat{a}) e^{c_2 \phi}. \qquad (B.1) $$

Here $\hat{a} = \hat{a}(\phi)$ is the half-distance between the dorsal and ventral side. If the low curvature side of the elliptical shell margin is oriented along the lateral side, we identify $\hat{a}$ with the reference major axis.

To characterise to what degree the coiling shell overlaps with itself, we introduce a tightness-of-coiling parameter $\tau$, following the ideas of [4] (see Figure 7 of that reference). For coiling in perfect contact, in which the ventral side
coincides with the dorsal side of the previous whorl, we have $\tau = 1$. Overlap between each whorl and the previous whorls corresponds to a value of $\tau > 1$. For an illustration see Figure B.12.

From the definition of $\tau$, we have that $\rho_v(\phi) = \tau\rho_d(\phi + 2\pi)$ which implies

$$e^{2\pi c_2} = \frac{1}{\tau} \left( \frac{c_1 + \hat{a}_0}{c_1 - \hat{a}_0} \right). \qquad (B.2)$$

To relate the coiling to the expansion rate $g_a$, we must first compute the arclength of the centreline spiral $\rho_m$. The centreline is described by the curve $\mathbf{r}_m(\phi) = \rho_m(\phi)\mathbf{e}_r(\phi)$ where $\mathbf{e}_r$ is the radial direction which in a Cartesian basis $\{\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z\}$ is $\mathbf{e}_r(\phi) = \cos\phi \mathbf{e}_x + \sin\phi \mathbf{e}_y$. Its arclength is

$$z(\phi) = \int_{0}^{\phi} \left| \frac{d\mathbf{r}_{m}(\varphi)}{d\varphi} \right| d\varphi = \frac{\sqrt{1+c_{2}^{2}}}{c_{2}} [\rho_{m}(\phi) - c_{1}] \qquad (B.3)$$

The reference length of the major axis satisfies by

$$\rho_v(\phi) - \rho_m(\phi) = \hat{a}(z(\phi)), \qquad (B.4)$$

where $\hat{a}(z(\phi)) = \hat{a}_0 + g_a z(\phi)$ in the case of linear expansion. Evaluating (B.4) provides the following relationship between $g_a$ and the coiling parameters:

$$\hat{a}_0 c_2 = c_1 g_a \sqrt{1 + c_2^2}. \qquad (B.5)$$

In order to obtain the coiling parameters $c_1$, $c_2$ in terms of $g$ and $\tau$, the equations (B.2) and (B.5) must be solved simultaneously, which can be achieved numerically.

## Appendix B.2. Parametric form of shell surface

The coiling formulation above enables to find the coiling parameters $c_1$ and $c_2$ in terms of the parameters $g$ and $\tau$, while the morphomechanical ribbing model presented in Sections 2 and 3 provides $a(z)$ and $b(z)$. The full 3D seashell surface is then parameterised by

$$\mathbf{r}(\theta, \phi) = \mathbf{r}_m(\phi) + b(z(\phi)) \cos \theta \mathbf{e}_z + a(z(\phi)) \sin \theta \mathbf{e}_r \qquad (B.6)$$

where $\theta$ is the cross-section ellipse parameter and $\phi$ is the coiling parameter. This surface parameterisation was used to plot all 3D seashells in this paper.

## Appendix B.3. Web simulation

To better illustrate the results of the model, we have created an interactive web simulation in which the user can examine the above described parametric surface as a function of shell parameters. The user can change the values of the stiffness parameter $k$, the tightness of coiling parameter $\tau$, the expansion rate $g_a$ of the ventral side, and the constant eccentricity $\hat{e}$. The seashell surface updates instantaneously and gives the option of producing an STL file for 3D printing of the shell. The simulation can be run from any browser on any desktop operating system at airlich.de/shell.
[1] Simkiss K, Wilbur KM. Biomineralization. Academic Press, New York; 1989.

[2] Moulton DE, Goriely A, Chirat R. Mechanical growth and morphogenesis of seashells. Journal of theoretical biology. 2012;311:69-79.

[3] Chirat R, Moulton DE, Goriely A. Mechanical basis of morphogenesis and convergent evolution of spiny seashells. Proceedings of the National Academy of Sciences of the United States of America. 2013;110(15):6015-6020.

[4] Moulton DE, Goriely A, Chirat R. The morpho-mechanical basis of ammonite form. Journal of theoretical biology. 2015 Jan;364(C):220-230.

[5] Erlich A, Moulton DE, Goriely A, Chirat R. Morphomechanics and Developmental Constraints in the Evolution of Ammonites Shell Form. Journal of Experimental Zoology Part B: Molecular and Developmental Evolution. 2016 Dec;326(7):437-450.

[6] Moulton DE, Lessinnes T, Goriely A. Morphoelastic rods Part 1: A single growing elastic rod. Journal of the Mechanics and Physics of Solids. 2012;61(2):398-427.

[7] Mammoto T, Mammoto A, Ingber DE. Mechanobiology and developmental control. Annual review of cell and developmental biology. 2013:29:27-61.

[8] Hutson MS, Ma X. Mechanical aspects of developmental biology: perspectives On Growth and Form in the (post)-genomic age. Physical biology. 2008;5(1):015001.

[9] Nelson CM, Jean RP, Tan JL, Liu WF, Sniadecki NJ, Spector AA, et al. Emergent patterns of growth controlled by multicellular form and mechanics. Proceedings of the National Academy of Sciences of the United States of America. 2005;102(33):11594-11599.

[10] Erlich A. Growth Laws in Morphoelasticity; 2017.

[11] Raup DM. Geometric analysis of shell coiling: coiling in ammonoids. Journal of Paleontology. 1967;p. 43-65.

[12] Dagys A, Bucher H, Weitschat W. Intraspecific variation of Parasibirites kolymensis Bychkov (Ammonoidea) from the Lower Triassic (Spathian) of Arctic Asia. Mitteilungen aus dem Geologisch-Paläontologischen der Institut Universität Hamburg. 1999;83:163-178.

[13] O'Keeffe SG, Moulton DE, Waters SL, Goriely A. International Journal of Non-Linear Mechanics. International Journal of Non-Linear Mechanics. 2013 Nov;56(C):94-104.

[14] Rodriguez EK, Hoger A, McCulloch AD. Stress-dependent finite growth in soft elastic tissues. Journal of biomechanics. 1994;27(4):455-467.
[15] Howell P, Kozyreff G, Ockendon J. Applied solid mechanics. vol. 43. Cambridge University Press; 2009.

[16] Checa A. A model for the morphogenesis of ribs in ammonites inferred from associated microsculptures. Palaeontology. 1994;37(4):863–888.