Mathematical models of neuronal growth

The establishment of a functioning neuronal network is a crucial step in neural development. During this process, neurons extend neurites—axons and dendrites—to meet other neurons and interconnect. Therefore, these neurites need to migrate, grow, branch and find the correct path to their target by processing sensory cues from their environment. These processes rely on many coupled biophysical effects including elasticity, viscosity, growth, active forces, chemical signaling, adhesion and cellular transport. Mathematical models offer a direct way to test hypotheses and understand the underlying mechanisms responsible for neuron development. Here, we critically review the main models of neurite growth and morphogenesis from a mathematical viewpoint. We present different models for growth, guidance and morphogenesis, with a particular emphasis on mechanics and mechanisms, and on simple mathematical models that can be partially treated analytically.


Introduction
The human brain contains an estimated ∼ 86 billion neurons (Azevedo et al. 2009) that interconnect through numerous slender processes, axons and dendrites, globally known as neurites. A typical nerve cell has one long axon, that transmits nervous signals, and multiple branching dendrites, that receive signals (Striedter 2016, p. 35), see Fig. 1. During early neuron development, multiple neurites sprout from the cell body (or soma). During this phase, typically, one incipient neurite comes to outgrow the others and differentiates into an axon (Kiryushko et al. 2004). This axon then migrates through the extracellular matrix in order to make synaptic connections with other distant neurons or tissues.
It is now understood that this fundamental and robust process depends on many different genetic, biochemical and physical factors. Mathematical modeling has become a prominent actor of developmental biology and is crucial to achieve a rational organization of the intricate physical processes involved in neuronal development. Here, we provide a critical, but not exhaustive, review of about three decades of mathematical modeling of axon and dendrite development. This paper is meant to complement existing reviews that emphasize physical and chemical mechanisms of growth and guidance (O'Donnell et al. 2009;Franze and Guck 2010;Franze 2020;Seo et al. 2020;Suter and Miller 2011;Miller and Suter 2018;Mortimer et al. 2008;McCormick and Gupton 2020;Abuwarda and Pathak 2020;Dickson 2002;Franze et al. 2013;Dent and Gertler 2003). computational and numerical models (Simpson et al. 2009;Kiddie et al. 2005;Graham and Van Ooyen 2006;Van Ooyen 2011Ascoli 2002;Maskery and Shinbrot 2005;Goodhill 2018), or mechanical models (Goriely et al. 2015a). Our particular emphasis is on analytically tractable mathematical models that provide key insights into possible mechanisms and properties emerging from nonlinear couplings, and the combination of different physical phenomena such as material transport, growth and mechanics.

Biological background
Neurite motility is mediated by the growth cone (Fig. 1), a highly dynamic actin-supported extension at the tip of the neurite that performs both sensory and locomotory functions. Within the growth cone, the cytoskeleton is mostly constituted of actin filaments (F-actin) that organize in a lamellipodium-like structure. Centrifugally arranged actin filaments polymerize at the periphery of the growth cone (P-domain) and undergo a treadmill-like retrograde flow mediated by myosin II motors in a transition zone (T-domain), where actin depolymerizes. Combined with focal adhesions that mechanically couple F-actin with the substrate (Chan and Odde 2008), this continuous movement generates forces that propagate through a connecting region (C-domain) to the trailing part of the neurite, the neurite shaft (Franze 2020), and mediate neurite motility.
The neurite shaft is a filament-like section that contains the neurite cytoplasmic compartment (Kevenaar and Hoogenraad 2015), composed mainly of stabilized microtubules cross-linked by microtubule-associated proteins (MAP) (Maccioni and Cambiazo 1995), F-actin and neurofilaments. Microtubules are vital for intracellular trafficking (Caviston and Holzbaur 2006) and for the neurite's structural integrity. The microtubule array is itself coated with an actomyosin cortical sheath generating contractile stresses (Franze 2020;Jülicher et al. 2007) that under certain circumstances may cause neurite retraction (Recho et al. 2016;Franze et al. 2009a). The inner and outer layers of the neurite shaft are connected through special proteins ensuring force transmission between the two structures (Coles and Bradke 2015).
We define neurite growth as the irreversible elongation of the neurite shaft supported by addition of new cellular material (Goriely 2017;Goriely et al. 2015a). The precise mechanisms underlying neurite growth have been debated for about 40 years and remain controversial. An important idea in the 80s-90s was that microtubules are stationary during neurite extension (Miller and Joshi 1996;Lim et al. 1990;Okabe and Hirokawa 1990;Takeda et al. 1995;Sabry et al. 1995) and assemble at the neurite distal tip to elongate the shaft (Dent and Gertler 2003;Goldberg and Burmeister 1986;Bamburg et al. 1986), a process regulated through growth cone traction forces Heidemann 1988, 1992;Lamoureux et al. 1989;). Thus, a vast literature has focused on slow axonal transport (Miller and Heidemann 2008;Roy 2014Roy , 2020, the general process by which cytoskeletal proteins, such as neurofilaments, actin and tubulin (the building unit of microtubules), transit along the neurite. The mechanisms of slow axonal transport have been debated for decades and remain relatively poorly understood (for a review of models, see Bressloff and Newby 2013).
To date, direct observation of this tip growth process is lacking. Moreover, this hypothesis seems in apparent contradiction with the finding that neurites isolated from a substrate may stretch rapidly when subject to externally applied deformation or force, a process which does not seemingly involve the growth cone (Dennerll et al. 1989;Bray 1984;Rajagopalan 2010;Zheng et al. 1991;Lamoureux et al. 2010;Heidemann et al. 1997).
The idea that the shaft is stationary has been dramatically challenged by more recent investigations which have clearly evidenced that the neurite shaft stretches during growth cone migration (Miller and Suter 2018;Athamneh et al. 2017;Lamoureux et al. 2010;O'Toole et al. 2008a;Miller and Sheetz 2006). These observations are partly consistent with a viscoelastic creep of the shaft under prolonged traction, mediated by adhesion (O'Toole et al. 2008a). It is now generally understood that neurite growth is largely controlled by Fig. 1 Schematic of a typical neuron showing the soma, the dendrites and the axon with its growth cone, prior to synaptogenesis. Close-up shows a schematic representation of the cytoskeleton in a migrating neurite. The shaft (on the left) is essentially composed of cross-linked microtubules surrounded by an actomyosin sheath. The growth cone is a lamellipodium-like structure that confers motility to the neurite, probes its environment using sensory filopodial protrusions and produces traction forces applied to the shaft promoting neurite extension mechanics and depends on the complex interplay of active pushing and pulling forces mediated by growth cone traction, microtubule assembly, shaft viscoelasticity and contractility, and adhesion (Recho et al. 2016;Franze 2020;Miller and Suter 2018;Suter and Miller 2011).
To reach their functional targets, axons must extend along precise paths. This early phase, known as axon guidance (or pathfinding), crucially relies on cues from the environment which are processed to instruct growth cone trajectory. This includes chemotaxis primarily, i.e., guidance of axons by gradients of diffusing chemicals like Slits or Netrins (Mortimer et al. 2008;Plachez and Richards 2005;Chilton 2006); but also haptotaxis, i.e., guidance by gradients in adhesion or substrate-bound chemicals (Sundararaghavan et al. 2011); durotaxis (or mechanotaxis), i.e., guidance by gradients in substrate stiffness (Abuwarda and Pathak 2020;Koser et al. 2016;Espina et al. 2021); electrotaxis (or galvanotaxis), i.e., guidance by electric field (Yao and Li 2016;Hamid and Hayek 2008;Gokoffski et al. 2019;Shapiro et al. 2005); curvotaxis, i.e., guidance by substrate curvature (Smeal et al. 2005); or guidance assisted by guidepost cells such as radial glial cells (Franze et al. 2009b;Rakic 1972) or Schwann cells (Thompson and Buettner 2006). The growth cone is the main sensory structure at play in all these guidance modalities.
Once an axon has reached its target, typically neurons or other target effector cells, it forms synapses (synaptogenesis) and is therefore bound to these cells. Then, the axon shaft must extend to accommodate the subsequent growth of the surrounding tissue. In this stretch growth phase, spectacularly fast axon extension may take place, dictated by the growth of the animal's body (Smith 2009) and governed by forces (Pfister et al. 2004). It is now understood that maintaining a proper tension along the neurite is important in normal function (Siechen et al. 2009). Body growth imposes strong kinematic constraints on the neuron, which must then sustain sufficient cell material supply while performing fast remodeling along the stretched neurite in order to maintain its structural integrity and avoid traumatic levels of strain. Several authors have modeled and studied this problem, e.g., O'Toole and O'Toole et al. (2008b); Purohit and Smith (2016), who have examined the role of stretching in controlling the flux of material along the shaft (see Sect. 3.3). However, this problem remains poorly understood; yet, it is nonetheless crucial for instance in the design of nerve regeneration therapies.

Modeling neuritic growth
Neurite growth involves different coupled mechanisms that have not been fully resolved. Initially, growth was approached from the point of view of transport, mass supply and addition, i.e., as a process entirely governed by the availability of given substances required for growth. In particular, many authors have focused extensively on the transport and assembly of tubulin at the tip of the neurite, long viewed as the main growth determinant. Alternatively, experiments have clearly shown that neurons grow in response to applied mechanical forces. It is then natural to model growth in terms of forces, strains, stresses and rheology. Here, we consider both approaches, starting with transport models.

Transport-limited growth
By definition, growth is related to the notion of mass uptake (Goriely 2017). One of the main growth determinants is soluble tubulin that was assumed to assemble at the end of the neurite. In young neurites, the new cell material is supplied by the soma. Therefore, the cell must first produce new tubulin units (and other cell constituents) and then carry them along the shaft, using diffusion and active transport, to reach the end of the neurite where growth takes place. This process is at the core of numerous mathematical works (see the reviews by Kiddie et al. 2005;Graham and Van Ooyen 2006;Van Ooyen 2011.

Zero-dimensional models
A first approach consists in modeling the neurite by a small number of compartments (Janulevicius et al. 2006;Samuels et al. 1996;Purohit and Smith 2016;Van Veen and Van Pelt 1994;Lin et al. 2020). Each compartment is associated with one or more chemical concentrations and exchanges material with other compartments via diffusion and/or active transport by motor proteins. For instance, in a model such as the one proposed by Van Veen and Van Pelt (1994), the neurite is modeled by two compartments, the soma, where the tubulin dimer concentration c 0 is defined, and the growth cone with tubulin concentration c 1 (Fig. 2a). The two compartments are virtually separated by a distance that represents the neurite length.
Growth relies on assembly and disassembly of tubulin at the neurite tip, following a polymerization reaction of the type where T n denotes a microtubule of length n for n > 1, or a tubulin dimer if n = 1 ; and k + and k − are rate constants. This process implies the kinetic law where e represents the typical elongation due to the addition of one tubulin dimer (in unit length per mole of tubulin). Growth is coupled to c 1 since it consumes tubulin ( k + ), and conversely, disassembled tubulin is reintroduced ( k − ) in the available pool of tubulin at the tip. In addition, tubulin is supplied by the soma through diffusion. A Fickian flux J is assumed, where tubulin diffuses with diffusion constant D as with A the cross-sectional area of the neurite's cytoplasmic compartment. Summing up all contributions, we obtain with V a typical volume where reaction takes place. To close the system, we take into account the production of tubulin by the soma with rate S: where we have used the same volume V introduced previously. Choosing the time, length and concentration units to be, respectively, A/D, SAe/D and SA/VD, we obtain the non-dimensionalization where = AD∕Se , = k − ∕S and = Ak + ∕VD are dimensionless. For simplicity, we henceforth drop the tildes and work with dimensionless variables. We stress that the length of the neurite is a dynamical variable, and not a spatial dimension of the problem. Hence, we have an autonomous system of three nonlinear ordinary differential equations for the variables c 0 , c 1 and , which can be analyzed using dynamical systems theory. This system has an exact solution given by with see dashed lines in Fig. 2b-d. This solution also fully captures the asymptotic dynamics, and for positive initial conditions we have that as t → ∞ , c 0 ∼ c 01 t , ∼ vt and c 1 → c 00 . Another interesting feature of the model that can be observed from (8) is that for small initial concentrations c 1 (0) < ∕ and c 0 (0) > c 1 (0) , i.e., when disassembly dominates over assembly, we have ̇ < 0 , and hence, we have a transient neurite collapse. However, under the same conditions, we have ̇c 1 > 0 . Hence, c 1 will increase in time until ̇ ≥ 0 at which time the growth of the neurite resumes.
Other parameter values: 0 = 0.1 , = 10 , = 1 and = 1 A drawback of the solution (9) is that both soma concentration c 0 and length increase without bound. In reality, tubulin dimers may inhibit the translation of new tubulin (Gay et al. 1989), and thus, cells have a limit on c 0 and may generally use a feedback loop to control it. It is also clear that the transported proteins have finite half-lives, which imposes a theoretical limit on neurite length (Miller and Samuels 1997;).
An interesting extension of the two-compartment model is the application to multiple neurites Samuels et al. 1996;Toriyama et al. 2010;Fivaz et al. 2008) or neurites with multiple branches (Li et al. 1992;Van Veen and Van Pelt 1994;Hjorth et al. 2014). For the sake of illustration, consider N > 1 neurites growing from the same soma, as illustrated in Fig. 2a. Defining c i to be the concentration at the tip of the ith neurite of length i , the flux for neurite i is defined following Samuels et al. (1996) as where we have added a growth-independent transport controlled by a parameter , and a growth-dependent transport controlled by a parameter . The last term models a growthsensitive feedback loop in which a faster-growing neurite receives a larger influx of tubulin. Hence, the N-compartment model leads to a system of 2N + 1 differential equations for c 0 , c i and i . The last term in (13) introduces nonlinearity and possible instabilities. Indeed, this new augmented system is very sensitive to the choice of initial lengths (but relatively insensitive to differences in initial concentrations). Since neurites arise from the soma at different times, we initiate the simulation with different initial lengths i (0) . For large values of the feedback parameter (here = 100 ) and neglecting growth-independent transport ( = 0 ), a winnertakes-all instability is observed, where one neurite comes to outgrow the others by monopolizing tubulin resources. An example of such phenomenon is shown in Fig. 3b where we have simulated N = 4 neurites with initial lengths 0 = 1 , 0.99, 0.98 and 0.97. This mechanism provides insight into the mechanism by which a neurite differentiates into an axon during early neuronal development. Note, however, that as soon as > 0 , the shortest neurites, whose growth is initially inhibited by the dominant neurite, catch up at later times as concentrations become homogeneous across the N tips, as shown in Fig. 3c The prediction that the initially longest neurite almost systematically becomes the axon is probably inaccurate in reality, as exemplified by situations where an initially shorter neurite outgrows longer ones (see, for instance, Fig. 6 in Fanti et al. 2008). In fact, the initial length of neurites is one among many other endogenous and exogenous factors that may participate in breaking the symmetry. Nevertheless, this simple model, based on the classic concept of patterning by autocatalysis and lateral inhibition (Turing 1952;Gierer and Meinhardt 1972;Meinhardt and Gierer 2000), shows how a simple competition mechanism can trigger symmetry breaking and prevent the formation of multiple axons (an idea further supported by later investigations, e.g., Fivaz et al. 2008;Toriyama et al. 2010;Takano et al. 2015Takano et al. , 2017Inagaki et al. 2011).

One-dimensional models
The main postulate of compartmental approaches is that the concentration profile along a neurite is linear and defined by only two quantities c i and c 0 at both ends of the neurite. For a long neurite such as an axon, diffusion and transport are not instantaneous, and proteins may degrade significantly before reaching the tip. Hence, delay and spatial heterogeneity in cellular material concentrations will appear along the neurite. Thus, the coarse two-compartment approximation (3) needs to be refined by considering the value of the concentration at every point along the neurite. A possible generalization of the compartmental approach consists in considering neurites made up of many concatenated compartments representing short contiguous portions of the neurite shaft, where each segment exchanges solutes with its adjacent neighbors (Hely et al. 2001;Kiddie et al. 2005;Van Ooyen 2001, 2004;Hjorth et al. 2014). Bulk growth can be simulated by growing and subdividing each compartment. Alternatively, new compartments may be added at the tip, to simulate tip growth. Mathematically, it is interesting to consider the continuum limit of these approaches through a proper partial differential equations formulation of the fundamental convection-diffusion problem, and then to use numerical discretization methods when needed (Smith and Simmons 2001;Graham 2004, 2006;Diehl et al. 2014Diehl et al. , 2016Zadeh and Shah 2010). Figure 4, reproduced from , gives the geometry of the problem and the principal mechanisms involved. The continuity equations express conservation of each substance considered (e.g., proteins, vesicles, mitochondria). Here, we consider one diffusing substance (e.g., tubulin) with concentration c obeying In this expression J is the anterograde flux of material (from the soma and oriented toward increasing x), and S represents any source or sink. The flux J = J d + J a includes the regular diffusion J d = D c∕ x with D the molecular diffusivity, and a kinesin-mediated flow J a = ac with a an effective active transport velocity. Typically, transported proteins may also decay with rate constant , which is included into a sink term S = − c . In 1 3 the absence of diffusion ( D = 0 ), (14) supports a steady state solution c = c 0 exp(− x∕a) , where c 0 is the concentration of tubulin in the soma. This solution identifies a characteristic length of a∕ that may be interpreted as a maximal length for the neurite (Miller and Samuels 1997).
Several boundary constraints have been used. For the proximal boundary conditions, we can specify at the soma, either the concentration c| 0,t =c 0 (t) Diehl et al. 2014Diehl et al. , 2016   where and are constants related to assembly and disassembly of tubulin, respectively. Growth is then modeled by adapting (2) in terms of the Stefan condition, which expresses the displacement of the distal boundary of the domain, depending on the availability of tubulin at the neurite tip. Moving-boundary problems are notoriously difficult to solve. However, in one dimension, the domain size can be fixed through the change of variable y(x, t) = x∕ (t) ∈ [0, 1] . Equation (14) is then recast as with c(y, t)∶=c( y, t) for all y ∈ [0, 1] , and where (t) is governed by (16). The full time-dependent problem is analytically intractable in general, and numerical approaches are necessary (see, for instance, the simulation by Diehl et al. 2016, shown in Fig. 5). However, focusing on the steady regime, ; ;  show that the neurite reaches a uniquely defined final length when tubulin degradation is considered (see also Miller and Samuels 1997) and grows without bound otherwise (as in Van Veen and Van Pelt 1994). Overall, these 0D and 1D models generally predict a decrease in tubulin concentration toward the tip, as in the 0D Van Veen-Van Pelt model or the 1D McLean-Graham approach. It can be noticed, however, that the boundary condition (15) used by  actually neglects the moving boundary and is not consistent with the expression of the flux. To resolve this issue, Diehl et al. (2014) represent the growth cone as a finite volume with tubulin concentration c gc , where tubulin polymerization is modeled separately through an ordinary differential equation. Therefore, the distal boundary constraint consists in enforcing the continuity of c at the interface, i.e., c( (t), t) = c gc (rather than a flux). This modification of the original model results in very different, non-monotonic, concentration profiles that show an accumulation of tubulin near the tip when transport is faster than elongation (Fig. 5). Experimentally, this is reminiscent of the distal accumulation of MAPs reported by Black et al. (1994), but is, however, inconsistent with the nearly uniform concentration of mitochondria observed by O'Toole et al. (2008b) (see Sect. 3.3) that may be explained by a slowdown of transport toward the tip (as also suggested by Hoffman et al. 1985;Watson et al. 1989).
The general approach presented here is based on the idea that elongation is governed by microtubule assembly at the tip, and that high tubulin concentration causes faster growth (16). This assumption has not been directly demonstrated experimentally and is partly challenged by recent observations, such as those made by Ren and Suter (2016) in pausing Aplysia neurites (see also the comment by Miller and Suter 2018). When the growth cone pauses, it rapidly increases in cell content due to transport, which correlates with a rapid increase in growth cone size. Yet, the growth cone neck does not advance, indicating that neurite growth is limited by growth cone advancement, rather than transport and microtubule assembly. In addition, this view neglects other rate-limiting constituents such as MAPs, neurofilaments, mitochondria or membrane lipids, as well as mechanical forces. Several models have addressed the modulation of tubulin assembly by forces, such as Heidemann (1988, 1992), or, more recently, Purohit (2015), who modified the Van Veen-Van Pelt model to include mechanical regulation of tip growth. In the next section, the focus is on the mechanical models of neurite growth.

Mechanically mediated growth
An important aspect of neuronal development involves regulation through physical forces (Mutalik and Ghose 2020; Franze 2020; Suter and Miller 2011;Franze et al. 2013;Miller and Suter 2018). In particular, the initiation and growth of neurites may be induced directly by forces applied on the neurite Bray 1984;Zheng et al. 1991;Dennerll et al. 1989). For example, Zheng et al. (1991) showed that the initiation of axons of chick sensory neurons can be triggered experimentally by the application of tension on the surface of the cell body. Then, further elongation can be induced by towing the neuron with a glass needle. The deflection of the needle tip can be used as a control to apply constant force. Strikingly, these experiments demonstrate that growth rates increase with applied forces. While it remains unclear whether neurites undergo tip growth, an alternative paradigm has emerged more recently, where forces drive the fluid-like expansion of the neurite shaft, while transport provides the material necessary to support elongation through mass addition (Miller and Suter 2018;Suter and Miller 2011).
To model the role of forces acting on the neuron and its internal components, we start with a system composed of a single isolated unit (zero-dimensional models) subject to forces, before considering more complex models with multiple components and with spatial variations in one dimension. In contrast to Sect. 3.1, the problem of material supply will be considered here as non-limiting.

Zero-dimensional models
It is well appreciated that isolated neurites respond elastically under rapid plucking (Heidemann et al. 1997), and creep irreversibly without thinning when subject to prolonged overstretch (Dennerll et al. 1989;Bray 1984;Heidemann et al. 1997;Pfister et al. 2004). Typical experimental studies of mechanically mediated neurite growth consider isolated neurites subject to imposed displacement or force. To model the mechanical response of neurites in such condition, spring-anddashpot-type models have been used (Dennerll et al. 1989;Bernal et al. 2007;Li and Qin 1996;O'Toole et al. 2015;Lin et al. 2020) where the entire neurite is represented as a single unit with no spatial variation. A more detailed approach emphasized here consists in viewing the neurite as a morphoelastic tubular compartment (Anthonisen and Grütter 2019; Goriely et al. 2015a;Goriely 2017;Wang and Kuhl 2019;Holland et al. 2015). Morphoelasticity (Goriely 2017) is an extension of nonlinear elasticity which provides a general mechanical description of growth based on a multiplicative decomposition of the deformation gradient (Rodriguez et al. 1994). It is a natural framework to represent nonlinear deformations along the neurite, as well as more complex multidimensional processes such as axon beading (Riccobelli 2021) or turning (García-Grajales et al. 2017). In particular, morphoelasticity naturally takes into account changes in crosssectional area, which has been used experimentally as an indicator of growth (O'Toole et al. 2008a;Lamoureux et al. 2010), as well as arbitrary nonlinear elastic response functions. The application of morphoelasticity to axon growth, elasticity and damage is extensively discussed in Goriely et al. (2015a). In particular, one core idea of morphoelastic modeling that differs in spirit with modeling through spring and dashpots is the particular focus on biological mechanisms. The extra variables introduced in the theory represent effects such as remodeling and addition of mass. The overall behavior of such systems under loads, namely its rheology, is then an emergent property rather than a choice of constitutive law based on experimental observations.
A neurite is represented initially by a homogeneous stress-free cylinder 0 of cross-sectional area A 0 and length L 0 (Fig. 6a). We assume that the neurite is fixed at one end (defined as the soma). At a later time t, the current (observed) configuration of the neurite is another cylinder with cross section a(t) and length (t) . We denote by 0 = X 0 , Y 0 , Z 0 ∈ 0 and = (x, y, z) ∈ the positional vectors in the initial and current configuration, respectively. Since we do not allow spatial variations, the deformation gradient ∶= ∕ 0 is uniform across the whole domain and is given by the tensor = diag , ⟂ , ⟂ , with and ⟂ the longitudinal and transverse stretches.
The main postulate of morphoelasticity consists in expressing the deformation gradient in terms of two tensors: a growth tensors due to the change of volume or internal reconfiguration, which is stress-free, and an elastic tensor , that characterizes the elastic response, and generates stress (Rodriguez et al. 1994). Constitutively, these two tensors are multiplicatively coupled via where = diag , ⟂ , ⟂ and = diag , ⟂ , ⟂ . The growth deformation transforms the initial cylinder 0 into a new intermediate cylinder g (t) of length L(t) = (t)L 0 and cross-sectional area A(t) = ⟂ (t)A 0 , whereas the elastic part of the deformation transforms the grown cylinder g (t) into the observed cylinder (t) with = ⟂ (t)L(t) and a(t) = ⟂ (t)A(t) (Fig. 6a).
A possible kinetic law that captures the stretch-induced growth of the neurite is the Bingham-type relation where k is a kinetic constant, and c ≥ 1 is a critical stretch above which growth may occur, as enforced by the ramp function (⋅) + = max (0, ⋅) (Anthonisen and Grütter 2019; Goriely 2017). Below this threshold, deformation is purely elastic. For simplicity, we further assume that the neurite maintains a homeostatic cross-sectional area A = A 0 in the grown configuration, which implies ⟂ = 1 . Kinematically, ̇∕ is a natural definition for the rate of growth, related to the spatial velocity gradient of continuum kinematics Holzapfel 2000), and is the natural choice for growth processes that occur uniformly in a region. For growth processes that are restricted to a physical region the choice of ̇ rather than ̇∕ in (19) would be a better modeling choice. Physically, (19) captures the stretch-and-intercalation process that has been proposed as a possible mechanism for bulk growth (Suter and Miller 2011;Holland et al. 2015), where the product k − c + phenomenologically encompasses thermodynamical processes and cellular regulation pathways involved in intercalation which we ignore here. We note that similar laws have been used to describe the straininduced growth of plant cell walls (Lockhart 1965;Boudon et al. 2015;Goriely et al. 2008;Ambrosi et al. 2019;Ali et al. 2014;Bozorg et al. 2016;Zhao et al. 2020). Next, we need to specify the geometric constraint or load applied to the system. This choice is motivated by the experimental setting or biological scenario considered. Experimentally, various methods have been used to control either the length of the neurite (Bray 1984;Purohit and Smith 2016) or the applied tension (Dennerll et al. 1989;Heidemann et al. 1997). In vivo, we may assume that integrated neurites bound to the body are subject to an imposed displacement, insofar as their internal tension remains small and does not itself affect the growth of the embedding medium. In this case, (t) is a prescribed function, and we may solve (19) using the identity = ∕ : where ∶=1∕k c defines a timescale for growth, and where we have taken (0) = 1 and > c . For a ramped traction of the form (t) = 0 + t∕ � , with 0 a constant and ′ a characteristic time of traction, (20) can be integrated to obtain 6 Zero-dimensional morphoelastic model of neurite growth a Multiplicative decomposition of the stretch from initial (stressfree) configuration (with length L 0 ), to intermediate, stress-free configuration (with length L), to current configuration (with length ). The multiplier accounts for the anelastic deformation of the neurite, while quantifies the elastic stretch that results in mechanical stress. b Cartoon illustrating the localized growth model where growth takes place only in the axon hillock. The basal zone (red) is the only part that grows. Growth is frozen when the considered mate-rial point exists the hillock. c Stretch-and-hold experiment (t) = 0 . Plot of (t) , (t) and (t) versus normalized time k c t . All stretches are normalized as x → (x − 1)∕ 0 − 1 . d, e Speed-controlled traction experiment (t) = 1 + t∕ � simulated with both d the exponential growth model (19) and e the linear growth model (30). Log-log plot of (t) − 1 , (t) − 1 and (t) − 1 vs. normalized time k c t . Dashed lines show the asymptotic approximations. Parameters: c = 1.1 , In a stretch-and-hold experiment with 0 ≥ c and � → ∞ , the axon is instantly stretched from its stress-free configuration L 0 to a fixed length 0 L 0 . In this case, the neurite relaxes exponentially (Dennerll et al. 1989), see Fig. 6c, and reaches a homeostatic stretch given by the critical yield stretch c .
Alternatively, in a speed-controlled experiment, 0 = 1 with ′ finite, the axon is pulled with speed L 0 ∕ � , as simulated in Fig. 6d. In this case, growth is asymptotically linear (t) ∼ t∕ c � , and as before, the excess stretch ( − c ) vanishes as t → ∞ . Note, however, that for c > 1 , the onset of growth is not at t = 0 but is delayed to t = c − 1 � .
Next, we consider a force-controlled experiment where the neurite is pulled with a given longitudinal external force F. In that case, we need to relate the stresses to the strains. To do so, we assume that, instantaneously, the material is hyperelastic and incompressible. It is then described via a Gibbs free energy density ( ) (in unit energy per intermediate volume). Incompressibility implies ⟂ = 1∕ √ , from which we introduce the auxiliary energy function In normal growth conditions, the timescale of elastic relaxation ( ∼minute) is much shorter than that of irreversible expansion, which depends on relatively slow biophysical processes ( ∼hour). Therefore, the system can be considered at quasi-static equilibrium at all times. The equilibrium condition follows from the principle of virtual work where W ext = F = FL is the virtual work due to the traction force F > 0 , and W int = − (LÂ ) is the virtual work due to the reaction force. Thus, the equilibrium configuration corresponds to a root of A variant of the stretch-and-hold experiment consists in pulling the neurite tip with a micro needle perpendicular to the neurite and whose base is fixed at a position x = L n ≥ L 0 .
Modeling the needle as a cantilever beam with length D and flexural rigidity B, we obtain the condition Taking B → ∞ gives the stretch-and-hold scenario seen before. However, the small deflection of the needle here provides a measure of the neurite tension, given by 3B∕D 3 L n − (Zheng et al. 1991). Another popular experiment consists in pulling the neurite laterally (Bernal et al. 2007;Rajagopalan 2010;Dennerll et al. 1989). For the sake of brevity, we ignore this last case (see Goriely et al. 2015a, for details).
Lastly, we can also include the active properties of the neurite due to the actomyosin sheath that generates contractile forces. Several methods allow to model active effects in elasticity (Goriely 2018). For our problem, a possible approach consists in augmenting the balance equation (23) with an active stress T a < 0 so that Thus, upon removal of the load, i.e., when F = 0 , the neurite is at equilibrium in a contracted configuration < 1.
Several forms of the strain-energy function ̂ have been discussed in Goriely et al. (2015a). For example, we consider a neo-Hookean material with shear modulus . The energy function is ̂ ( ) = ∕2 2 + 2∕ − 3 , and (23) reduces to the cubic 3 − F∕A 0 2 − 1 = 0 which supports only one real solution for the stretch entering (19). For a constant F, is constant and growth is exponential.
In reality, towed neurites grow at a nearly constant speed depending on the applied force (Bray 1984), but this remains to be assessed on longer time scales. To explain this possible discrepancy, we stress that, in postulating (19), we assumed that growth takes place at all points of the shaft. The precise location of growth in stretched axons has not been fully elucidated (Futerman and Banker 1996;Smith 2009). For instance, in the compartmental model developed by Purohit and Smith (2016) (discussed more in details in Sect. 3.3), growth of stretched axons is assumed to take place via polymerization reactions localized in the axon hillock, i.e., the proximal segment of the neurite close to its junction with the soma. To describe this mechanism, we assume that growth happens only in a physical zone of constant length L h in the intermediate configuration, located near the proximal boundary, as shown in Fig. 6b. In this formulation, the hillock is not a Lagrangian material volume in the sense that its length changes with respect to the reference configuration 0 . We modify (19) to account for this new assumption, and we posit a localized exponential growth law: where X depicts the axial coordinate in the intermediate configuration, and H is the Heaviside step function enforcing that growth only occurs for X ≤ L h . Note that, in contrast to , is still uniform along the whole shaft provided a uniformly-defined strain energy function. Thus, the resting length of the neurite L(t) obeys which is independent of L(t) . For a speed-controlled experiment, with pulling velocity d ∕dt = L 0 ∕ � , the elastic stretch tends to a small residual value of order for t ≫ ′ ≫ , i.e., when pulling is slow; and becomes large, when t ≫ ≫ ′ , i.e., when pulling is fast, which may lead to disconnection. This is consistent with Purohit and Smith (2016) (see Sect. 3.3) in which axon tension plateaus at a finite value that depends on the pulling speed.
In a force-controlled experiment with tension F, the system develops a constant stretch (23), and its size increases with constant speed ̇ = L h k − c . Although, locally, the growth mechanism itself is the same, this regime is very different from the exponential growth predicted by the previous growth model described by (19). This is due to the fact that the expanding zone L h remains constant (see Goriely 2017, pp. 69-70).
Interestingly, a popular growth law (Wang and Kuhl 2019;Holland et al. 2015;Goriely et al. 2015a) to describe stretched axons is An example simulation for a speed-controlled experiment is shown in Fig. 6e. Since is here uniform, the resting length increases as which is equivalent to (27) up to a rescaling. Thus, (30) rheologically captures the mechanical response of a neurite undergoing accretive or localized growth. Note, however, that, to our knowledge, there is no clear proof of the existence of such growth regime. In fact, every segment of the axon has the potential to undergo growth (Lamoureux et al. 2010). The hillock growth model proposed by Purohit and (27) Smith (2016) is based on the observation that axons elongate mostly near the cell body, when the latter is towed, while the growth cone is held fixed with respect to the substrate (reverse towing experiment detailed in Lamoureux et al. 2010). However, this argument is not fully satisfactory, as Lamoureux et al. (2010) considered neurites embedded on a substrate, therefore not mechanically isolated. As shown in the next section, this second scenario is different, and, in this case, mechanics demands that the axon grows faster near the soma even when the whole shaft is assumed to expand. In conclusion, the exact location of growth in stretched axons and a correct morphoelastic representation for it remain an open problem.

One-dimensional models
The 0D models are suitable for neurites grown in controlled experimental conditions where longitudinal tension is homogeneous. However, in elongation mediated by the growth cone, this assumption will not hold due to the presence of dissipative forces such as cell-substrate adhesion, and one needs to consider variations of both stress and velocity along the neurite. In fact, it has been clearly shown that migrating presynaptic neurite stretch non-uniformly while elongating (Miller and Suter 2018;Athamneh et al. 2017;Lamoureux et al. 2010;O'Toole et al. 2008a).
A simple and conceptually interesting model has been first introduced by O' Toole et al. (2008a) (and adapted to morphoelasticity by Goriely 2017, see pp. 80-85). In this model, the growing axon is seen as a simple fluid structure in contact with a substrate. Here, we focus on the permanent self-similar regime, for which the stress and velocity profile are constant in time in the vicinity of the neurite tip (Oliveri et al. 2021). Therefore, it is convenient to work in the advected frame attached to the tip, and to parameterize the neurite by the coordinate x oriented toward the soma and originated at the neurite tip (Fig. 7, inset).
At any position x, we may virtually cut the neurite to define n(x, t) the force applied by the cross section x + to the cross section x − . The classic procedure to derive the local balance equation based on the balance of linear momentum (Goriely 2017) provides an equation for n: where f depicts the lineal density of tangential body force applied along the shaft. In order to model the interaction between the shaft and the substrate, we assume a body force due to damping of the form with v the anterograde velocity in the laboratory frame, and a friction coefficient. In addition, we postulate the constitutive relation expressing the stress-induced growth of the shaft, with a bulk growth parameter analogous to a viscosity. This constitutive law is akin to the exponential growth law seen previously (Sect. 3.2.1) where plays a role similar to ∕k . Note that the minus sign results from the choice of orientation for the x-axis, opposite to the anterograde motion.
On differentiating (32) and plugging (33, 34), we obtain: where h = √ A ∕ defines a characteristic length for the problem. The distal end ( x = 0 ) of the neurite is subject to a tensile load F, while tension vanishes for x ≫ h , as energy dissipates through friction. This dictates the boundary conditions n| (t),t = F, n| ∞,t = 0, and the solution for the tension and the velocity profiles shown in Fig. 7, which illustrates the effects of both parameters and . We have in particular the two limits that, respectively, correspond to a stalling regime, when neurite viscosity is high, and a rupture regime, when viscosity is low, with a strain at the tip dv∕dx → ∞ . The adhesion both reduces the velocity of the neurite and the propagation of forces along the shaft. We stress that the tip speed v( (t), t) = F∕ √ A is constant in time since growth happens in an effectively finite zone a, as discussed in Sect. 3.2.1, which for a small, results effectively in a tip growth regime. Experimentally, the tension profile can be measured from the spatial variation of the cross-sectional area. Indeed, assuming an incompressible neurite, we compute An extension of this model was developed by Recho et al. (2016) using a similar formalism. In this extended model, the neurite is composed of three compartments, shown in Fig. 8a. The microtubule core is represented by an elastic solid undergoing mass depletion and mass addition that aim to maintain a density at chemical equilibrium. Namely, a compressed core will have an excess density, compensated by a mass removal. Conversely, mass addition allows to compensate for the decrease in density due to an excessive traction. This mechanism is coupled with an elastic constitutive behavior that provides the pressure of the core. The core is itself embedded inside a viscocontractile actomyosin cortex. The two compartments are mechanically coupled via frictional interactions. A no-flux boundary condition at the junction between the microtubule core and the growth cone implies the absence of tip growth. The structure formed by these two compartments is pulled at its distal end by a third compartment that represents the advancing growth cone made up of actomyosin and pulling on the substrate. Finally, the dynamics of the structure is obtained from first physical principles, namely the balance of mass and the balance of momentum.
The resulting model supports three states, shown in Fig. 8b. A retraction regime occurs when microtubule pressure and growth cone traction are insufficient to overcome the actomyosin contractile forces. Conversely, relatively In this state, the structure neither grows nor retracts unless the parameters promoting either regimes exceed a threshold value. This scenario is reminiscent of the seminal "rubber band gun" model proposed by Buxbaum and Heidemann (1988), who considered microtubules subject to longitudinal compression due to the front membrane of the neurite. Traction forces produced by the growth cone reduce compression, decreasing the free energy required for microtubule assembly and growth. The power of Recho et al.'s model is that it addresses each structural compartment of the neurite separately and relates their interactions to parameters that can be modified experimentally. Indeed, it is able to reproduce in details the observed qualitative effects induced by drugs, such as blebbistatin, that reduces myosin II motor contractility; cytochalasin, that inhibits actin polymerization; trypsin, that reduces substrate adhesion; nocodazole, that destabilizes microtubules; or taxol, that stabilizes microtubules. These drugs affect distinct and well-identified target parameters of the models, and their respective effects can be simulated specifically. This also explains coupled effects that can appear as paradoxical when looked at in isolation.
The models presented so far are based upon continuum formulations, which allows us to make use of well-known mathematical techniques in order to derive analytical predictions, and extend these predictions numerically. Note that the neurite cytoplasm is in reality composed of several polymers that interact and contribute to global mechanical, geometrical and dynamical properties in a non-trivial way. Thus, continuum models are an approximation that must be used with proper caution, and that may be unsuited to representing subtle small-scale properties and processes. To capture molecular details and their emergent properties, many authors have recently examined neurites at a discrete level, for instance, by representing microtubules as discrete structures connected via cross-linking proteins ( . This sophisticated and realistic type of approach comes naturally with increased model complexity and computational difficulty, but allows for addressing finer questions on the effects of particular proteins or specific mechanisms.

Coupling mechanics and transport
Mechanics and transport are two coupled facets of neurite growth. The classic theory of transport-mediated growth assumes that neurite elongation is controlled by the concentration of material, supplied through a production rate that is generally taken to be constant. Clearly, this view cannot capture scenarios where axon length is controlled externally and not by tubulin concentration, like in stretch growth. In contrast, mechanical models describe the expansion rate as a rheological process akin to plasticity and involve the internal tension of the neurite as the main growth determinant. However, as a neurite elongates and thickens, it naturally increases its demand in material (O'Toole et al. 2008b). In particular, to sustain rapid elongation and adapt to changes in stretch rates, the cell body potentially actively modulates the production and transport of material to support growth and avoid disconnection or injury (Suter and Miller 2011;Athamneh and Suter 2015); however, little is known on this mechanism. For instance, Ahmed and Saif (2014) showed that increased tension correlates with faster vesicle transport, whereas Loverde et al. (2011) reported reduced fast axonal transport due to strain. Mechanisms underlying a stretch-mediated regulation of transport remain elusive and might involve mechanosensitive ion channels (Franze 2020); measures the ratio between the actomyosin viscosity and the effective microtubule core viscosity, and Q measures the excess of growth cone traction with respect to actomyosin contraction. We clearly see that growth is possible only if traction is sufficiently larger than contraction. Dashed line is an analytical estimate of the boundary between the collapse and motile states, obtained from a simplification of the model however, this has not been demonstrated systematically. Generally speaking, the rheology of neurites is a complex and dynamically regulated process that is not completely captured by passive mechanical models with constant kinetic properties. For example, it is plausible that parameters such as k in (19), or in (34) will depend on material availability, especially in long and fast-stretching neurites.
In contrast to the models detailed in Section 3.1.2 that focus on tip growth, O'Toole and O'Toole et al. (2008b) investigate the contribution of stretching in axonal transport of mitochondria. Here, we consider an isolated axon fixed at its base and stretched with speed v, with length (t) = 0 + vt . It can be experimentally observed that mitochondrial density is more or less uniform along the axon and increases with rate as time progresses. Thus, we assume a lineal density of the form c(x, t) = c 0 + t . From the balance of mass, we obtain the flux along the axon: where ln (2) × defines the half-life of a mitochondrion, and J 0 is the flux entering the axon at x = 0 . A no-flux boundary condition J( (t), t) = vc( (t), t) at the tip provides that grows quadratically in time. Along the axon, stretch results in a purely advective flux, which under uniform deformation is given by J stretch (x, t) = xvc(t)∕ (t) . Thus, the rest of the flux J other (x, t) = J(x, t) − J stretch (x, t) , due to cellular transport, decreases linearly along the axon as Figure 9 shows the spatial temporal evolution of both J stretch and J other . Note that the authors also address the case where the axon is attached to its substrate via adhesion (O'Toole et al. 2008a), in which case J other becomes the dominating contribution to the flux (since advection is localized at the tip, see Sect. 3.2.2).
By contrast to classic 1D models of transport (Sect. 3.1.2), the cellular contribution to the flux here decreases along the axon and depends on v and in a non-trivial way. This strongly suggests the existence of a mechanism allowing the axon to sense its own length and elongation rate, in order to achieve uniform mitochondrial density during stretch growth. However, this approach, which consists in weighing the various contributions to the flux using balance arguments, does not allow us to precisely infer the mechanisms involved.
To model mechanotransduction and growth in axons, Purohit and Smith (2016) reexamined Samuels et al.'s model to address the case of integrated axons subject to externally applied stretch, coupling mechanics to transport via a hypothesized stretch-activated ion entry (Fig. 10a). Here, the pulling speed is imposed, and therefore, transport processes mostly contribute to changing the internal state of the neurite, i.e., its cross-sectional area and tension. As in Samuels et al. (1996) and (13), the transport rate is assumed to be proportional to growth speed. To predict disconnection, an increase in tubulin production at the soma as a response to ion entry is hypothesized. The latter is promoted by axon tension, which decreases the free energy required for opening mechanosensitive ion channels along the shaft. The entry of ions in turn triggers an increase in cellular material production rate (S with our notation) which models the response of the axon to fast stretching. Finally, the resulting tension is compared with a critical disconnection tension, as shown in Fig. 10b, c, which allows to quantitatively evaluate the risk of nerve rupture, which is of interest, for instance, in the design of regeneration procedures.
In the last few decades, the active regulation of growth and cellular processes by mechanical forces has become a landmark of developmental biology (Hamant 2017). From a modeling point of view, the nonlinear coupling between forces and geometry can shed light on rich and non-intuitive behaviors. In neurons, this type of regulation remains poorly understood. Arguably, the development of mathematical models integrating active properties such as mechanosensing will be instrumental to building a unified theory of neurite growth.

Modeling the neurites in their environment
We have discussed models of neurites growing along a line with a focus on the fundamental mechanisms of elongation. However, the path and pattern formation of dendrites and axons growing in 2D or 3D is also of fundamental importance to understand the development of the brain architecture. Thus, we next describe the morphogenesis of neurons embedded in the multidimensional space and in relationship with their physical environment.

Guidance
During guidance, the axon tip perceives and transduces information from its local environment to find its path. Historically, our understanding of this process was chiefly predicated upon chemical signaling, in particular chemotaxis (Mortimer et al. 2008), which is the ubiquitous and predominant modality of growth cone pathfinding. For chemotactic information to be exploitable by a cell, a number of physical constraints must be simultaneously satisfied (see the classic article by Berg and Purcell (1977)). In a sequel of papers, Goodhill and coworkers developed seminal ideas to address the constraints of growth cone chemotaxis, using arguments from diffusion theory (Goodhill 1997;Goodhill and Urbach 2003;Goodhill 1998;Goodhill and Baier 1998;Goodhill and Urbach 1999). A core idea is that for a difference C in chemical concentration C to be detectable by a growth cone across its length scale r , several conditions must be simultaneously satisfied. First, the overall concentration must not be too high (saturated receptors) or too low (insufficient binding) (Goodhill 1998). Second, the gradient itself, e.g., the fractional difference p = C∕C ≈ ( r∕C) C∕ r , must be sufficiently large to overcome noise and provide an informative signal. The concentration C is obtained by solving the diffusion profile associated with a point source of ligand. It is well known that there is no physical (positive) solution to the diffusion equation in one or two dimensions, which can be alleviated by taking into account the natural extinction of the ligand (Krottje and Van Ooyen 2007). However, in three dimensions, a ligand diffusing from a steady point source has a concentration C given by with r the distance from the source; erfc the complementary error function; D the diffusion coefficient; and q the rate of ligand production at the source. Therefore, combining (43) with the constraints on detection provides the range of radii r ∈ r min (t), r max (t) within which detection is possible at a given time t, as shown in Fig. 11. While chemotaxis is the most studied and main modality of guidance, other less conspicuous guidance mechanisms also are at play. For example, there is now considerable evidence that the mechanical stiffness of the tissues is used as a regulatory cue during the formation of the line shows the limiting (t) curve above which the axon disconnects. Dashed lines show different stretching trajectories with piecewise constant growth rate. c Corresponding tension normalized by failure tension L . The black and blue curves illustrate two different disconnection scenarios associated with different traction parameters. Images reproduced from Purohit and Smith (2016) with permission from Elsevier nervous system (Abuwarda and Pathak 2020;Koch et al. 2012;Chan and Odde 2008;Koser et al. 2016;Franze 2020). Chan and Odde (2008) proposed a stochastic model of growth cone locomotion, modeling the actin motorclutch system at the edge of the growth cone, based on attachment and breakage of focal adhesion bonds depending on substrate stiffness (see also the mean-field treatment by Bangasser and Odde 2013). The authors predicted that the growth cone traction, which depends on its ability to pull its substrate, depends on medium compliance in a nonlinear fashion. On soft substrates, the growth cone develops traction by a so-called load-and-fail dynamics where the adhesion is maintained for a sufficient amount of time. However, on excessively stiff media the same mechanism results in a catastrophic loss of adhesion leading to a frictional slippage regime, in which the growth cone cannot progress efficiently. While these works discuss mechanisms and constraints at play in chemotaxis or mechanical guidance, they do not seek to describe the path of the growth cone as a response to these stimuli. To that end, trajectories of growth cones seen as a random motion on a 2D substrate have been studied by multiple authors (Katz et al. 1984;Mortimer et al. 2010;Maskery et al. 2004;Van Ooyen 1999, 2000;Krottje and Van Ooyen 2007;Borisyuk et al. 2008;Buettner et al. 1994;Wang et al. 2003;Pearson et al. 2011;Ben-Jacob 2000, 2001;Basso et al. 2019;Yurchenko et al. 2019Yurchenko et al. , 2021Davis et al. 2017;Roberts et al. 2014); see also the review by Maskery and Shinbrot (2005). The main focus in these works is not the mechanism of elongation, as treated before, but rather the shape and statistical properties of growth cone trajectories, as a function of external cues.
A simple model for the motion of a growth cone (t) consists in considering the competition between stochastic, random walk events, and deterministic events based on guidance cues. Within this framework, a general probabilistic process describing growth cone motion is where is a deterministic drift that depends on sensory cues; is the Wiener stochastic process; and D is a diffusion constant (Maskery and Shinbrot 2005). For example, the case where is constant (e.g., when the growth cone perceives a constant chemical cue) can be decomposed as (i) a translation t combined with (ii) a diffusion with typical travel length √ Dt . Note that (44) is not specific to growth cones and in particular does not capture the directional persistence of growing neurites (Katz 1985). An alternative to (44) is the Langevin equation for Brownian motion that includes inertial forces which can mimic this effect. Curvature effects can also be included by assuming that guidance cues affect the direction of the tip, rather than its position, which allows to dampen abrupt turns (Krottje and Van Ooyen 2007;Borisyuk et al. 2008;Pearson et al. 2011).
To model chemotaxis, Hentschel and Van Ooyen (1999), for instance, considered a version of (44) including multiple chemical agents , and assuming that growth cones follow the gradient of each concentration field c as where the coefficients quantify the strength of the response elicited by each agent , which can be either repulsive, < 0 , or attractive with > 0 . Mathematically, the chemotactic force assumes a conservative structure where the terms − c relate to chemotactic potentials.
These previous approaches are simple and relatively easy to compare with experiments. However, they are rudimentary from the point of view of cellular mechanisms. Several authors have represented the cellular effects involved in growth cone motion by using more detailed representations of the growth cone sensory and locomotory apparatus. For instance, Aletti et al. (2008); Goodhill et al. (2004); Xu et al. (2005) represent the growth cone as a semicircle covered
Parameter values are carefully discussed in the original article. For a given distance from the source, we compute the time interval during which detection is possible, represented by the green domain. The red domain represents the times when chemotaxis is not possible yet due to insufficient ligand concentration (minimum concentration C) or no longer possible due to a shallow gradient (constraint on p). For the parameter values used by the authors, the upper constraint on the concentration (saturated receptors) concerns points that are very close to the origin and is therefore ignored here. We see that for r ≲ 1500 μm , detection is always possible after a given time. Image adapted from Goodhill (1998) (with permission from Elsevier) with binding sites. These sites are engaged with a probability that depends on the local value of ligand concentration that may vary across the growth cone diameter. The new orientation of the growth cone is then modified according to the direction of maximum binding. An extended model accounting for locomotion was used by Betz et al. (2009) who also represented the random fluctuations of the growth cone edge governed by a Langevin-like dynamics controlled by ligand binding, from which they deduce growth cone trajectory. Haptotactic guidance has been addressed as well, for instance, by Van Veen and Van Pelt (1992), who modeled the contacts between the growth cone and discrete adhesive loci. When one adhesion site is detected within the growth cone, a force is generated and redirects the growth cone toward it. A similar approach was also employed by Li et al. (1995).
Aeschlimann and Tettoni made a step toward an integrated mechanistic approach of chemotaxis and proposed a feedback mechanism where a chemical gradient is first detected by the growth cone filopodia, which in turn triggers the entry of calcium ions at their base (Aeschlimann 2000;Aeschlimann and Tettoni 2001). The influx and diffusion of calcium along the growth cone periphery stimulate further protrusion of filopodia. (More complex computational models for calcium signaling have also been investigated, see Roccasalvo et al. 2015;Forbes et al. 2012;Sutherland et al. 2014;Hely et al. 1998.) The feedback loop proposed by Aeschlimann and Tettoni is combined with a mechanical model for filopodium retraction. This retraction produces a resultant that pulls the growth cone toward the stimulus. However, to actually operate a turn, the growth cone must bend the microtubule bundle (Franze 2020). To model this effect, the authors represent the trailing shaft as a chain of viscoelastic springs connected by angular springs. (This approach may be seen as the discretization of a morphoelastic rod, see Moulton et al. 2013.) In line with Aeschlimann and Tettoni's model of the morphoviscoelastic properties of the shaft, Zubler and Douglas (2009) proposed a detailed computational models of neurite migration in 3D, which includes steric hindrance due to other cells. Recently, García-Grajales et al. (2017) also proposed a cell-level morphoelastic approach by extending Recho et al. (2016) to a computational framework. In this approach, the authors use a multidimensional model of the shaft, including differential growth, contraction, bending, torsion or shear, and an advanced finite-element method is proposed for simulations over long distances, as shown in Fig. 12.
Overall, these models for axon growth and guidance describe, with various level of details, the behavior of a single axon in response to a stimulus field. They can be linked to studies of the growth cone to understand how a signal is sampled and integrated. Independently of these microscopic mechanisms, the overall macroscopic behavior of single axons is very much as expected. However, in many settings axons tend to interact and bundle together. Hence modeling these interactions is key to understand the formation of neuronal networks.

Collective migration and fasciculation
It is observed that axons often migrate and bundle together in a process called fasciculation which is believed to improve the accuracy of axonal projections toward their target (Van Vactor 1998). An initial fasciculation is often terminated by defasciculation to innervate the target zone and establish precise connections. Theoretically, this process of fasciculation and defasciculation is an example of collective behavior where different biological units come together and interact, creating in the process emergent structures and non-intuitive dynamics (Vicsek and Zafeiris 2012).
The Hentschel-Van Ooyen model given by equations (44-45) can be extended to include collective effects by considering growth-cone-secreted ligands, namely, chemical fields c i generated by multiple moving growth cones i (t) . Unsurprisingly, the net effect is the formation of densely packed bundles. Those can separate in response to a new regulatory cue secreted by the targets, and instructing defasciculation. Note that the authors ignored the transient effects due to diffusion of ligands from a moving source and focused on quasi-steady regimes. (Non-steady solutions were studied numerically by Krottje and Van Ooyen (2007).) In another version of the model, the authors introduce fasciculation due to short-range contacts combined with noise in the growth cone trajectories. Interestingly, this accounts for the emergence of so-called pioneer axons (Fig. 13), namely axons that grow first toward a target region, and that have been proposed to lead the way for the rest of the bundle. This behavior can be explained by the fact that, in the model, unbundled axons show more random motions and therefore explore more of their surroundings. Therefore, those axons are more likely to reach areas of the domain where Fig. 12 A simulation of axon growth and chemotactic guidance (the red fields denotes repellent and the green one the attractant). The color within the axon shaft gives the xx component of the stress within the structure, the x-axis being along the horizontal (based on the method given in García-Grajales et al. 2017) the chemoattractant gradient is steeper and chemoattraction is stronger, which leads them to further venture ahead of the group and migrate on their own. The track formed by pioneer growth cones is then followed by the rest of the bundle.
Analytical predictions are difficult to derive from these off-lattice multiagent models. By simplifying the geometry of the problem using a discrete random walk approach, Chaudhuri et al. (2011) examined the statistical properties of fasciculation. In their model, individual growth cones progress along the edges of a rhombic grid. At each discrete step of the automaton, each growth cone makes one step on the lattice and steers right or left according to a probability that depends on the presence of neighboring axons governing random fasciculation and defasciculation events. This model does not include chemotactic signals and focuses on axon-axon interactions. However, it includes other interesting processes such as the removal of existing axons with a given typical lifetime, or the difference in bundling affinities between two types of axon with different chemical signatures, resulting in sorting dynamics. This approach allows for the use of mean-field arguments, combined with relatively straightforward numerical simulations, to build insight into the non-equilibrium statistical mechanics aspects of fasciculation.
Another classic approach to collective cell migration and chemotaxis consists in assuming that a population of cells can be modeled as a continuous entity, which naturally lends itself to a continuum treatment in terms of partial differential equations (Murray 1993;Painter 2019). This approach has been employed to study the formation of neural crest (Giniūnaité et al. 2020), angiogenesis (Pillay et al. 2017) and epidermal wound healing (Sherratt and Murray 1990), for instance. Surprisingly, it has not been applied systematically to nerve cells, arguably due to their large and peculiar morphological structure that differs from most other cells. In a relatively unnoticed paper, the late Nobel Prize laureate P.-G. de Gennes proposed a simple mean-field model for collective axon migration (essentially in 1D), where a continuous population of growth cones is described by a density (De Gennes 2007). The flux J encompasses the isotropic random motion of growth cones, which is analogous to diffusion (with constant D 0 ), and a chemotactic flux J chemo (Murray 1993). In particular, the author postulates that a chemorepellent c is itself secreted by the growth cones. Upon fast ligand diffusion and decay, he assumes that c = C with C a constant. Thus, assuming that growth cones are repelled with velocity v = − c c∕ x as in (45), with c > 0 , one derives the advective flux J chemo = −C c ∕ x . The standard continuity equation then contains a diffusion parameter D = C c + D 0 that depends on . In contrast to regular diffusion, that would disperse the axons and preclude cohesive migration, this nonlinear equation supports solutions with a relatively sharp front (as illustrated in Fig. 14), suggesting a mechanism helping leading growth cones travel together. From a biological standpoint, the main assumptions of the model have not been directly justified and are greatly idealized. However, this general mathematical approach provides an interesting foundation for a mean-field treatment of axon migration.
Rather than focusing on the detailed physical interactions between individual axons, it is of interest to study the motion of an entire bundle of axons viewed as a single structure. For example, a bundle was recently modeled as a tip-growing morphoelastic filament in Oliveri et al. (2021). In particular, this model focuses on durotaxis (though it is also applicable to haptotaxis) and assumes that each individual axon has a migration speed that depends on local substrate stiffness, following Chan and Odde (2008); Koch et al. (2012). Based on  Oliveri et al. (2021), this effect is modeled by assuming that each growth cone composing the bundle's front produces a different force depending on local substrate stiffness, which varies across the bundle's finite width. Collectively, these forces create a torque that deflects the bundle. Surprisingly, this hypothesis leads to a behavior analogous to optic rays. For instance, at a sharp stiffness interface, a bundle is deflected according to a Snell-type law n 1 sin 1 = n 2 sin 2 linking the incident angles 1 and 2 via refractive indices n i that relate to medium compliance, as shown in Fig. 15a. Figure 15b also shows a simulation of the growing xenopus tadpole optic tract based on experimental brain stiffness data from Thompson et al. (2019). This simulation is in qualitative agreement with the observed trajectory, but does not include important factors such as chemorepulsion (Campbell et al. 2001;Piper et al. 2006;Atkinson-Leadbeater et al. 2010) or steric hindrance with other cell bodies (Koser et al. 2016).
Quantitative experimental validation remains necessary to estimate the relative importance of durotaxis, among the many other cues experienced by axon tracts in vivo.
Note that this model represents bundles of axons growing cohesively, as seen, for example, in the xenopus optic tract (Fig. 15), but not systematically in other types of nerve. Moreover, the emergence of the exact Snell law in the model is predicated upon simple assumptions on the bundle's mechanics (unshearable rod). In general, the bending properties of axons-a fortiori bundles of axons-remain poorly characterized. Nevertheless, although the model cannot capture guidance in fine details, its generality and simplicity might reveal essential properties of axon migration. In addition, the general result is potentially extendable to single axons that might follow similar optic-like behaviors at stiffness or adhesion interfaces. . We show two types of exact similarity solution, respectively, associated with the linear regime (red) where regular diffusion ( D 0 ) dominates, and in the idealized strong nonlinear regime (blue) where diffusion is small compared to growth-cone selfrepulsion ( C c ). In the second scenario, a sharp linear front of slope ∼ t −2∕3 advances slowly like ∼ t 1∕3 , whereas in the regular diffusion, the front flattens like ∼ exp Fig. 15 Optic-like reflection refraction of an axon bundle due to durotaxis as modeled in Oliveri et al. (2021). a Simulated trajectories of two bundles crossing a straight interface between a soft zone and a stiff zone. The curved portion of the trajectory represents the phase where the bundle's cross section (not represented) has only partially crossed the interface and is therefore subject to a torque. Outside of this zone of influence, the motion is rectilinear and the incident angles 1 and 2 are related via a Snell-type law, which accounts for both a reflection and a refraction regime that depend on the arrival angle. b Deflection of the xenopus tadpole optic tract in the middiencephalon where we can observe a stereotypical caudal turn of the tract toward the tectum. The tract originates in the right retina and then crosses the midline in the optic chiasm, to follow the contralateral brain surface toward the left optic tectum. Inset shows the trajectory of the axon bundle (blue) and the measured stiffness of the brain surface (heat map) obtained by atomic force microscopy (AFM) on the tadpole's brain surface (adapted from Thompson et al. 2019, under the terms of the CC BY 4.0 license). Note the contrast in stiffness between the front and rear parts of the brain surface, potentially responsible for the axon deflection (Koser et al. 2016). R.h.s. panel shows five representative simulated trajectories, with addition of noise and using a stiffness map obtained from AFM to define the stiffness field. The model, that includes only durotaxis as a guidance cue, reproduces qualitatively the trajectory of the actual optic tract, in particular the caudal turn toward the tectum that is represented by a quadrant in the top right corner

Branching and morphogenesis
As they develop, neuronal structures, axons and dendrites, may branch and turn into more or less complex tree-like structures. Different geometries and topologies of the neuritic trees are associated with different neuron types that have been linked with different neuronal activities (Masland 2004). A striking example is provided by cerebellar Purkinje cells, shown in Fig. 16a, that exhibit remarkably large and convoluted dendritic tree structures (Kapfhammer 2004). Understanding the processes that govern neurite branching and tree formation is important to build a picture of how the neuronal circuits and nervous functions are established. In particular, the failure of these processes has been associated with developmental abnormalities and cognitive impairment as observed, for example, in Down syndrome (Haydar and Reeves 2012;Mrak and Griffin 2004). Many generative models have been proposed to characterize and mimic the morphological diversity of neuritic trees (Ascoli 2002), using formalism such as Lindenmeyer grammars (Hamilton 1993;Ascoli 2002;Donohue et al. 2002) and other stochastic branching systems (Torben-Nielsen and De Schutter 2014;Carriquiry et al. 1991;Nowakowski et al. 1992;Uemura et al. 1995;Dityatev et al. 1995;Van Pelt et al. 1997;Van Pelt and Uylings 2005;Villacorta et al. 2007;Fujishima et al. 2012); or cellular automata (Luczak 2006;Albinet and Pelce 1996). A parsimonious example is the BESTL model which describes branching probability as a function of node depth (its centrifugal order) and the current number of terminal nodes in the tree (Van Pelt and Uylings 2005;Van Ooyen 2003). These approaches can reproduce observed neurons with remarkable fidelity and can be useful in applications such as classification. However, for most of them, the connection to intracellular physical mechanisms is not clearly made.
Mechanistically, neurite bifurcation occurs when the growth cone cytoskeleton splits to form two daughter branches. This event depends on several exogenous and endogenous factors (Bilimoria and Bonni 2013;Gibson and Ma 2011;Kalil et al. 2000;Kalil and Dent 2014) that have been modeled generally independently.
Several authors modeled branching as a process regulated intracellularly, and independently of the neurite's actual spatial embedding. Graham and Van Ooyen (2004) modeled the spatially dependent concentration of tubulin  with e a real exponent (Hillman 1979;Rall 1959), relating the area of the parent branch ( a p ) to those of its two children ( a c 1 and a c 2 ). Other authors have modeled the role of the MAPs (Kiddie et al. 2005;Hely et al. 2001;Graham and Van Ooyen 2004;Van Pelt et al. 2003) that are associated with increased neurite branching (Audesirk et al. 1997). Phosphorylated MAPs bind microtubules more loosely, allowing them to split. Kiddie et al. (2005); Hely et al. (2001) modeled MAP phosphorylation processes controlled by calcium in growing trees, using discrete multicompartmental approaches. For instance, Kiddie et al represented both MAP-dependent branching and tubulin assembly within the same framework.
Conversely, other authors have explored exogenous cues as the main mechanism for growth cone splitting. In this approach, it is assumed that branching results from conflicting forces applied to the growth cone (Van Veen and Van Pelt 1992; Li et al. 1995;Li and Qin 1996). Assuming that forces are applied at discrete points on the periphery of the growth cone parameterized by an angle, they define an angular distribution whose variance can be computed. In particular, a large variance signifies conflicting forces, e.g., a bimodal distribution. Branching occurs whenever this variance exceeds a given threshold, generally independently of past branching events, branch length, or current tree shape. These approaches are simple but remain phenomenological in the sense that the actual mechanical effects underlying branching are not explicitly represented and therefore are not connected to actual physical parameters such as the stiffness and viscosity of the growth cone, the adhesion with the substrate or the density of the cytoskeleton. Thus, there is an opportunity for a more detailed approach including a discussion of the actual mechanics at play in growth cone splitting.
Another remarkable phenomenon found in the morphology of neurons is the space-filling behavior of some cells, for instance, Purkinje cells (Fig. 16a). In these neurons, dendrites originating in the same cell body innervate the plane without self-overlapping (Fiala et al. 2008). This property probably facilitates efficient coverage of the space and prevents distinct branches to respond to the same input signal. Mathematically, this may be viewed as a problem of patterning with mutual avoidance, suggesting an approach based on reaction-diffusion equations. Similar to Hentschel and Fine (1996) who modeled dendritic growth by morphogen diffusion within the cytoplasm, Sugimura et al. (2007) proposed a compartmental activator-inhibitor model for dendritic self-avoidance, in which a hypothetical activator (u) is secreted and diffuses within the cell intracellular space. This activator then activates the production of an inhibitor (v), which conversely, suppresses the activator. The cell compartment is modeled by means of a density c(x, y, t) obeying the local growth law with a(u) < 1 , a decreasing function of u. A phase portrait analysis shows that growth occurs only when a(u) becomes negative, i.e., when the activator u becomes sufficiently concentrated; otherwise c = 0 is stable. We see that, in general, the only possible stable points are c = 0 and c = 1 , which allows us to outline the intracellular compartment c .
Contrary to classic reaction-diffusion, the activator diffuses and reacts only within c . By contrast, the inhibitor diffuses in the plane allowing for extracellular interaction between distinct branches. The general mechanism is illustrated in Fig. 16b. The two quantities u and v are coupled via Turing-like equations where f and g are polynomial reaction functions that satisfy basic conditions for the existence of patterns. Two simulations are shown in Fig. 16c, d. We see that the inhibitor produced by a branch diffuses in the medium and blocks the growth of neighboring branches by inhibiting the activator, thus precluding dendrite overlapping. Various parameters of the model account for qualitatively different types of dendritic trees, where the branching behavior (straight or wavy) relates to the different types of possible Turing patterns (dots or stripes). Despite its theoretical interest, this approach is relatively phenomenological and does not capture finer properties of the dendrites, such as variations in cross-sectional areas at branching loci (Hillman 1979;Tamori 1993;Rall 1959), growth inhibition between distant branches (as seen in Sect. 3.1.1), or other non-local and memory effects, such as the dependency of growth dynamics on the current topology and geometry of the tree. Furthermore, such activator-inhibitor pair, which provides a parsimonious and general mechanism for dendritic tiling, has not been identified yet in vivo, to our knowledge, and could likely emerge from a combination of multiple factors. Finally, contact-mediated rather than diffusion-mediated inhibition appears to be a dominant mechanism at play in space filling (Grueber and (47) Sagasti 2010;Matthews et al. 2007;Lefebvre et al. 2012;Soba et al. 2007).

Modeling the neurites in their environment
Our understanding of axonal growth and guidance is mostly based on controlled in vitro experiments where axons and dendrites evolve on a two-dimensional substrate. However, there is mounting evidence that the growth cone is a complex dynamic structure that can rapidly adapt its behavior depending on the chemical or mechanical properties of its environment. In particular, the extracellular matrix is a rich viscoelastic, three-dimensional environment that is very different from the typical substrates used in in vitro experiments, and where axon locomotion might be modified, as proposed by Santos et al. (2020). To understand how neurons develop, we must obtain a better picture of how neurites maneuver in complex environment, and how they interact with other cells during their development. Moreover, neurites receive multiple cues for guidance that are ultimately translated into mechanical forces for motion similar, in essence, to the problem of plant tropism (Moulton et al. 2020). Forces are probably key in orchestrating all these tropic responses. Indeed, there is mounting evidence that mechanics does not simply regulates the rate of neuronal growth, but also finely controls the nature and intensity of the response to chemical signals, suggesting that forces, mechanosensing and chemical processes are inextricably linked (see Franze 2020, and references therein). Challenge: Develop multiscale and multiphysics mathematical models to understand the relationship between signal transduction, cue integration, force generation and mechanical environment in a complex three-dimensional milieu.

Neuronal axon regeneration
Axonal injuries often result in loss of function due to the disconnection of neurons. Whereas in peripheral nerves, injury may be followed by regeneration and recovery of functions, in the central nervous system of mammals, rupture often leads to permanent disabilities as these neurons fail to regenerate and reconnect to their original target in order to rebuild a proper neuronal network (Mahar and Cavalli 2018). Strangely, the situation is quite different in fish and salamanders that have the capacity for long-distance axon regeneration and functional recovery after spinal cord injury. For instance, in the zebrafish, a spinal cord injury is followed by an immediate loss of function at the injury site, but in the long run, recovery is observed, associated with a complex response involving multiple cellular pathways and morphological changes (Tsata and Wehner 2021). Clearly, this is an extremely complex process, and there is an increasing understanding that mechanics has a role to play in it. Moreover, aside from the obvious fact that forces need to be generated for motion, elongation and reconnection, the mechanical environment is also affected as the surrounding tissue has been shown to stiffen during repair (Fig. 17), as evidenced by Möllmert et al. (2020). Hence, the authors argue that mechanosensing may be important in the regeneration process and it opens the tantalizing possibility that it can be manipulated to assist spinal cord repair.
Lastly, the extraordinarily high rate at which a single cell like a neuron extend during development remains a fascinating and open problem (Smith 2009). Understanding how an axon can sustain such extreme elongation without breaking and, conversely identifying the physical processes that limit stretch induced axon elongation will be crucial in the future of regeneration therapies.
Challenge: Develop mathematical models that include both biochemical and mechanical responses during axonal injury and neuronal regeneration to uncover the physical mechanisms contributing to successful spinal cord repair and nerve regeneration in general.

Fig. 17
Atomic force microscopy-based indentation measurements on acute zebrafish spinal cord sections showing the spatial distribution of apparent Young's moduli of an entire spinal cord cross section obtained from the level of the 12th vertebra. Image reproduced from Möllmert et al. (2020) with permission from the Biophysical Society

Growth of axons and tissues
In the brain, the white matter (inner layer of the cortex) is mostly made up of myelinated axon tracts that connect different areas of the outer cortex, rich in neuronal bodies and part of the gray matter (Goriely et al. 2015a). Therefore, the human brain embeds an intricate and highly anisotropic arrangement of fibers, illustrated in Fig. 18. A distinct morphological trait of large mammals is a folded-or gyrified-brain surface. Cortical folding-gyrification-is a crucial and widely studied aspect of brain morphogenesis. In particular, it has received a marked interest from the biomechanics and mathematics communities (Goriely et al. 2015a, b;Garcia et al. 2018;Greiner et al. 2021;Budday et al. 2014b).
Current experimental and theoretical evidence support the hypothesis that gyrification is primarily due to a mechanical instability. The latter is triggered by the relative difference between tangential growth rates in the gray and white matter, where the fast-growing cortex comes to buckle and fold to accommodate the slow-growing subjacent tissue (Ronan and Fletcher 2015;Goriely et al. 2015b). The potential contribution of axons in controlling this instability is possible but has remained elusive despite the popularity of the idea in neurosciences (Xu et al. 2010;Bayly et al. 2014). Most mechanical models (Ben Amar and Bordner 2017;Toro and Burnod 2005;Budday et al. 2014aBudday et al. , 2015Tallinen et al. 2016;Holland et al. 2015;Bayly et al. 2013;Wang et al. 2021;Budday et al. 2014b;Holland et al. 2018) are based on some version of morphoelasticity, which aims at representing the residual stresses caused by differential growth and therefore provides a natural paradigm for gyrogenesis. However, morphoelasticity is predicated upon a continuum formulation and therefore does not naturally take into account the microstructure such as axon tracts in any detail. The role of axons in gyrification has been addressed by Holland et al. (2015) who represented the preexistent pattern of nerve fibers as a field of directors defined at all points of the white matter. Based on the fact that axons grow under stretch and that they make up most of the white matter, the 3D growth tensor is assumed to be anisotropically biased to occur only along the fibers, all orthogonal stretches being purely elastic. The authors propose that axon tension, albeit not being the primary cause of folding, pilots the symmetry breaking and resulting folding pattern. However, in reality, the causal link is not clear as axon paths could be partially determined by the folding pattern itself. Modern imaging techniques and detailed connectome data (Fig. 18), in tandem with modeling, might prove valuable in studying the precise timing of these correlated events.
Challenge: Develop mathematical methods to explore the trajectories of migrating axons in the entire brain cortex, and the possible links between gyrification and fiber network geometry.

Conclusion
We have reviewed key mathematical models aimed at answering questions regarding neurite growth, pathfinding and patterning. Many of the fundamental mechanisms underlying these processes are only partially understood or remain elusive, and this is where mathematical modeling can play a key role. It allows for systematic testing of hypotheses and, with the emergence of quantitative biophysical models, can be used to identify fundamental mechanisms. In this regard, we want to emphasize the importance of mechanics in this process. Growth, morphogenesis, guidance all rely on the generation and sensing of forces that need to be appropriately taken into account. While simple zero-dimensional models are suitable to build intuition, we now understand that, from a mechanical point of view, these structures are not rigid but active morphoelastic materials that combine properties of fluids and solids together with active force generation, internal remodeling, and growth, requiring advanced mathematical models.
Many other models focus on transport processes or biochemical patterning. These different approaches offer snapshots of properties that must be combined in order to obtain a global picture. Indeed neurites are complex, dynamic and pluriform structures. Thus, future models must integrate the many rules that govern neurite protrusion, extension, sensing, collapse, stalling, steering, exploration, branching, fasciculation and defasciculation. The standard experimental Fig. 18 Diffusion spectrum imaging detects the movement of water molecules that flow along nerve fibers in the brain. The result is a map of the brain's neuronal network approach is to isolate different effects to understand a particular mechanism. By contrast, mathematical modeling, combined with numerical simulations and systematic validation, has the potential to integrate these effects in order to understand the emergence of complex morphologies and behaviors.
"Upstairs, my little noseminers! Go! Flee before me! Onward and upward! Go pump some neurons. Expand your craniums." Robin Williams, Mrs. Doubtfire, 1993. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.