Keywords: Galileo Galilei, Il Saggiatore, Rome, 1623 modeling developmental biology . . . . . . . . . .

intendere, se prima non il sapere a intender la lingua, e conoscer i caratteri ne quali è scritto. Egli è scritto in lingua matematica, e i caratteri son triangoli, cerchi ed altre figure geometriche, senza i quali mezzi è impossibile intenderne umanamente parola; senza questi è un aggirarsi vanamente per un oscuro labirinto.

How living objects acquire their shape has intrigued inquisitive minds for centuries and is still an open question (see [START_REF] Coen | The art of genes: How organisms make themselves[END_REF], for a layman's historical introduction to developmental biology). Incontrovertibly, this crucial aspect of biological development (known as morphogenesis) is for great part regulated by genes. Within a tissue, those may be expressed non-homogeneously, instructing each cell, more or less directly, with a specific behavior. For this reason, the study of genes and their translation into proteins is a cornerstone of developmental biology, which has seen increasing scientific and technological investments in the last decades.

However, living things are often more than the sum of their parts: knowing the individual behavior of each individual cell (or molecule) does not completely inform us on the whole body behavior. This remarkable feature of complex -in particular livingsystems (known as emergence) was well illustrated by various attempts to simulate life in silico, for example in cellular automata like Conway's Game of Life [START_REF] Berlekamp | Winning Ways for Your Mathematical Plays[END_REF]) or Langton's ant [START_REF] Langton | Studying artificial life with cellular automata1[END_REF]. Relying on minimal local rules governing their temporal evolution, these artificial systems may display rich and highly unpredictable large-scale behaviors.

Genetic processes undoubtedly play a central role in providing cells with properties (their local "rules"). However, they probably do not directly hard code the global behavior of a developing organism, which mostly emerges from more or less local cell-cell interactions. In particular the role of genes on the establishment of global shape is most plausibly indirect, and our understanding of it remains fragmented.

Introduction

This kind of nontrivial process, in general, is the object of the theory of morphoelasticity (exposed later on), that aims at mathematically characterizing the link between the shape and the mechanics of biological objects. This theory is in particular based on a two time scale resolution of growth and on a separation between the elastic and plastic components of deformation. In this thesis, I use a computational approach based on these general principles to model the force-based control of plant tissue growth.

Chapter 1

State of the art 1.1 Morphogenesis at the shoot apical meristem

Unlike animals, most plants continuously generate organs throughout their lifespan. In higher plants, the formation of those organs (shoots, roots, flowers, leaves, etc.), termed organogenesis, occurs at different precise zones of the plant named meristems.

The aboveground organs (e.g. leaves, flowers) in particular, are formed at the shoot apical meristem (SAM)1 (Fig. 1.1(a)). The SAM is a dome-shaped pool of stem cells, located at the tip of the plant axis [START_REF] Ha | Shoot apical meristem form and function[END_REF][START_REF] Steeves | Patterns in plant development[END_REF][START_REF] Traas | Cellular basis of shoot apical meristem development[END_REF]. It produces stem tissues (resulting in the axial expansion of the shoot) and initiates lateral organs (branches, leaves, flowers). Organ initiation occurs periodically, at the periphery of the meristem's central zone. As organs are initiated and expand, the meristem tip progresses due to the axial growth of the meristem, making room for next organs to appear (Fig. 1.1(b)).

Organogenesis at the SAM is spatially and temporally regular (following a given phyllotactic spatiotemporal program, Traas, 2013) and organ shape is highly conserved. This yet relies on the collective behavior of many cells that, individually, have only a local "knowledge" of the global picture. To solve this apparent paradox, a high level of regulation involving some type of cell-cell communication is necessary, guaranteeing a globally cohesive cell behavior. Supporting this view, several mechanisms have been identified as crucial regulators of morphogenesis -e.g. auxin transport (see [START_REF] Sassi | Auxin and self-organization at the shoot apical meristem[END_REF], for a review), architect genes [START_REF] Weigel | The ABCs of floral homeotic genes[END_REF] or mechanosensing [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] -without which the regularity of organogenesis (a) Cartoon of Arabidopsis thaliana's inflorescence (Image: www.pflanzenforschung.de) and close-up on the inflorescence meristem (scanning electron microscopy, scale bar = 100 µm, Image courtesy of J. Traas). (b) General structure of the shoot apical meristem (top view). Organs form at the periphery of the meristem at a regular pace, as the meristem tip progresses. (c) General structure of the animal and plant cells. Unlike animal cells, plant cells are embedded in a stiff wall. This cell wall prevents relative cell displacements, as two cells cannot detach from each other. By contrast, animal cells, that have their membranes in contact, may slide against one another more freely. Redrawn from [START_REF] Coen | The art of genes: How organisms make themselves[END_REF].

is dramatically impaired. To provide a rational, integrated view of how these mechanisms control morphogenesis in time and space, mathematical and computational models have been designed, in the form of virtual tissues coupling various biological features [START_REF] Chickarmane | Computational morphodynamics: a modeling framework to understand plant growth[END_REF]. In particular, one essential challenge in constructing these models is how to take growth into account.

In plants, cells are embedded within a stiff and sticky exoskeleton, the cell wall, that prevents them from moving relative to one another (Fig. 1.1(c)). In walled cells, osmosis results in a balance between cell turgor pressure and cell wall elastic stretch. Under sufficient turgor-induced tension, the cell wall may expand irreversibly, giving the cell room to absorb more water and augment its volume (see [START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF]Traas, 2010, for reviews). This irreversible process defines cell growth. In a multicellular context, due to cell-cell, short-range interactions, pressure forces not only affect cells autonomously, but also propagate across the tissue (see Fig. 1.2 and Boudon et al., 2015).

Local growth then results from the integration of both long and short range mechanical influences. To understand how those mechanical interactions integrate at the scale of an entire tissue, a number of biomechanical approaches have been proposed.. In the sequel, I critically review the various models published. I first present a kinematic description of plant growth, i.e. the set of variables that allow to describe the motion, independent of its causes, and second, an overview of the kinetics of growthi.e. the link between motions/equilibria and forces -that in plants has been described in various manners.

Kinematics of growth 1.2.1 Kinematic parameters of tissue motion

Meristem development excludes cell migration and apoptosis (i.e. programmed cell death) 2 . In this context, topological changes in the tissue result exclusively from cell divisions, and organ shape changes can only emerge from individual cell dilation. It was then proposed that complex morphogenesis can be achieved through spatial heterogeneities and anisotropies in the individual cell dilation speeds, i.e. through the fact that some cells (or groups of cells) grow faster than their neighbors and with a rate that depends on direction [START_REF] Coen | The genetics of geometry[END_REF].

Authors have stressed that, although overall shape deformation may be complex, it can be decomposed into a limited set of kinematic parameters (i.e. descriptors of the motion, independent of its causes), if taken at a sufficiently small scale [START_REF] Avery | Structure and development of the tobacco leaf[END_REF][START_REF] Coen | The genetics of geometry[END_REF][START_REF] Erickson | Modeling of plant growth[END_REF]. Those parameters are the main direction of the deformation, its degree of anisotropy, its volumetric expansion rate and its rotation rate [START_REF] Coen | The genetics of geometry[END_REF].

Mathematically, this can be formalized in the framework of continuum kinematics (see for instance Holzapfel, 2000, for an introduction). Assume that at a given time t, each point x of the tissue continuum moves relative to the lab frame with an instantaneous velocity v (x, t) (Fig. 1.3(a)). Clearly, v (x, t) gives little information on the growth at x, since it may encompass all kinds of motion, including, in particular, rigid transformations of the whole domain. A more meaningful descriptor of shape change would focus on the relative -rather than absolute -motion of the domain points. This relative motion is quantified by the velocity gradient (often noted L), that compares the instantaneous velocity at the control point x to that at neighboring positions:

L (x, t) def.
= ∂v ∂x (x, t) i.e. L ij = ∂v i ∂x j , i, j = 1, 2, 3 .

(1.1)

The velocity gradient L is a second order tensor field that describes the instantaneous linear transformation that affects an infinitesimal volume at the vicinity of point x.

In 1D, the velocity gradient describes the instantaneous lengthening/shortening of a given length element (see Fig. 1.3(b)). In general, this transformation may affect the size and shape of this volume, and its orientation (while differentiation eliminates uniform components, i.e. translations). This can be evidenced by splitting L into its symmetric (Σ) and skew-symmetric (Ω) parts:

L = Σ + Ω, with Σ def. = 1 2 L + L T and Ω def. = 1 2 L -L T , (1.2)
respectively termed Eulerian strain rate tensor and spin tensor (Fig. 1.3(c)).

As Σ is by construction a symmetric tensor (and positive since growth is an expansion process), it defines a direction of main extension (its first eigenspace), a dilation amplitude (its trace) and three degrees of dilation anisotropy (e.g. the pairwise ratios between the three eigenvalues). These anisotropies reduce to a single one, if we restrict ourselves to observe deformations in two dimensions. Note that tensor Σ vanishes in the case of a rigid motion.

As for the spin tensor Ω, it quantifies the rotational motion of particles and is dual to the instantaneous angular velocity vector ω:

Ω • h = ω × h where 2ω = curlv (1.3)
for any vector h (operator '×' depicts the vector cross product). Equation (1.3) makes the link between Ω and rotation speed more apparent. Further details may be found in [START_REF] Holzapfel | Nonlinear solid mechanics[END_REF] (see chapters Algebra of tensors and Rates of deformation tensors). Growth of a 1D rod. Between t and t + dt, the infinitesimal portion of length dx as moved and increased in length. The new length is obtained by substracting the position of its right extremity x + dx + v (x + dx) dt and that of its left extremity x + v (x) dt. The line dilation rate of the element is equal to the ratio of the new length to the previous length, or 1 + L (x) dt. (c) Schematic representations of the various kind of deformations applied to an infinitesimal square volume. Pure strain, pure rotation, and general deformation.

General principles of morphoelasticity

To model the physics of biological growth, authors stressed that part of the observed deformation is reversible and governed through elastic processes, and the other part is irreversible, governed through plastic processes. Mathematically, this idea was expressed in a continuum formalization originally proposed in [START_REF] Rodriguez | Stress-dependent finite growth in soft elastic tissues[END_REF] and generally referred as morphoelasticity (see also [START_REF] Garikipati | The kinematics of biological growth[END_REF]. This approach, summarized hereafter, was subsequently developed in [START_REF] Goriely | On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity[END_REF] which discusses a general procedure to represent the continuous growth of an elastic domain in the formalism proposed by Rodriguez and colleagues (see also [START_REF] Goriely | Morphoelasticity: a theory of elastic growth[END_REF]).

Let's assume that at initial time t 0 (resp. at time t ≥ t 0 ) the tissue is defined by an initial reference configuration B (t 0 ) (resp. current configuration B (t)). We consider the mapping x : X ∈ B (t 0 ) → B (t) that associates each material particle of the domain in the configuration B (t 0 ) with its position at t in B (t) (Fig. 1.4(a)). As in Section 1.2, the mapping transforms the infinitesimal volumes of B (t 0 ) according to a linear operation defined by the gradient of the mapping, that is termed deformation gradient and noted F :

F def. = ∂x ∂X .
(1.4)

Remark that, by virtue of the chain rule 3 , the previously-introduced velocity gradient L (Eq. (1.1)) can be related to the time derivative of F according to:

L = ∂F ∂t • F -1 .
(1.5)

The observed configuration B (t) generally contains mechanical stresses (residual or due to external forces). The fundamental postulate of morphoelasticity -formulated by Rodriguez and coworkers -is that, at any time, it exists a hidden stress-free state of the system, obtained by virtually suppressing external loading and "detaching" the tissue pieces apart from each other (Fig. 1.4(b)). By definition, growth is the change of this stress-free state with time, quantified by the growth tensor F g (the subscript 'g' stands for growth), that maps the initial (stress-free) configuration B (t 0 ) onto the current stressfree configuration B ′ (t). Note that there is no requirement for the deformation field F g to be compatible, i.e. to derive from a continuous displacement field. Indeed, as the control of growth is supposedly very local, the tissue pieces may readily grow discontinuously, hence no such compatible displacement exists in general. Compatibility is assured by means of the elastic deformation F e 4 (the subscript 'e' stands for elastic), at the price of introducing mechanical stresses. The total deformation is the composition of both deformations, which translates into the following multiplicative decomposition:

F B (t 0 ) → B (t) = F e B ′ (t) → B (t) • F g B (t 0 ) → B ′ (t)
.

(1.6)

We assume that only the elastic contribution F e affects the strain energy density Ψ of the system:

Ψ = Ψ (F e ) .
(1.7) 4 Again, F e does not define a compatible displacement if F g does not.

3 L = ∂v ∂x = ∂v ∂X • ∂X ∂x = ∂ ∂X ∂x ∂t • ∂X ∂x = ∂ ∂t ∂x ∂X • ∂X ∂x = ∂F ∂t • F -1 .
The variation of F g may be quantified by the following measure of the rate of growth:

L g def.

=

∂F g ∂t • F -1 g .

(1.8) that quantifies the growth-related part of the total velocity gradient (Figs. 1.4(b) and 1.4(c)).

Note that the decomposition expressed by Eq. (1.6) applies continuously in time [START_REF] Goriely | On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity[END_REF]. In the case of biological growth, the time scale associated with the elastic and viscous effects is generally several orders of magnitude smaller than that of the biological effects governing growth. Thereby, the transient viscoelastic regime is generally neglected and the system is considered in a state of quasi-static elastic equilibrium. Based on this assumption, Goriely and Ben Amar have discussed a general algorithmic method to compute growth, consisting of a succession of small incremental elastic relaxations, and small incremental irreversible deformations (see Goriely and Ben Amar, 2007, and Fig. 1.4(b)).

This general view is frequent in the computational modeling of plant growth, although originally, morphoelasticity describes the growth of elastic continua specifically (which does not apply to all models). We shall see that plant growth models often feature the most essential properties of morphoelastic materials, that are i/ the separation of irreversible and reversible deformations and ii/ the separation of elastic and plastic time scales. In an effort to uniformize concepts and notations, I will use the previous formalism as a conceptual backbone in what follows.

So far in this chapter, no requirement on the biophysics governing both F g and F e in time has been given. In particular, F g may depend on various quantities, like mechanical quantities (strain or stress), genes, temperature or nutrient concentration for example. In the sequel, I detail how different choices in implementing the elastic and plastic regimes have resulted in different models. I will mostly focus on morphoelastictype modeling approaches, from a very general standpoint. Accordingly, I chose not to comment agent-based models like J önsson et al. (2006), cellular Potts models like [START_REF] Grieneisen | Auxin transport is sufficient to generate a maximum and gradient guiding root growth[END_REF], or geometric (non mechanical) models like Holloway andHarrison (1999, 2007), [START_REF] Smith | A plausible model of phyllotaxis[END_REF]. A specificity of the morphoelastic approach is that the energetic state of the system is given by a displacement field, relative to a reference configuration, in contrast to viscous approaches (like [START_REF] Cheddadi | Coupling water fluxes with cell wall mechanics in a multicellular model of plant development[END_REF][START_REF] Dupuy | A biomechanical model for the study of plant morphogenesis: Coleocheate orbicularis, a 2d study species[END_REF], that express forces as a function of the rate of deformation. For the sake of simplicity and conciseness, some licence will be taken in the treatment of these cited models, as they will not systematically be the object of a separate review.

(a)

Observed config. at t Observed config. at t+dt Initial config.

Grown config. at t

Grown config. at t+dt Virtual configurations

(b)

Mechanical equilibrium

Rest configuration

Rest configuration Mechanical equilibrium (c)

Figure 1.4: (a) Cartoon of the mapping from initial configuration (at t 0 ) to current configuration (at t). (b) Schematic of the growth of a morphoelastic material. The current configuration is obtained as the combination of both a plastic deformation and an elastic deformation. Growth consists of a modification of the stress-free configuration of the system, that is virtual and does define a compatible displacement field in general.

(c) One-dimensional schematic of the growth of a linear spring. At time t the spring has a rest length l 0 F g (t) (where l 0 depicts the initial rest length), and an observed length l (t) = l 0 F (t) > l 0 (t) that results from the loading f (and depends on the spring constant k). The irreversible expansion of the spring is accounted for by an increase of its rest length, that results in a new equilibrium.

Kinetics of growth 1.3.1 Computing static equilibrium

The notion of mechanical equilibrium may be defined for a conservative system through the minimum potential energy principle (see [START_REF] Holzapfel | Nonlinear solid mechanics[END_REF][START_REF] Lanczos | The variational principles of mechanics[END_REF], for an introduction to the variational principles of mechanics). In principle, there exists a potential energy functional Π (F e ) for both stress and loads, associated with the system's configuration. The minimum potential energy principle dictates that an equilibrium of the system corresponds to a state for which the total potential Π (F e ) is stationary. If an explicit form of Π can be found, one can consider the problem of finding an equilibrium as an optimization problem.

In a number of cases however, this formalism cannot be employed. This is the case for example when the system is non-conservative, or when the external work depends on the movement of the system in a nontrivial manner (e.g. work of pressure forces when displacement is large). In such cases, a Newtonian, force-based formulation can be employed, consisting of solving a balance of force. In the computational biomechanical models of plants, both energy-based and force-based approaches have been employed.

Several authors have developed 2D force-based vertex models, representing the tissue in the form of a network of 1D mechanical elements such as linear springs [START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF][START_REF] Cheddadi | Coupling water fluxes with cell wall mechanics in a multicellular model of plant development[END_REF][START_REF] Rudge | A computational model of cellular morphogenesis in plants[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] possibly coupled with torsion springs [START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF]; or beam finite elements representing both tension and bending [START_REF] Dupuy | Coordination of plant cell division and expansion in a simple morphogenetic system[END_REF][START_REF] Dupuy | A biomechanical model for the study of plant morphogenesis: Coleocheate orbicularis, a 2d study species[END_REF]. In these models, each cell was represented as a polygon composed of discrete mass nodes connected through one-dimensional arcs (Fig. 1.5).

The equilibrium is obtained by balancing the pressure forces with the internal tissue forces. This was done either by integrating the second order equations of motion at each node, until all the node forces fall below a given threshold (for example in [START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF]; or through a first order approximation of the equation of motion, that assumes a strong viscous drag allowing to neglect inertia [START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF].

Energy-based formulations using similar descriptions of the tissue topology were also proposed; in 2D [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF][START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF] or 2D embedded in 3D [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. In [START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF], [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], the potential energy was expressed as the sum of the elastic energy (associated with edge deformation) and an external potential explicitly due to pressure. More abstract formulations were chosen in [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF], [START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF]. These rely on an approach often termed generalized Hamiltonian formulation, that allows to easily implement specific cell behaviors (e.g. anisotropic growth), e.g. by penalizing deviations from target cell shape [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF][START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF] or area [START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF].

Hamiltonian-based biomechanical modeling generally permits a broad variety of more or less arbitrary energetic contributions, sometimes at the expense of a clear connection with biophysical reality.

From a numerical point of view, Hamiltonian-based methods consider the optimization problem consisting of minimizing the total potential energy Π, that is explicitly formulated. One common heuristic is stochastic optimization, like Metropolis dynamics [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], that is used in [START_REF] Cieslak | Integrating physiology and architecture in models of fruit expansion[END_REF], [START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF]. In a nutshell, this algorithm explores the space of configurations (i.e. the possible moves) in a pseudo-random manner, that consists of looking for energetically more favorable moves, while occasionally authorizing slight increases of energy with a certain probability (in order to escape a non globally-optimal minimum). By contrast, the model detailed in [START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF] was based on a quasi-Newton minimization (the Broyden-Fletcher-Goldfarb-Shanno method, see [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] that is deterministic, and as the classical Newton-Raphson minimization method, consists of looking for a root of the energy gradient5 . To my knowledge, the authors of other cited works (i.e. [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] did not make mention of the optimization method they used.

In general, vertex-based methods like mass-spring models are relatively easy to implement, benefit from low computational costs, and provide tractable modeling frameworks to represent complex geometries and cell division. Nonetheless, they lack accuracy in approximating some specific properties of continuum multidimensional mechanics, such as shear or incompressibility [START_REF] Gelder | Approximate Simulation of Elastic Membranes by Triangulated Spring Meshes[END_REF]. In particular, the mechanics of in-plane walls is in general not represented, or at most, in an abstract manner. In addition, these models are restricted to 2D analysis (or 2D embedded in 3D in the only case of [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF].

To alleviate these limitations, higher-dimensional growth models have been developed.

In the line of energy-based approaches like [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF], [START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF], [START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF], Cieslak and coworkers have used a generalized Hamiltonian formulation in a 3D model, where the tissue is represented as a collection of connected compartments [START_REF] Cieslak | Integrating physiology and architecture in models of fruit expansion[END_REF]. In this model, the potential energy stored by each compartment is a function of its volume only. This is a coarse approximation, and this approach will fail in capturing the elasticity of the tissue more in general. For this purpose, a multidimensional formulation of the strain energy functional is required.

Figure 1.5: Cartoon of a mass-spring system modeling a multicellular tissue.

To provide a more proper physical description of the elastic regime, multidimensional continuum approaches have been proposed to describe tissue elasticity. These generally consider the strain energy density as a function of several variables, depending on multidimensional strain. The finite element method (FEM) has become overwhelmingly popular in simulating this kind of system (for an introduction to the general FEM, see for example [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF], or [START_REF] Zienkiewicz | The finite element method: solid mechanics[END_REF], for applications to solid mechanics). The FEM is a mathematically-sound, generic and efficient method, often used for structural mechanics, and implemented by a myriad of open source or proprietary packages 6 .

In modeling plant morphogenesis, authors have benefited from these convenient properties in modeling tissues in various forms and at different scales, yet with roughly comparable implementations. This was done for example in [START_REF] Fayant | Finite element model of polar growth in pollen tubes[END_REF], detailing a single-cell model where the cell wall is represented with 2D shell finite elements, which take into account both tangential stress and bending moment. Larger multicellular tissues have also been represented. Fozard and coworkers, for example have extended the previously cited vertex-based models by simulating the mechanics of in-plane walls using planar membrane finite element, that only take into account tangential strains [START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF]. Others have modeled the tissue in 3D, as an assembly of 2D meshed surfaces mimicking the network of cell walls forming the tissue, using 2D shell elements (Robinson andKuhlemeier, 2018, Sampathkumar et al., 2014a) or simpler membrane elements [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF][START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. Other models have considered larger-scale tissues seen as continuum domains (without cells), in 2D in [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF], [START_REF] Hong | Variable cell growth yields reproducible organ development through spatiotemporal averaging[END_REF] (that employ a model originally developed in [START_REF] Bonazzi | Symmetry breaking in spore germination relies on an interplay between polar cap stability and spore wall mechanics[END_REF], for fission yeast) and [START_REF] Coen | The genetics of geometry[END_REF], [START_REF] Kuchen | Generation of leaf shape through early patterns of growth and tissue polarity[END_REF], 2D embedded in 3D [START_REF] Kierzkowski | Elastic domains regulate growth and organogenesis in the plant shoot apical meristem[END_REF], or 3D [START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF].

A different approach was chosen by Bozorg and colleagues who used the triangular biquadratic springs (developed in [START_REF] Delingette | Triangular springs for modeling nonlinear membranes[END_REF]. Triangular biquadratic springs build a formal equivalence between spring-based modeling and membrane mechanics. This approach was used to model the tissue either as a pressurized membrane [START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF][START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF]; or as a multicellular complex [START_REF] Daher | Anisotropic growth is achieved through the additive mechanical effect of material anisotropy and elastic asymmetry[END_REF] similar to that modeled in [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF], [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], [START_REF] Mosca | On the micro-indentation of plant cells in a tissue context[END_REF], [START_REF] Robinson | Global compression reorients cortical microtubules in Arabidopsis hypocotyl epidermis and promotes growth[END_REF].

Generally, the cited authors have considered elements under finite (by contrast with infinitesimal) deformations, which is usually associated with a force-based resolution of the equilibrium, computed by solving the second [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF] or first order [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF][START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Bozorg et al., 2016) equation of motion, through various forward or backward integration procedures. At each integration step, external forces are computed from pressure and current geometry, while (internal) reaction forces are obtained as the negative gradient of the strain energy (itself often expressed based on Saint-Venant Kirchhoff-like hyperelastic formulations, see [START_REF] Holzapfel | Nonlinear solid mechanics[END_REF].

By contrast, in small deformations (e.g. in [START_REF] Hervieux | Mechanical shielding of rapidly growing cells buffers growth heterogeneity and contributes to organ shape reproducibility[END_REF][START_REF] Hong | Variable cell growth yields reproducible organ development through spatiotemporal averaging[END_REF][START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF], equilibrium can be directly found using a static FEM resolution, that consists of considering the equilibrium as the solution of a linear system, itself derived from a variational formulation [START_REF] Zienkiewicz | The finite element method: solid mechanics[END_REF]. When coupled with growth, this last solution generally requires a very small growth time step size to be accurate.

While more or less complex models have been proposed to express the stretch of the wall itself, the loading forces and their source (turgor pressure) were generally represented in a simpler manner. In the majority of the cited tissue models explicitly based on turgor (like [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF][START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF][START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF][START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF][START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Hervieux et al., 2016[START_REF] Hong | Variable cell growth yields reproducible organ development through spatiotemporal averaging[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF][START_REF] Robinson | Global compression reorients cortical microtubules in Arabidopsis hypocotyl epidermis and promotes growth[END_REF][START_REF] Rudge | A computational model of cellular morphogenesis in plants[END_REF], Sampathkumar et al., 2014a), steady and/or uniform cell turgor pressure was assumed, to simplify numerical analysis, and in absence of clear evidence of heterogeneities in cell pressures [START_REF] Long | Cellular heterogeneity in pressure and growth emerges from tissue topology and geometry[END_REF]. Typically, the cited authors have neglected the incompressibility of the liquid phase, which numerically would add a stiffening constraint that can be delicate to deal with. Very recently, authors have started reexamining the steady and uniform pressure hypothesis on both experimental [START_REF] Long | Cellular heterogeneity in pressure and growth emerges from tissue topology and geometry[END_REF] and theoretical [START_REF] Cheddadi | Coupling water fluxes with cell wall mechanics in a multicellular model of plant development[END_REF] grounds, using a 2D vertex-based model (and extending the discussion initiated by [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF][START_REF] Ortega | Augmented Growth Equation for Cell Wall Expansion[END_REF], to a multicellular context). This work is a novel step towards a multiphysical understanding of plant tissue mechanics, coupling both wall viscoelastoplasticity and water fluxes, and by doing so, allowing a more realistic treatment of pressure and incompressibility.

Computing growth

General implementation of growth

As soon as equilibrium is reached (see previous section), growth may be accounted for by modifying the rest configuration of the walls, via a modification of the growth tensor F g (Section 1.2.2). Due to Equations (1.6) and (1.7), this operation modifies the elastic deformation F e which generally results in a loss of equilibrium.

In one dimension, growth is defined as an increase of rest length as illustrated in Fig. 1.4(c). Such an approach has been used in several aforementioned vertex-based computational frameworks (i.e. [START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF][START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Rudge | A computational model of cellular morphogenesis in plants[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] where growth is represented as an increase of the rest length of cell edges. Note that the implementation of growth in [START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] does not completely comply with a morphoelastic description, as growth is not explicitly implemented as a modification of the rest lengths, but primarily as a displacement (rest length is redefined by the current lengths after displacement is applied, then elastic equilibrium due to pressure is computed). This geometric description is rather ad hoc and not based on physical considerations.

Other more abstract 2D and 3D Hamiltonian approaches have represented growth as an increase of cell (or tissue compartment) reference volume or shape [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Cieslak | Integrating physiology and architecture in models of fruit expansion[END_REF][START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF].

As an alternative to vertex-based models, other authors have considered the growth of a multidimensional material, in general in 2D (or 3D in [START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF]. To do so, Bassel and coworkers have extended 1D edge growth to 2D finite elements [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF]. Their solution consists of formulating growth as an increase of the rest length of the element edges, resulting in a redefinition of the element rest shape. Such a computational solution does not derive from an analytic continuum formulation and it precludes analytic treatments. In addition its physical and geometric interpretations are unobvious. This is because growth is here still seen as a 1D process, although the elastic regime is described through a proper multidimensional finite element formulation.

Bozorg and coworkers also pointed out the inadequacy of this approach in producing mesh-independent solutions and in properly representing mechanical anisotropies [START_REF] Bozorg | A continuous growth model for plant tissue[END_REF].

By contrast, the tensorial formulation of morphoelasticity (exposed in Section 1.2.2) provides a mathematically sounder, genuinely multidimensional extension of the 1D notion of rest length. This method has seen many applications to the growth of biological organisms (see Menzel and Kuhl, 2012, and references therein), and plant growth in particular, like in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], [START_REF] Bozorg | A continuous growth model for plant tissue[END_REF].

In this paragraph, I have discussed the procedures that are employed to define growth, from a strictly algorithmic point of view. In the sequel, I discuss the kinetic and physical assumptions that are made to govern the evolution of the rest configuration. These can be classified in two main views. In a first view, growth is assumed to be driven by morphogens (defined later). In another view, growth is seen as a mechanical relaxation process fueled by pressure-induced loading.

Growth as an explicit function of morphogens

How the rest configuration changes as time progresses depends upon the genetic regulatory network. To express this dependency, Coen and colleagues have proposed a parsimonious mechanism, suggesting that genes directly modulate the kinematic parameters of tissue growth [START_REF] Coen | The genetics of geometry[END_REF][START_REF] Green | Genetic Control of Organ Shape and Tissue Polarity[END_REF][START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF][START_REF] Kuchen | Generation of leaf shape through early patterns of growth and tissue polarity[END_REF]. This is done by expressing the growth tensor as a function of local gene expression and other abstract polarizing factors (that provide directional information):

L g = L g (genes, polarizing factors) .

(1.9) More generally, this approach consists in seeing growth as a direct function of morphogens7 , that are abstractly defined as any substance that locally, and more or less directly, controls growth.

Similarly, the growth of the element edges in [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF] was specified as a function of some morphogen concentration. Analogous solutions were used in [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF], [START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF], [START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] where growth was expressed as a function of abstract kinematic parameters, yet without explicitly invoking a morphogen.

The models cited in this section neither completely model the mechanical forces that cause tissue dilation, nor the rheological response of the tissue to these forces. The control of growth either relies on abstract biochemical parameters [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF][START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF] or phenomenological evolution laws [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Louveaux | Cell division plane orientation based on tensile stress in Arabidopsis thaliana[END_REF][START_REF] Merks | VirtualLeaf: an opensource framework for cell-based modeling of plant tissue growth and development[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF]. It was yet established, that plant cell expansion is due to pressure-induced forces, that cause the wall to expand, with a speed that depends upon cell rheological properties (see [START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF]Traas, 2010, andreferences therein).

Turgor-driven growth

To mechanistically bridge the gap between genes and shape change, another approach was proposed, where genetic and metabolic activities regulate cell expansion indirectly, by modulating the biomechanical parameters of the cells (e.g. pressure, stiffness, etc.).

In this approach, wall growth is explicitly seen as a response to forces, which are the product of turgor pressure.

In the approaches cited in Section 1.3.2.2, growth was decoupled from the elastic regime, and occurred independent of strain. This modality of growth can be physically seen as the result of some insertion of matter, that cause the tissue to dilate (which, in plants, does not occur independent of strain, [START_REF] Cosgrove | Growth of the plant cell wall[END_REF]. By contrast, turgorinduced growth is closer to a creep process, in which loading causes irreversible expansion provided a sufficient stretch is reached. At the scale of the cell wall, these two views are clearly in conflict. However, one cannot easily exclude the possibility that morphogen-driven growth at a continuous, supracellular scale may be a homogenization of pressure-driven growth at wall scale. For example, a morphogen in such an approach could be the local pressure itself.

More generally, these two views may in principle give rise to utterly different behaviors. In particular, their behavior with respect to residual stresses may differ, as proposed by Bozorg et al., who suggested that morphogen-driven growth theoretically favors large tensional and compressive residual stresses, compared to strain or stress-based laws [START_REF] Bozorg | A continuous growth model for plant tissue[END_REF]. Figure 1.6 illustrates this principle considering a 1D beam undergoing either strain-independent (Fig. 1.6(a)) or strain-dependent (Fig. 1.6(b)) growth. In this example, strain-based growth favors a decrease of the elastic strain and a dissipation of strain energy. By contrast, strain-independent growth induces a compression, possibly resulting in a mechanical instability (buckling). Such compression cannot occur in the first case, as plastic expansion will not occur in absence of tension. In (a), the growth may favor the accumulation of compressive stress, that eventually may result in mechanical instabilities. The beam may buckle under compression, reaching a more energetically-favorable state. This does not occur in case (b), as growth stops in absence of strain.

The accumulation of residual stresses was dodged -rather artificially -in several strainindependent approaches. In [START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF] the stresses ensuring tissue connectivity were canceled after the application of growth. In the 2D model detailed in [START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF], edge stretch was canceled after growth displacement was applied, leading to a new, stress-free configuration. This can be summarized by the following operation:

∀n > 0 : F n g ← F n-1 (1.10)
that, by virtue of Eq. (1.6), is equivalent to a cancellation of the elastic deformation (n refers to the growth iteration).

In [START_REF] Bonazzi | Symmetry breaking in spore germination relies on an interplay between polar cap stability and spore wall mechanics[END_REF], [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF], [START_REF] Hong | Variable cell growth yields reproducible organ development through spatiotemporal averaging[END_REF], growth was modeled through a similar stress cancellation, although the model used is these three works is explicitly strain-based. In fact, it is actually the simplest instance of a strain-based growth, which assumes fast and total relaxation of the strain at each growth step.

By construction, such a stress cancellation operation -which is justified neither theoretically nor experimentally -will fail in capturing the residual stresses that are experienced by the walls. In reality, as growth occurs nonhomogeneously and more or less fast, residual stresses may exist over long time periods, all the more that strainbased growth occurs only above a sufficient stretch level [START_REF] Cleland | Effect of osmotic concentration on auxin-action and on irreversible and reversible expansion of the avena coleoptile[END_REF][START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF][START_REF] Probine | Cell growth and the structure and mechanical properties of the wall in internodal cells of nitella opaca: Ii. mechanical properties of the walls[END_REF]. Below, the wall remains under elastic deformation. The plant tissue as a whole is therefore more plausibly in a homeostatic stressed configuration, as evidenced by cutting experiments [START_REF] Dumais | New evidence for the role of mechanical forces in the shoot apical meristem[END_REF] or using an adhesion-deficient mutant (Verger et al., 2018) for example. Models not relying on a cancellation of stress -typically morphoelastic approaches, in the case of continua -naturally provide more appropriate tools to model the possible effect of residual stresses on growth.

Several rheological models of strain-based growth have been developed in the form of plastic rheological laws. In particular, the fact that cell growth occurs only above a sufficient stretch level has been captured through various threshold-based formulations.

This approach was pioneered by Lockhart in the 60s, who proposed a simple piece-wise linear model of one-dimensional cell growth [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF]. This was later on completed by Ortega, who extended Lockhart's work to transient elastic regime [START_REF] Ortega | Augmented Growth Equation for Cell Wall Expansion[END_REF].

Lockhart's model describes the uniform elongation of a single turgid, cylindrical cell, resulting from the quasi-static balance between constant turgor pressure and wall stretch.

The formulation of wall creep is equivalent to that of a Bingham plastic, i.e. a plastic material yielding linearly above a threshold tension (see Fig. 1.7(a) and [START_REF] Bingham | An investigation of the laws of plastic flow[END_REF][START_REF] Steffe | Rheological methods in food process engineering[END_REF]:

L g = Φ σ -σ y + . (1.11)
Here, L g depicts the 1D rate of growth (Eq. (1.8)); Φ depicts the cell extensibility (the "rate of irreversible flow", in Lockhart's terminology), that can be interpreted as the inverse of the characteristic time of growth; σ depicts the wall stress (that, at steady regime, is proportional to turgor pressure); σ y is the yield stress, above which cell lengthening occurs; ( • ) + = max ( • , 0) is the ramp function.

Equation (1.11) provides a nonlinear rheological law capturing the two-regime cell behavior previously observed in vivo (Cleland, 1959, Probine and[START_REF] Probine | Cell growth and the structure and mechanical properties of the wall in internodal cells of nitella opaca: Ii. mechanical properties of the walls[END_REF]. Note that such a phenomenological abstraction disregards the molecular processes underlying growth. From a structural point of view, the wall material is an entanglement of various polysaccharides. In growing, meristematic or differentiating cells, the walls contain mainly cellulose, pectins and hemicelluloses that combined together form a more or less elastic solid medium [START_REF] Ali | Physical Models of Plant Development[END_REF][START_REF] Cosgrove | Growth of the plant cell wall[END_REF][START_REF] Wolf | Growth Control and Cell Wall Signaling in Plants[END_REF]. Several hundreds of proteins involved in wall modification have been identified. Many enzymes like pectin methyl-estherases or xyloglucan transgycosylase/hydrolases [START_REF] Armezzani | Transcriptional induction of cell wall remodelling genes is coupled to microtubule-driven growth isotropy at the shoot apex in arabidopsis[END_REF][START_REF] Cosgrove | Enzymes and Other Agents That Enhance Cell Wall Extensibility[END_REF], or proteins like expansins [START_REF] Cosgrove | Loosening of plant cell walls by expansins[END_REF], may disrupt the wall constituents, causing the wall to yield irreversibly. In growing cells, this yield is simultaneously compensated by synthesis and secretion of new components into the wall, maintaining its mechanical and structural integrity [START_REF] Cosgrove | Growth of the plant cell wall[END_REF].

To provide a mechanistic link between macroscopic rheology and microscopic wall biochemistry, several molecular-scale models have been proposed, mechanistically accounting for Lockhart-like behaviors (Ali andTraas, 2016, Dyson et al., 2012).

Essentially, Lockhart's model, which is one-dimensional, does not capture the various degrees of possible growth anisotropy. Mechanically, anisotropic growth may result from the fact i) that the forces applied to an infinitesimal wall parcel are anisotropically distributed and ii) that the mechanical response of the wall to these forces is itself anisotropic, due to anisotropies in the wall structural properties (Boudaoud, 2010, Hamant andTraas, 2010).

Several of the above-cited planar (or surface) models of growth have used viscous laws [START_REF] Dupuy | A biomechanical model for the study of plant morphogenesis: Coleocheate orbicularis, a 2d study species[END_REF][START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF], or Lockhart-like nonlinear laws [START_REF] Bessonov | Deformable cell model and its application to growth of plant meristem[END_REF][START_REF] Cheddadi | Coupling water fluxes with cell wall mechanics in a multicellular model of plant development[END_REF][START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Rudge | A computational model of cellular morphogenesis in plants[END_REF] to express the growth of the cell edges. In this kind of approach, growth anisotropy can in principle be achieved by i/ affecting different rheological parameters to each cell edge (as performed in Hamant et al., 2008, Rudge and[START_REF] Rudge | A computational model of cellular morphogenesis in plants[END_REF]; or ii/ by applying an anisotropic distribution of forces to the nodes.

Two-dimensional continuum modeling approaches, that more properly capture the elastic regime, have extended Lockhart-like approaches to multidimensional expansion, using a tensorial formalism [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Bozorg et al., 2016[START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF]. Using a mathematical expression reminiscent of Lockhart's formulation (Eq. (1.11)),

Boudon and coworkers have written the rate of growth L g as a piece-wise linear function of the (reversible) Green-Lagrangian strain E:

L g = Φ E -E y + , (1.12)
where E def.

= 1 2 F T e • F e -I ; E y def.

= E y I is a positive spherical tensor defining the strain threshold in all directions; ( • ) + depicts the tensor ramp function defined in Appendix 1.A and illustrated by Fig. 1.7(b) [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF]. This approach was extended in [START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF], where L g was a piece-wise linear function of strain, stress or morphogen concentration.

The growth law proposed in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF] (Eq. (1.12)) is formulated in terms of strain. Hence, it differs from Lockhart's original formulation, that was expressed in term of stress (Eq. (1.11)). This subtlety is not very meaningful in 1D models with linear constitutive relations, but becomes relevant considering anisotropic multidimensional materials like plant walls (see next Section 1.4), where strain and stress eigendirections do not coincide in general. The choice of a strain-based rather than stress-based growth law is motivated by a number of observations suggesting that the main orientation of the elastic strain correlates with that of growth, although, to my knowledge, this was never assessed in a quantitative manner. In organs like stems, for example, growth is almost exclusively longitudinal. Yet, due to the cylindrical shape of the stem, tangential stress is most probably mainly circumferential (from theoretical considerations, see Landau and Lifshitz, 1986), while, anisotropic stiffness properties supposedly promote main strain in the longitudinal direction, i.e. orthogonal to main stress.

As Φ is a scalar, Eq. (1.12) states that the direction of main growth and elastic strain coincide, and that anisotropic growth mostly results from anisotropic elastic strain. Elastic strain in turn depends on stiffness properties and mechanical loading. To model this dependency, Boudon and colleagues have chosen a linear coupling between the Green-Lagrangian strain and the second Piola-Kirchhoff stress -that describes the internal forces in the reference configuration (S = ∂Ψ ∂E ) -through the fourth order elasticity tensor (C w , the subscript stands for 'wall'):

E = C -1 w : S (1.13)
or equivalently S = C w : E (where ':' depicts the tensor double-dot product). Tensor C -1 w is often referred as the compliance tensor.

We see next how Eqs. (1.12) and (1.13) mechanistically account for anisotropic growth, and how its parameters can be controlled.

Origin and regulation of anisotropic growth

In this section we ask where anisotropic growth originates from and what levers are accessible to a single cell to regulate it.

We may first remark that the right-hand side of Eq. (1.12) displays several parameters that may potentially be affected by the cell. Rheological parameters such as Φ and E y emerge from structural and physiological properties of the cell. These can be more or less directly controlled through genes. Yet, as they are scalars, they mostly affect the growth rate, and offer very little control upon growth anisotropy. From this statement, and by virtue of Eq. (1.13), we understand that growth anisotropy mostly depends upon the possible anisotropy of tensor quantities like C w and/or S.

It is generally assumed that stress largely results from tissue shape [START_REF] Boudaoud | An introduction to the mechanics of morphogenesis for plant biologists[END_REF], which affects the distribution of pressure forces. This idea is based on the common assumption that the plant tissue globally behaves like a pressurized shell (see Beauzamy et al., 2015, and references therein), indicating that local stress and geometry/curvature are coupled through a Young-Laplace type relationship. This assumption largely prevails in the literature, as it provides a tractable conceptual framework for inferring mechanical stress, and a convenient simplifying assumption for mathematical and computational modeling (e.g. in [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Kierzkowski | Elastic domains regulate growth and organogenesis in the plant shoot apical meristem[END_REF]. The validity of this hypothesis however, was seldom assessed in a quantitative manner. Recently, computational modeling of multi-layered tissues led to suggest that the exact pattern of stress may in fact be not as simple as that expected from a simple shell-based consideration [START_REF] Ali | Simulations and analysis of turgorinduced stress patterns in multi-layered plant tissues[END_REF]. In addition, the assumption that shape prescribes stress, no longer holds as soon as one considers the residual stresses arising from differential growth rates. Most plausibly, these stresses may be dominant in zones of highly heterogeneous growth rate, like the crease between the central zone of the SAM, and fast-growing incipient organs [START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF].

In principle, each individual cell is able to regulate its osmolarity, and by this means, to modify its internal pressure and the stress experienced by its walls. Yet, this kind of possible regulation offers limited control on the anisotropy of the cell wall stress, as it will mostly affect the amplitude of stress (that can be quantified by the trace of S). In addition, as mentioned before, an important part of cell stress is due to forces applied by the rest of the tissue [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], which are transferred through cell adhesion (Verger et al., 2018). This leaves little room for a cell to autonomously regulate the anisotropy of its growth through stress.

By contrast, stiffness, which depends on the local structural and chemical properties of the wall, may be regulated by the cell in a more autonomous and direct manner.

In particular, the main load-bearing element of the cell wall is its network of cellulose microfibrils (CMFs), which are stiff strand-like polymers (for reviews, see for example [START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF]Traas, 2010). These may be polymerized in disorganized directions, conferring isotropic stiffness to the wall (Fig. 1.9(a)); or around a preferential axis, resulting in more or less anisotropic stiffness (Fig. 1.9(b)). As captured by Figure 1.8: In the pressure vessel assumption, the tissue is equivalent to a pressurized membrane, on which one can define a tangent stress field. The situation is as if the whole tissue behaved like a gigantic pressurized cell [START_REF] Kutschera | Regulation of cell expansion. The cytoskeletal basis of plant growth and form[END_REF]. Reprinted from [START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF].

Eqs. (1.12) and (1.13) and Fig. 1.9(c), this anisotropic stiffness may result in anisotropic elastic strain, and then anisotropic growth. CMFs are polymerized by transmembrane cellulose synthase complexes (CSCs) which move through the cell plasma membrane and synthesize cellulose in their wake [START_REF] Kumar | Plant cellulose synthesis: CESA proteins crossing kingdoms[END_REF]. In turn, this movement has been suggested to be guided by et al., 2012[START_REF] Heath | A unified hypothesis for the role of membrane bound enzyme complexes and microtubules in plant cell wall synthesis[END_REF][START_REF] Paredez | Visualization of Cellulose Synthase Demonstrates Functional Association with Microtubules[END_REF]. Microtubules are highly dynamic filamentous cytoplasmic polymers (Fig. 1.10(a)) that, collectively may come in disorganized arrangements, or in highly aligned bundles. This, and the functional association between CMTs and CSCs, suggest that CMT organization prescribes CMF organization. What directional cue regulates the coalignment of CMTs?

At first thought, it seems logical to propose that this regulation depends on local gene expression, through the control of some chemical concentration. As a concentration is a scalar quantity, potentially only a chemical gradient could serve as a directional cue.

Yet, to my knowledge, no scenario going in that direction was identified.

Alternatively, some intrinsic chemical polarity of cells would be a natural candidate for orienting CMTs. Several proteinse.g. PIN1 (Barbier de [START_REF] Barbier De Reuille | A protocol to analyse cellular dynamics during plant development[END_REF][START_REF] Heisler | An auxin-driven polarized transport model for phyllotaxis[END_REF][START_REF] Smith | A plausible model of phyllotaxis[END_REF] or BASL [START_REF] Mansfield | Ectopic basl reveals tissue cell polarity throughout leaf development in Arabidopsis thaliana[END_REF] -are asymmetrically distributed within a given cell, defining a polarity axis. At global scale, cells may also display cohesive polarity, making cell polarity a good candidate for an indicator of growth direction. Nonetheless, to date, no clear functional association was identified between these proteins and CMTs. In particular, while the alignment of CMT arrays may strongly correlate with PIN1 distribution, experimental results suggest the absence of a causal dependency between these, opening up the possibility of a third upstream factor [START_REF] Heisler | Alignment between pin1 polarity and microtubule orientation in the shoot apical meristem reveals a tight coupling between morphogenesis and auxin transport[END_REF].

This factor could be mechanical stress [START_REF] Heisler | Alignment between pin1 polarity and microtubule orientation in the shoot apical meristem reveals a tight coupling between morphogenesis and auxin transport[END_REF], Sampathkumar et al., 2014b), that in principle may provide a rich, multidimensional information field (three-component information: direction, anisotropy, amplitude). In vivo, CMT alignment strongly correlates with the theoretically predicted main axis of stress [START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Hejnowicz | Tensile tissue stress affects the orientation of cortical microtubules in the epidermis of sunflower hypocotyl[END_REF][START_REF] Robinson | Global compression reorients cortical microtubules in Arabidopsis hypocotyl epidermis and promotes growth[END_REF], Sampathkumar et al., 2014a[START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF][START_REF] Williamson | Alignment of Cortical Microtubules by Anisotropic Wall Stresses[END_REF]. This was evidenced at different scales, on normal tissuese.g. at the SAM epidermis [START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] or on pavement cell periclinal membranes (Sampathkumar et al., 2014a[START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] -as well as through more invasive experiments like cell ablation [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Sampathkumar et al., 2014a[START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF], Verger et al., 2018) or tissue compression [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Robinson | Global compression reorients cortical microtubules in Arabidopsis hypocotyl epidermis and promotes growth[END_REF][START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF].

If mechanical stress feeds back on CMTs, then in turn, as CMTs affect CMF distribution in the wall, stress indirectly instructs the local growth pattern and controls morphogenesis [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. 1.5 Understanding the emergent behavior of stress-adaptive growing tissues

Literature review

From an abstract point of view, this kind of control typically characterizes a nonlinear dynamical system in the form of a feedback loop (illustrated in Fig. 1.11). This so-called stress feedback hypothesis consists of a local rule governing cell behavior. In turn, such local behavior integrates at larger scale into a global emergent behavior. To date, this emergent behavior, that is highly nonintuitive, is not fully characterized.

The stress feedback hypothesis is still under debate as the precise mechanism remains elusive, and because experimental proofs are often indirect (see [START_REF] Hamant | Are microtubules tension sensors? Under review[END_REF].

Modeling approaches, integrating the stress feedback as a working hypothesis, may allow assessing its theoretical plausibility and implications on morphogenesis, notably by yielding testable predictions.

Several growth models have integrated a stress-based control of the anisotropy of stiffness. This was done through 2D spring-based modeling [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF], in 2D generalized Hamiltonian approaches [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF], as well as in continuum approaches [START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF], Hervieux et al., 2016).

Other authors have used a static model (without growth) to model the stress feedback in the elastic regime [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF]. In [START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF], CMTs were modeled in a discrete manner, by representing each CMT as a spring connecting two nodes of a given 2D cell polygon. Essentially, at each iteration, each node is considered as the potential starting point of a CMT, that is inserted only if a number of geometric conditions are verified. This method is relatively simple but rather ad hoc; in particular, the number of CMTs (therefore global stiffness) may depend on mesh resolution. In [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], where cells were modeled as a polygon made of springs, the spring constants were modulated by projecting the main stress axis (obtained as the circular mean of the node forces, see [START_REF] Mardia | Statistics of directional data[END_REF] on each edge. In the 2D Hamiltonian approach detailed in [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF], the rate of growth defined for each cell is a function of the cell stress tensor, hence growth does not explicitly result from a change in the cell rheological properties. Finally, in other continuum approaches such as [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF], [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF] stiffness coefficients directly depend upon various observables of the stress tensor (i.e. trace, anisotropy, main orientation).

In [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF], [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF], [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF], the link between stiffness and stress is modeled in a more or less phenomenological manner. These works in particular do not consider the dynamics and organization of underlying subcellular actors such as CMTs or CMFs. Therefore, the parameters (stiffness, sensitivity of the stress feedback) of these models are abstract, and not connected to the chemical constants that govern the reorganization of CMTs and CMFs. To advance towards a more quantitative understanding of plant growth, finer, multiscale approaches are required.

Agent-based models have also been developed, where microtubules are represented as individual agents [START_REF] Allard | Mechanisms of Self-Organization of Cortical Microtubules in Plants Revealed by Computational Simulations[END_REF][START_REF] Chakrabortty | A computational framework for cortical microtubule dynamics in realistically shaped plant cells[END_REF][START_REF] Dixit | Encounters between Dynamic Cortical Microtubules Promote Ordering of the Cortical Array through Angle-Dependent Modifications of Microtubule Behavior[END_REF], Mirabet et al., 2018[START_REF] Muratov | Mechanism of dynamic reorientation of cortical microtubules due to mechanical stress[END_REF][START_REF] Nédélec | Computer simulations reveal motor properties generating stable antiparallel microtubule interactions[END_REF][START_REF] Tindemans | Survival of the Aligned: Ordering of the Plant Cortical Microtubule Array[END_REF].

This kind of approach allows to finely study the dynamics of self-organization of microtubules, that rely mostly on mutual short-range microtubule-microtubule interactions.

However, as they focus on complex behaviors at molecular scale, involving an important number of microtubules, extending these models to virtual tissues composed of several hundreds of cells would be impracticable.

Objectives of the thesis

Continuum approaches provide a mathematically tractable framework that benefits from a huge arsenal of analytic and numerical techniques. However integrating small scale molecular features within such framework can be challenging. In this thesis I couple continuum mechanics with small scale molecular dynamics, proposing a parsimonious, coarse-grained model of the wall structure and stress-based CMT organization.

To some extent, this thesis was influenced by the reflections developed in [START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF],

where a pioneering approach of the problem is proposed, interrogating the macroscopic properties of the stress feedback and its potential role in the large context of morphogenesis. Along a similar line, viewing biological development as a nonlinear dynamical system, I here ask a number of general questions. What are the emergent properties of a growing material whose growth is itself controlled through the mechanical stress field? What stable shapes are allowed by this system? Under what conditions? Would growth be globally over-constrained if all walls stiffened along the main stress?

To answer these questions, I first developed a parsimonious, multiscale theoretical model of a stress-responsive tissue (Chapter 2). This model establishes a dependency between the chemical kinetics of CMTs/CMFs, and continuum tissue mechanics, in a more explicit manner than in previous approaches. Analytic and numerical tissuescale FEM studies were performed, in static regime (without growth). In this model, stiffness is a function of the CMF distribution, that itself results from CMT-guided cellulose deposition. CMT distribution in turn depends upon the stress that modulates their chemical stability, in a phenomenological manner (since no mechanism has been actually identified).

In Chapter 3, the model is coupled with growth, previously implemented in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF]. This model provides a computational framework in which growth behavior is connected to physical quantities and chemical constants. In particular, growth affects the density and orientation of the fibers (that are transported by the growing medium), and therefore feeds back on stiffness. This growth model is then applied to a practical biological problem in Chapter 4, suggesting that stress feedback may be a mechanism for the amplification and maintenance of incipient leaf flatness.

Appendices 1.A The tensor ramp function

Here, the scalar ramp function

( • ) + def.
= max ( • , 0) is extended to the set of the symmetric second order tensors of dimension d. Let T be such a tensor; T admits a unique rank decomposition, that reads:

T = d n=1 λ n t n ⊗ t n (1.14)
where {λ n } n∈[1,d] ⊂ R depicts the set of the eigenvalues of T , and

{t n } n∈[1,d] ⊂ R d is an
orthonormal family of corresponding eigenvectors. We define the tensor ramp function as:

(T ) + def. = d n=1 (λ n ) + t n ⊗ t n . (1.15)
This definition assures that (T ) + is a positive symmetric tensor, since all its eigenvalues are positive. As the first author of this publication, I am the main contributor to the work hereafter exposed, and wrote most part of the manuscript.

Abstract

A crucial question in developmental biology is how cell growth is coordinated in living tissue to generate complex and reproducible shapes. We address this issue here in plants, where stiff extracellular walls prevent cell migration and morphogenesis mostly results from growth driven by turgor pressure. How cells grow in response to pressure partly depends on the mechanical properties of their walls, which are generally heterogeneous, anisotropic and dynamic. The active control of these properties is therefore a cornerstone of plant morphogenesis.

Here, we focus on wall stiffness, which is under the control of both molecular and mechanical signaling. Indeed, in plant tissues, the balance between turgor and cell wall elasticity generates a tissue-wide stress field. Within cells, mechano-sensitive structures, such as cortical microtubules, adapt their behavior accordingly and locally influence cell wall remodeling dynamics. To fully apprehend the properties of this feedback loop, modeling approaches are indispensable. To that end, several modeling tools in the form of virtual tissues have been developed. However, these models often relate mechanical stress and cell wall stiffness in relatively abstract manners, where the molecular specificities of the various actors are not fully captured.

In this paper, we propose to refine this approach by including parsimonious biochemical and biomechanical properties of the main molecular actors involved. Through a coarse-grained approach and through finite element simulations, we study the role of stress-sensing microtubules on organ-scale mechanics.

Main quantities

k ′ φ (θ) = k ′ φ e -γf (θ) Angular microtubule depolymerization coefficient c 0 Total concentration of tubulin K ρ , K φ
Equilibrium constants of the CMF/CMT kinetics η Measure of the relative stiffness between the gel and the fiber

Introduction

In the absence of cell migration, plant development occurs through inhomogeneous and anisotropic irreversible deformations of the tissue, driven by turgor-induced forces.

It has been proposed that global shape changes are controlled locally, at the cellular level, through three geometric parameters: growth rate, growth direction, and the degree of growth anisotropy [START_REF] Coen | The genetics of geometry[END_REF][START_REF] Erickson | Modeling of plant growth[END_REF][START_REF] Kennaway | Generation of Diverse Biological Forms through Combinatorial Interactions between Tissue Polarity and Growth[END_REF].

How individual cells regulate these parameters is fundamental to understanding how plants control their shape.

Growth results from turgor-induced forces, that are counterbalanced by the cell wall, a stiff exoskeleton surrounding the cells and preventing them from bursting [START_REF] Hamant | The mechanics behind plant development[END_REF]. The wall is a complex assembly of cellulose microfibrils cross-linked by a network of polysaccharides, including hemicelluloses and pectins [START_REF] Cosgrove | Wall Structure and Wall Loosening. A Look Backwards and Forwards[END_REF].

When stretched above a certain threshold these cross-links may break, allowing microfibrils to slide (as modeled in [START_REF] Dyson | A model of crosslink kinetics in the expanding plant cell wall: Yield stress and enzyme action[END_REF] thus causing the irreversible expansion of the wall. New material is simultaneously added to maintain wall integrity [START_REF] Cosgrove | Growth of the plant cell wall[END_REF]. This irreversible process defines cell growth. Throughout development, cells may have different growth rates at different positions and times. In addition, since microfibrils can be deposited in more or less parallel arrays, and/or because the wall may be mechanically stressed unequally in various directions, cell growth may be anisotropic (Fig. 2.1(a), [START_REF] Baskin | Anisotropic expansion of the plant cell wall[END_REF].

The central role of wall mechanics in growth has motivated the development of a rheological formalization of tissue expansion, initiated by Lockhart on one-dimensional (1D) cell expansion [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF][START_REF] Ortega | Augmented Growth Equation for Cell Wall Expansion[END_REF]. In brief, Lockhart modeled the cell wall as a Bingham plastic undergoing linear irreversible expansion above a certain yield pressure. In recent years, computational modeling has been used to represent higherdimensional tissue expansion, in 2D or 3D (for a review, see [START_REF] Ali | Physical Models of Plant Development[END_REF]. Several studies were based on discrete tissue descriptions, where the tissue was modeled as a network of 1D elements [START_REF] Corson | Turning a plant tissue into a living cell froth through isotropic growth[END_REF][START_REF] Dupuy | A biomechanical model for the study of plant morphogenesis: Coleocheate orbicularis, a 2d study species[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. This allows a reasonable trade-off between modeling expressiveness and computational complexity, but fails at representing efficiently a number of specific properties of higher-dimensional mechanics [START_REF] Gelder | Approximate Simulation of Elastic Membranes by Triangulated Spring Meshes[END_REF], such as shear, incompressibility, and anisotropy. To alleviate this difficulty, other authors have adapted the formalism of continuous media to growth [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Bozorg et al., 2016[START_REF] Dumais | An anisotropicviscoplastic model of plant cell morphogenesis by tip growth[END_REF][START_REF] Dyson | A fibre-reinforced fluid model of anisotropic plant cell growth[END_REF][START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF][START_REF] Goriely | On the definition and modeling of incremental, cumulative, and continuous growth laws in morphoelasticity[END_REF], Hervieux et al., 2016[START_REF] Rojas | Chemically Mediated Mechanical Expansion of the Pollen Tube Cell Wall[END_REF][START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF]. Boudon et al. and Bozorg et al. developed a Lockhart-like 2D description of growth using a second-order growth tensor (L g ), expressed as a function of the elastic (reversible) deformation, represented by the elastic strain tensor (E):

L g = Φ (E -τ I) + (2.1)
where Φ is the wall extensibility, τ is a threshold strain, ( • ) + is the tensor ramp func- tion (defined in Appendix 1.A) that allows growth only in directions where deformation is higher than the threshold τ ; and I is the second-order identity tensor [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Bozorg et al., 2016). For the sake of simplicity, it is generally assumed that the cell wall extensibility (Φ) is quantified by a scalar coefficient or, equivalently, that growth and strain directions coincide. In the linear elasticity assumption (Hooke's law), strain and stress (depicted by tensor S) are linearly coupled via the fourth-order elasticity tensor (C w ):

E = C -1 w : S (2.2)
(':' is the tensor double-dot product). Combining the two previous expressions shows that anisotropic growth results from an anisotropy in stress and/or in stiffness. There is a reduced potential for cells embedded in a tissue, to autonomously regulate their local stress, which is mostly imposed by the rest of the tissue [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF]. By contrast, by modifying locally the wall's structure, cells can affect growth anisotropy in a more direct manner.

Such a local regulation does occur in plant tissues and has been proposed to rely partly on the activity of cortical microtubules. Microtubules participate in locally structuring the cell wall, by guiding the trajectories of the transmembrane complexes synthesizing the microfibrils (Fig. 2.1(b), [START_REF] Paredez | Visualization of Cellulose Synthase Demonstrates Functional Association with Microtubules[END_REF]. There are important indications, that microtubules in many cells align along the main stress direction [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Sampathkumar et al., 2014a,b, Williamson, 1990). This has led to the hypothesis that forces, generated in the tissues by turgor and differential growth rates, feed back into the structure of the cell wall via microtubule organization.

This so-called stress feedback has been integrated into tissue-scale models. Hamant et al. used a spring-based cell discretization of the tissue [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. In this approach, stress feedback is taken into account by modulating spring constants according to the angle between the spring and the cell's average loading axis. Other authors have instead used a continuous formalism, expressing a coupling between multidimensional stiffness and stress, in [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF], through triangular biquadratic springs [START_REF] Delingette | Triangular springs for modeling nonlinear membranes[END_REF] and in [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF], through the finite element method (FEM). In the last two implementations, stress is considered as an input-parameter that modulates the main axis of stiffness and its anisotropy.

These studies express an abstract phenomenological link between stress and stiffness, not explicitly taking into account the underlying molecular mechanisms. Here, we describe a modeling approach where subcellular processes involved in the wall synthesis are more explicitly described. Based on a parsimonious kinetic model of microtubule reorganization and microfibril synthesis, we propose a homogenized expression of the wall's dynamical elastic properties. We present a general expression of force-induced microtubule stabilization in the direction of higher tension. This simple model leads to subtle predictions on the dynamics of microtubule arrays and gives a quantitative description of wall stiffness. We derive an analytical, quantitative behavior of the microtubule cortex and the material mechanics. Tissue-scale FEM simulations of the system reproduced several macroscale behaviors previously observed in vivo.

Model

In this section, we introduce our model theoretically, by first detailing the elastic behavior of the wall considered as a fiber-reinforced medium, and then the kinetics of CMTs and CMFs.

Mechanical model

Mechanical equilibrium

The meristem is generally assimilated to a thin pressurized shell with uniform thickness (Fig. 2.2(a), [START_REF] Beauzamy | Mechanically, the Shoot Apical Meristem of Arabidopsis Behaves like a Shell Inflated by a Pressure of About 1 MPa[END_REF][START_REF] Kutschera | Regulation of cell expansion. The cytoskeletal basis of plant growth and form[END_REF]. This shell defines a continuum on which stiffness, strain and stress can be defined as tensor fields evolving in space and time, respectively C w (x, t), E (x, t) and S (x, t) (NB, since the next developments are made at a fixed position x, we get rid of the spatial parameter x in the following notations). Due to large radius of curvature (compared to thickness), we assume plane stress within the surface [START_REF] Landau | Theory of Elasticity[END_REF].

Because biological processes occur at a time scale much smaller than elastic relaxation, we assume quasi-static mechanical equilibrium at any time:

∀t ∇ • σ (t) + f ext (t) = 0 (2.3)
where ∇ • σ is the divergence of the tangential Cauchy stress and f ext is the density of pressure force applied normally to the surface at x.

Constitutive relation of the material

In higher plants, the wall is largely composed of three polysaccharides: pectins, hemicelluloses and cellulose [START_REF] Cosgrove | Wall Structure and Wall Loosening. A Look Backwards and Forwards[END_REF]. The former is polymerized as stiff, linear strands, the CMFs. They are embedded in a matrix of pectins and hemicelluloses that form an isotropic, soft hydrogel. In addition, other polysaccharides, e.g. hemicellulose, as well as numerous enzymes constantly interact and remodel these two main structural elements [START_REF] Wolf | Growth Control and Cell Wall Signaling in Plants[END_REF].

In the sequel, it is assumed that this fiber-reinforced hydrogel behaves as a linear deformable solid (Eq. (2.2)), composed of an isotropic homogeneous matrix and a fiber network (Fig. 2.2(b)). We assume that fibers are much stiffer than the rest of the material. For that, we neglect the mechanical effect of the possible interaction between the fibers and the matrix. The whole system (i.e. fibers + matrix) can thereby be seen as the parallel association of both components. The resulting elasticity tensor C w is therefore the sum of the fiber-related C f and the matrix C g :

C w = C g + C f . (2.4)
The gel is assimilated to an isotropic deformable solid, described by a Young's modulus and a Poisson's ratio. Tensor C g can be expressed in an arbitrary surface basis (e x , e y ), by employing the Voigt form1

C g = Y     1 ν 0 ν 1 0 0 0 1-ν 2     (2.5)
where ν and Y respectively stand for the Poisson's ratio and for the reduced Young's modulus2 .

Describing the fiber-related component, which contains the anisotropy of the whole material, requires a mechanical model of the fibers. Those are modeled by linear springs distributed in the surface's tangential plane. For the sake of simplicity, we assume that they do not mutually interact and dot not move relative to the medium3 .

Let ρ (θ) be the local angular density of fibers4 in an arbitrary tangential reference frame (e x , e y ) (Fig. 2.2(b)). Since mechanically, a fiber is apolar5 , function ρ is π-periodic. As one may expect, the elasticity tensor C f is a function of the density ρ. By assuming affine deformation of the fibers (namely that each individual fiber is deformed as the macroscopic medium), one shows that C f can be expressed according to the low-frequency Fourier coefficients ρ 0 , ρ 1 , ρ1 , ρ 2 and ρ2 (Appendix 2.A and [START_REF] Cox | The elasticity and strength of paper and other fibrous materials[END_REF]:

[C f ] = πY f ρ 0 16     3 + ρ 2 +4ρ 1 ρ 0 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 1 -ρ 2 ρ 0 3 + ρ 2 -4ρ 1 ρ 0 2ρ 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 2ρ 1 -ρ 2 ρ 0 1 -ρ 2 ρ 0     (2.6)
where Y f is the reduced Young's modulus of one fiber.

Assuming that the wall's thickness, the gel's stiffness and Y f are constant in time, Eqs. (2.4) to (2.6) imply that the evolution of the wall's mechanical properties only results from the dynamics of the CMFs distribution.

The dynamics of cortical microtubules and microfibrils

The orientation of cellulose deposition depends on the orientation of cortical microtubules that guide the trajectories of transmembrane cellulose synthase complexes (CSCs).

In our model, tubulin (namely, the elementary molecular sub-unit of CMTs) is described as a biphasic gel: it can be found either in a monomeric soluble state or in a polymerized form (the CMT). The distribution of CMTs is described through its πperiodic angular distribution φ.

In our model, CMT orientation influences only the orientation of cellulose deposition, and has no impact on the mass of cellulose polymerized per unit time. In brief, this relies on the hypothesis that a CSC protein participates to cellulose deposition regardless of whether it is bound to a CMT [START_REF] Emons | Microtubules and cellulose microfibrils: how intimate is their relationship?[END_REF], as detailed in Appendix 2.B. Hence, cellulose polymerization is expressed as a function of the normalized microtubule distribution φ * = φ/ π φ. In addition to that polymerization term, we introduce an isotropic decay that prevents unlimited accumulation of cellulose in the wall. At the considered time scale, we neglect growth, and then, the CMF advection that it would induce. Finally, we express the evolution of the angular density of CMFs as:

d t ρ (θ, t) = k ρ φ * (θ, t) -k ′ ρ ρ (θ, t) , (2.7)
(k ρ and k ′ ρ are kinetic constants).

The evolution of the microtubule angular density φ emerges from a description of the tubulin polymerization/depolymerization kinetics. Tubulin is present in the cytoplasm under two forms: a soluble one composed of free dimers (with concentration c free ), and a polymerized one depicted by φ. We assume that its total concentration c 0 (namely regardless of its form) is constant and uniform. This concentration verifies:

c 0 = c free + π φ. (2.8)
Microtubule dynamics is governed by the conservation equation of tubulin:

d t φ (θ, t) = k φ c free -k ′ φ (S, θ) φ (θ, t) .
(2.9)

The synthesis term (first term on the right-hand side) describes the isotropic CMT nucleation from the pool of free dimers (with kinetic constant k φ ). The second term models the microtubule depolymerization rate. Unlike in Eq. (2.7), the microtubule depolymerization rate varies with the orientation according to the angular function k ′ φ (S, θ), that models stress sensing. The last is described by assuming that microtubule stability increases with the tensile stress angular distribution f , reducing the degradation probability according to:

k ′ φ (S, θ) = k ′ φ e -γf (S,θ) , (2.10)
where k ′ φ stands for the inverse half-life of polymerized tubulin in absence of stress, γ is a (positive) coupling parameter, and f (S, θ) = e θ • S • e θ depicts the force per unit surface in direction θ (details on the establishment of Eq. (2.10) are provided in Appendix 2.C). This angular variability of the microtubule depolymerization probability over the local orientation of the stress field is the key feature of our modeling approach and accounts for the experimentally-observed tendency of microtubules to align in the direction of maximal tension [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. The general hypothesis behind this is that microtubules under tension undergo some conformational change making them less sensitive to katanin severing activity and therefore increasing their lifetime (Landrein and Hamant, 2013a, Sampathkumar et al., 2014a).

In the last two sections, we established a dynamical model of the wall's stiffness, by describing (i) the relation between the elasticity tensor C w and the microfibril bundle (Eqs. (2.4) to (2.6)), (ii) the microtubule-guided deposition of CMFs (Eq. (2.7)) and (iii) the stress-regulated kinetics of CMTs (Eqs. (2.8) to (2.10)). Importantly, Eqs. (2.4) to (2.7) and (2.11) can be conceptually combined into a single ordinary differential equation relating the time evolution of the stiffness elasticity C w to stress.

Tissue-scale numerical simulation procedure

Our model expresses the local behavior of the system, focusing on an infinitesimal portion of tissue in which stress is an exogenous parameter. We extended the analysis to a large-scale tissue, represented as an elastic surface, loaded with constant pressure (Fig. 2.2(a)). In this context, stress results from curvature and pressure; and may evolve in time as the stiffness properties vary.

Due to its nonuniform curvature, the tissue experiences nonuniform stress field. In practice, the stress and the elastic response of such a system (its displacement field) can generally only be determined by numerical approximation.

For that, we employed the FEM (in space) with triangle meshes, coupled with a finite difference time discretization. At initial time t = t 0 , the CMT/CMF system is considered at equilibrium (null-stress steady state, isotropic distribution, see Eq. (2.43)). The structures are not loaded and no residual stress is assumed.

Assuming that the CMTs' transient regime is much shorter than that of cellulose deposition, CMT distribution can be considered at equilibrium, which is described by the steady-state solution of Eq. (2.9):

φ ∞ (S, θ) = c 0 K φ e γf (S,θ) 1 + K φ π dθ ′ e γf (S,θ ′ ) ,
(2.11)

with K φ = k φ /k ′ φ , the equilibrium constant of the CMT stress-free polymerization process.

Under this assumption, only the CMF distribution can vary. It is iteratively updated, on each finite element, by considering the steady-state CMT distribution (Eq. (2.11)), and the CMF evolution equation (Eq. (2.7)) solved through Euler's method (Eq. (2.42)).

After this step, stiffness is updated according to Eq. (2.6), bringing the system out of equilibrium and providing an new initial state. For further details on the numerical pipeline, refer to Appendix 2.E.

Results

We detail hereafter the outputs of the model exposed in the previous section. We first expose how the CMT steady state accounts for the experienced stress field. We then show how the cell wall equivalent elasticity tensor is modified and how that leads to nontrivial strain response.

The steady CMT/CMF pattern reflects stress

This section exposes results on the steady-regime behavior derived from Eqs. (2.4) to (2.7) and (2.11) (refer to Appendix 2.D for detailed mathematical developments).

The steady-state expression of the CMT distribution φ ∞ (Eq. (2.11)) expresses their tendency to crowd in the direction of largest tensile stress (Fig. 2.3(a)). Quantification of the microtubule coalignment has been proposed by [START_REF] Boudaoud | FibrilTool, an ImageJ plug-in to quantify fibrillar structures in raw microscopy images[END_REF], based on the concept of nematic tensor (see Section 2.E.3). The authors define microtubule anisotropy as the ratio:

α ∞ φ = Q ∞ φ,1 -Q ∞ φ,2 Q ∞ φ,1 + Q ∞ φ,2 ∈ [0, 1] (2.12) where Q ∞ φ,1 ≥ Q ∞ φ,2
≥ 0 are the eigenvalues of the CMT nematic tensor. This definition of anisotropy can be expressed in the Fourier formalism, as it is equivalent to:

α ∞ φ = | φ∞ 1 | φ∞ 0   = φ 2 1 + φ2 1 φ 0   (2.13)
where { φ∞ n } n∈N depicts the Fourier sequence of φ ∞ (see Section 2.E.3 for a proof). Similarly, the anisotropy of stress can be defined as α f = (S 1 -S 2 ) / (S 1 + S 2 ) (where S 1 ≥ S 2 are the stress eigenvalues)6 . Again, by developing the expression f (S, θ) =

e θ • S • e θ , one shows that α f = 2| f1 |/ f0 (= 2 f 2 1 + f 2 1 /f 0 ).
As shown in Figs. 2.3(a) and 2.3(b) and Eq. (2.35), α ∞ φ is a growing function of the stress anisotropy, but also of the mean stress f 0 = Tr S/2. Similarly, the concentration of CMTs (πφ ∞ 0 /2) is a growing function of stress anisotropy and amplitude (Fig. 2.3(c) and Eq. (2.33)). Note that if high stress amplitude favors high CMT coalignment, we can expect that this effect also comes with a lengthening of the transient regime's duration, due to the augmented stability of CMTs that diminishes CMT turnover7 .

Considering infinitesimal displacements, one can assume that the global shape, and therefore the distribution of the pressure forces, do not vary during the evolution of the material stiffness. This allows to derive the steady-state CMF distribution from Eq. (2.11): ,θ) π dθ ′ e γf (S,θ ′ ) .

ρ ∞ (S, θ) = K ρ φ * ∞ (S, θ) = K ρ e γf (S
(2.14)

The anisotropy of the CMF steady distribution, defined as

α ∞ ρ = |ρ ∞ 1 |/ρ ∞ 0 , is equal to that of CMTs (α ∞ φ )
, but its concentration does not depend on stress (due to the normalization of φ in Eq. (2.7)). The CMF steady distribution is also independent of the CMT chemical kinetics.

The effect of stress feedback on elastic deformation

Under the previous assumption that stress does not vary much as stiffness evolves in time, the steady-regime elasticity tensor C ∞ w is an explicit function of the CMT distribution (Eqs. (2.4) to (2.7), (2.14) and (2.33)). In particular, one can choose (without loss of generality) the surface basis such that e x coincides with the main eigenvector of stress.

In this basis, the CMF distribution ρ ∞ is even (its odd Fourier coefficients φ∞

c 0 = 1 mol.µm -3 , γ = 10MPa -1 , K φ = 1. vanish) and C ∞ f becomes: [C ∞ f ] = K ρ Y f 8     3 + 4α ∞ φ + β ∞ φ 1 -β ∞ φ 0 1 -β ∞ φ 3 -4α ∞ φ + β ∞ φ 0 0 0 1 -β ∞ φ     (2.15)
where

β ∞ φ = | φ∞ 2 |/ φ∞ 0 and K ρ = k ρ /k ′ ρ
is the isotropic-stress equilibrium constant of cellulose deposition. By combining the previous equation with the expression of the Fourier coefficients of φ (Eq. (2.33)), tensor C ∞ w = C ∞ f + C g can be written as a function of stress, yielding the stress/strain nonlinear relationship illustrated in Fig. 2.4(a).

Eq. (2.15) states that the material is asymptotically more resistant to traction in the direc-

tion of main stress (C ∞ w,11 ≥ C ∞ w,22
). Moreover, it is orthotropic (symmetric) according to the stress axes (C ∞ w,13 = C ∞ w,23 = 0)8 . The expression of stiffness that was obtained is specific to fibrous materials and provides for that more realistic mechanical behavior. This is due to the asymmetry of the CMF distribution, due to a transient phase shift between the respective distributions of CMTs and CMFs, representing the inertia of the system.

In anisotropic elasticity, the respective directions of maximal stress and strain are in general different. In our specific case though, at steady-state, the orthotropy of the elasticity tensor in the stress eigenbasis imposes that both fields (stress and strain) share the same eigenvectors. An alternative remains: either directions of maximal (and minimal) eigenvalues of both fields are parallel or orthogonal. Those two regimes were notably characterized in [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF] according to the anisotropy of the material's stiffness. We show how the orthogonality between the main directions of stress and strain can be expressed as a condition on stress only, using the dependence of stiffness on stress. By comparing the ratio

E ∞ ⊥ E ∞ = S 2 C ∞ w,11 -S 1 C ∞ w,12 S 1 C ∞ w,22 -S 2 C ∞ w,12
(2.16) (obtained by inverting Eq. (2.2), and illustrated in Fig. 2.4(c)) between longitudinal strain E ∞ (collinear to main stress) and orthogonal strain E ∞ ⊥ (orthogonal to main stress), to 1, a necessary and sufficient condition for orthogonality can be derived:

α f < C ∞ w,11 -C ∞ w,22 C ∞ w,11 + C ∞ w,22 + 2C ∞ w,12
.

(2.17)

By expressing matrix coefficients C ∞ w,ij according to Eqs. (2.4), (2.5) and (2.15), the previous condition becomes:

α f < α ∞ φ 1 + 2η with η = (1 + ν) Y K ρ Y f , (2.18)
where η measures the relative stiffness between the gel and the fiber components; the quantity α ∞ φ (1 + 2η) -1 is a measure of the material's anisotropy. This strain orthogonalization directly relies on the sharpness of the CMT alignment, and therefore on the stress field (Eqs. (2.11) and (2.35) and Fig. 2.4(c)). Note that since the anisotropy ratio α ∞ φ lies in [0, 1], strain orthogonalization cannot occur for excessive stress anisotropy. In fact, Eq. (2.18) cannot be verified if α f ≥ (1 + 2η) -1 ; this limitation comes from the residual isotropy due to the isotropic component (Eq. (2.4)). In brief, we have shown how two mechanical regimes (orthogonality or collinearity of the stress and strain main directions) can be explicitly accounted for by the chemical properties of the CMT array (notably through the sensitivity parameter γ). In particular, strain orthogonalization cannot occur if microtubules are less anisotropic than stress (Eq. (2.18)).

Although the number of CMFs, regardless of their orientation, does not depend on stress, we observe differences in the material's areal dilation between the stress-responsive (γ > 0) and isotropic (γ = 0) materials. This difference is a general property of fibrous materials, for which areal dilation varies with material anisotropy. Fig. 2.4(b) illustrates the emergence of this effect from stress directly. It allows identifying a particular domain where stress feedback produces higher dilation than in the case of an isotropic, non stress-sensing material. In a biological context, that implies that anisotropic-walled cells could expand faster than isotropic ones, all other parameters being equal. In addition, due to the non-linearity in the strain-based growth equation (Eq. (2.1)), an anisotropic cell may experience one-directional growth, whilst its isotropic counterpart could be kept under the strain threshold in all directions, preventing its irreversible expansion.

Effect of strain orthogonalization on tissue-scale development

In order to characterize the impact of the feedback process on an entire structure, we compare the strain field of an inert (non-responsive) material (γ = 0) and a stressresponsive material (γ > 0), at steady regime. In fact, Eq. (2.1) formulated in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF] expresses the local growth rate as a piece-wise linear function of the strain, suggesting that the elastic response of the tissue directly allows computing its immediate growth. In our context, where no growth is considered, we infer qualitative properties of growth from elastic strain.

One of the simplest shapes, involving anisotropic stresses, is the cylinder. Similar geometries have been studied in the context of plant modeling as they can be considered as idealized stems [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF][START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF]. We employ here a cylinder closed by a hemispherical cap and fixed at its boundary (Fig. 2.5). In this case, stress is known from continuous mechanics: it is anisotropic in the cylindrical trunk where hoop stress S 1 equals twice axial stress S 2 , and isotropic on the cap (mean stress is equal to S 2 ), as confirmed by FEM (Figs. 2.5(a) and 2.5(b)).

In absence of feedback (γ = 0, then CMT distribution is isotropic, Fig. 2.5(c)), circumferential stress results in circumferential strain (Fig. 2.5(e)) that should tend to break the cylindrical shape. By activating the feedback process (γ > 0), CMTs align in a circumferential way on the trunk, and remain isotropic on the cap since stress provides no directional bias here (Fig. 2.5(d)). In absence of feedback, dilation is a linear function of f 0 , which is not the case with feedback. (b) Relative variation (∆ ON -∆ OFF ) /∆ OFF (in %) of areal dilation between a stress-responsive material (γ = 1) and a non-responsive material (γ = 0), as a function of f 0 (expressed in MPa) and stress anisotropy α f . Blue bold contour indicates level 0 (∆ OFF = ∆ ON ). (c) Ratio E ∞ ⊥ /E ∞ of strain (in the stress eigenbasis, where the basis vectors are sorted according to the descending order of eigenvalues) as a function of f 0 (expressed in MPa) and α f , with stress feedback. Main stress and strain are orthogonal iff E ∞ ⊥ /E ∞ > 1; domains labeled ⊥ and respectively indicate those where main stress and strain are orthogonal or collinear. Blue and red bold contours respectively indicate levels 0 (null orthogonal strain E ∞ ⊥ ) and 1 (isotropic strain

E ∞ ⊥ = E ∞ ). Parameters: c 0 (arbitrary), γ = 10 MPa -1 , K φ = 1 (arbitrary), K ρ = 1 mol.µm -3 , Y f = 2 GPa, ν = 0.2, Y = 0.2 GPa.
Enabling the feedback (γ > 0) may allow strain to become maximal in the axial direction (orthogonal to stress), provided a value of γ satisfying Eq. (2.18) (Figs. 2.5(b) and 2.5(f)). This results in a slimming and an axial stretching of the cylindrical trunk.

One also observes that the dilation of the finite elements (areal expansion, equal to Tr E for small deformation) is bigger with the feedback than without (Figs. 2.5(e) and 2.5(f)), corresponding to the positive domain in Fig. 2.4(b). The previous simulation involves a positively curved shape. Yet, plants, and in particular the meristem, display several negatively curved zones (for example at the frontier between organs). This scenario has been idealized by merging two sphere meshes, creating a negatively curved junction (Fig. 2.6).

In that case, the outgrowth generates a high circumferential stress at the concave bend between the two spheres, where the gaussian curvature is negative (Figs. 2.6(b) and 2.6(c)).

Consequently, the strain field is also highly circumferential in absence of feedback (Fig. 2.6(f)). This contributes to augment the size of the frontier, relative to the reference state.

In contrast, the stress feedback allows to compensate the circumferential stress, and develop axial strains (Figs. 2.6(c) and 2.6(g)). In this case, the decrease of the circumferential strain contributes to maintaining the size of the frontier, reinforcing the separation between the meristem and its primordia, for instance. This example reinforces the idea exposed in [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF], that the feedback promotes the anisotropy of shape, preserving non-spherical curvature.

As a supplemental layer of shape complexity, the effect of a topological change has been examined. We represented this scenario through a circular hole in a hemispherical cap (fixed at its boundary), mimicking the ablation experiments described in [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Sampathkumar et al. (2014a).

Here, highly anisotropic stresses are generated at the hole's vicinity (Figs. 2.7(b) and 2.7(c)).

In absence of feedback, this would promote an expansion of the hole's size (Fig. 2.7(f)).

With feedback, the steady-regime strain features radial orientation (Fig. 2.7(g)) or, at least, less anisotropic circumferential orientation (where strain orthogonalization fails, see Eq. (2.18)), allowing a decrease of the hole's circumference relative to the nonstress-sensing case.This comes with a relaxation of the circumferential stress amplitude (around 10% with the present parameters) at the close periphery of the hole (Figs. 2.7(b) and 2.7(c)). To sum this up, it seems that a strain-based growth, coupled with a stress feedback may tend to close topological defects such as holes.

Discussion

Existing modeling efforts assessing the feedback mechanism between mechanics and growth in plant tissues rely on a phenomenological black box to express the influence of stress on cell wall stiffness [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Hervieux et al., 2016). This work is an effort to open this black box by making a first step toward a more quantitative, mechanistic, multiscale analysis of stress-based feedback. This approach

, K φ = 1, K ρ = 1 mol.µm -3 , Y f = 20 GPa, ν = 0.2, Y = 0.3 GPa.
relies on a parsimonious description of the chemical behavior of cortical microtubules, coupled to cellulose microfibril deposition. A homogenization procedure allows us to derive from the angular distribution of these microfibrils a continuous expression of stiffness as a function of its low-frequency Fourier coefficients, as described in [START_REF] Cox | The elasticity and strength of paper and other fibrous materials[END_REF]. In this formalization, the material stiffening in the direction of main tension results from the anisotropization of the microtubule array due to the stress field. The coupling between a linear constitutive relation (Hooke's law, Eq. (2.2)) and a nonlinear dependency of the elastic properties upon stress yields an integrated formulation of the behavior, in which stiffness is expressed as a function of stress itself:

E ∞ = C ∞ w (S) -1 : S. (2.19)
This expression captures the nonlinear behavior of the wall depicted in Figs. 2.3 and 2.4.

Although growth is not included, the model can be useful in exploring the role of a stress-based feed-back on growth. Indeed, assuming that growth is a piece-wise linear function of strain (Eq. (2.1), [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Bozorg et al., 2016), the predicted strain patterns can be used as a proxy for the immediate growth.

In a situation where elasticity is isotropic or slightly anisotropic, strain and stress share the same main axis. Our FEM simulations suggest that in that case strain-based growth should favor spherical shapes, while cylinders or more complex shapes cannot be obtained (see also [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF]. Indeed, these shapes require that the main strain and stress directions are orthogonal. This is possible only if the elastic anisotropy of the material is higher than the stress anisotropy (Ineq. (2.18)). In vivo, expansion can occur in parallel or perpendicular to the maximal stress direction, depending on the type of tissue [START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF]. Our model captures this complexity thanks to the exponential form of the microtubule degradation rate (Eq. (2.10)). Indeed, this superlinear amplification of the differential of force in the microtubule distribution (Fig. 2.3(a)),

allows the system to satisfy the strain orthogonalization condition (Ineq. (2.18)) in specific cases, such as along an axis or around a hole or a bump. Fig. 2.4(c) shows that orthogonalization may happen, provided a sufficient stress amplitude f 0 , but fails with excessive stress anisotropy (2.18), Fig. 2.4(c)) due to the presence of an isotropic matrix. The regulation of the matrix through the dimensionless ratio η, should affect the mean stiffness and anisotropy of the material. Such a regulation, e.g. occurring through auxin-induced alteration of the pectin gel, has been experimentally reported (Braybrook andPeaucelle, 2013, Peaucelle et al., 2011).

Based on the strain orthogonalization condition (Ineq. (2.18)), our model allows to make predictions on the growth behavior in different biological scenarios. In the case of wound-induced stress, idealized in our simulations conducted on a punctured dome (Fig. 2.5, [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], Sampathkumar et al., 2014a), stiffening in the circumferential direction leads to maximal strain in the radial direction. Supposedly, this could enhance the efficiency of a wound healing process based on strain-based growth.

Similarly, our simulations predict that a bulging primordium generates a highly anisotropic stress pattern around itself. This, in turn, should sharply amplify growth heterogeneity between the outgrowth and its surroundings (Fig. 2.6). More generally, in the previously studied scenarios, stress feedback tends to modify the strain field, counteracting the natural tendency of pressure to induce a spherical shape, preserving the negative curvature at the periphery of an organ or promoting axial expansion (Fig. 2.5).

In addition, our model includes crucial molecular processes underlying the control of wall mechanics. The degradation term expressed by Eq. (2.10) can, in vivo, be related to katanin-based microtubule severing. Katanin is a major regulator of microtubule dynamics in both plants and animals and is involved in the formation of anisotropic microtubular arrays [START_REF] Burk | Alteration of Oriented Deposition of Cellulose Microfibrils by Mutation of a Katanin-Like Microtubule-Severing Protein[END_REF]. Although the precise mechanism is not understood, it facilitates microtubule ordering along the main stress direction [START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF]. At a molecular scale, the stabilization of CMTs by stress can be related to possible stress-induced conformational changes undergone by tubulin (Landrein andHamant, 2013a, Sampathkumar et al., 2014a), analogous to those occurring in animal cortical biopolymers [START_REF] Ghanti | Theory of strength and stability of kinetochore-microtubule attachments: collective effects of dynamic load-sharing[END_REF]. In the context of an anisotropic stress field, the tensile forces exerted on microtubules and, therefore, the efficiency of the katanin activity vary from one direction to another. We modeled this molecular mechanism as a stress-dependent expression of the microtubule depolymerization rate (Eq. (2.10)). This expression, inspired by studies on other force-sensitive biomolecular systems [START_REF] Bell | Models for the Specific Adhesion of Cells to Cells[END_REF][START_REF] Changeux | Allostery and the Monod-Wyman-Changeux Model After 50 Years[END_REF][START_REF] Nicolas | Cell mechanosensitivity controls the anisotropy of focal adhesions[END_REF], depicts the idea that the force-induced conformational change of microtubules modifies the enthalpy of the katanin-microtubule reaction. More specifically, its activation barrier is increased by the work of tensile forces exerted on microtubules, leading to their stabilization. The resulting exponential dependency of the degradation rate on tensile force enables the amplification of the microtubule response to mechanical loading. For the sake of simplicity, we neglected interactions between microtubules in the dynamics described by Eq. (2.9). In a finer model, those should be integrated as well, for they support the capacity of selfalignment of microtubule arrays [START_REF] Allard | Mechanisms of Self-Organization of Cortical Microtubules in Plants Revealed by Computational Simulations[END_REF][START_REF] Dixit | Encounters between Dynamic Cortical Microtubules Promote Ordering of the Cortical Array through Angle-Dependent Modifications of Microtubule Behavior[END_REF][START_REF] Tindemans | Survival of the Aligned: Ordering of the Plant Cortical Microtubule Array[END_REF] and could enhance the stress-sensitive mechanism we describe here.

Experimental tools to quantify microtubule distributions in plant tissues have been developed recently [START_REF] Boudaoud | FibrilTool, an ImageJ plug-in to quantify fibrillar structures in raw microscopy images[END_REF], Tsugawa et al., 2016). Combining these measurement techniques with our modeling approach, applied to high-quality digitalized images [START_REF] Cerutti | DRACO-STEM: An Automatic Tool to Generate High-Quality 3d Meshes of Shoot Apical Meristem Tissue at Cell Resolution[END_REF] would be a way to challenge some of the results exposed here, such as, the dependency of microtubule anisotropy α φ on the mean stress f 0 . It could also help estimate the transduction coefficient γ which plays a central role in the sensitivity of the feedback mechanism. This could constitute a fruitful line of future work.

Appendices 2.A The elasticity tensor as a function of the microfibril distribution

This appendix details Eq. (2.6), that establishes a relation between the fiber organization and the associated stiffness, as described in [START_REF] Cox | The elasticity and strength of paper and other fibrous materials[END_REF] .

We consider that CMFs behave like linear springs (of rest length l 0 and spring constant k). Moreover, by assuming affine deformation, their deformation reads:

∆l (θ) l 0 ≃ e θ • (E • e θ ) = E : e θ ⊗ e θ = E : Θ.
where Θ = e θ ⊗ e θ is the projector on direction θ ('⊗' depicts the tensor outer product).

The previous expression yields the total stretching energy density per unit volume of fibers (summed over all directions):

Ψ = 1 2 Y f π dθρ (θ) (E : Θ) 2 = 1 2 E : Y f π dθρ (θ) Θ ⊗ Θ : E (2.20)
with:

Y f = kl 2 0 ǫ
(where ǫ depicts wall thickness). Deriving twice the energy yields the elasticity tensor:

C f = ∂ 2 Ψ ∂ 2 E = Y f π dθρ (θ) Θ ⊗ Θ
or in the basis {e x , e y }:

(C f ) ijkl = Y f π dθρ (θ) cos (θ) ∆ ijkl sin (θ) 4-∆ ijkl with ∆ ijkl = δ i,1 + δ j,1 + δ k,1 + δ l,1
, where 'δ' stands for the Kronecker delta function.

Linearizing the cosine-sine products yields an expression of C f involving the 0 th , 1 st and 2 nd order Fourier coefficients:

[C f ] = πY f ρ 0 16     3 + ρ 2 +4ρ 1 ρ 0 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 1 -ρ 2 ρ 0 3 + ρ 2 -4ρ 1 ρ 0 2ρ 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 2ρ 1 -ρ 2 ρ 0 1 -ρ 2 ρ 0     (2.21)
where, for any natural number n:

ρn = π dθ π ρ (θ) e -2inθ and (ρ n , ρn ) = 2 (ℜ (ρ n ) , -ℑ (ρ n )) .
(2.22)

2.B Modeling cellulose deposition via CSC trajectories

We here develop a kinetic model for CSC-CMT binding, providing a chemical background to the cellulose deposition model. Here, cellulose deposition is controlled through the orientation of the CSC trajectories depicted by an a priori non-uniform π-periodic distribution µ (θ) that depicts the trajectories of proteins (regardless of their sense).

d t ρ (θ) = k 0 ρ µ (θ) -k ′ ρ ρ (θ) , (2.23)
We assume that a constant pool of CSCs is available, characterized by a constant concentration µ tot . Proteins can be either free or bound to a CMT. The free proteins are assumed to follow individual Brownian trajectories within the membrane bi-layer; for that the related cellulose deposition is isotropic, characterized by a concentration µ free .

The contribution of the bound proteins is characterized by the distribution µ bound (θ)

that depicts the trajectories of proteins on CMTs. We assume that CSCs are active regardless of their binding state, implying that cellulose deposition occurs even in the absence of CMTs [START_REF] Emons | Microtubules and cellulose microfibrils: how intimate is their relationship?[END_REF]. Microtubules solely introduce a bias in the distribution of trajectories µ. The total concentration of CSCs reads:

π µ = µ free + π µ bound = µ tot .
(2.24)

The kinetics of CSC-CMT binding is modeled by a multiplicative coupling between the concentration of CMTs in a given direction θ and the concentration of available free proteins, along with an isotropic liberation rate:

d t µ bound (θ) = k µ µ free φ ∞ (S, θ) -k ′ µ µ bound (θ) (2.25)
where k µ and k ′ µ are kinetic constants. Eqs. (2.24) and (2.25) lead to the steady-state expression:

µ ∞ bound (θ) = µ tot • K µ φ ∞ (S, θ) 1 + K µ π φ ∞ (2.26)
where K µ = k µ /k ′ µ depicts the affinity of CSCs for CMTs. By combining the previous steady-state distribution of CSCs with the expression of cellulose deposition (Eq.

(2.23)), we obtain:

d t ρ(θ) = k 0 ρ µ tot • 1 + K µ φ ∞ (S, θ) 1 + K µ π φ ∞ -k ′ ρ ρ(θ) (2.27)
that can be roughly simplified to obtain Eq. (2.7) by assuming that for all θ, K µ φ ∞ (S, θ) ≫ 1 (strong affinity):

d t ρ(θ) = k 0 ρ µ tot • φ ∞ (S, θ) π φ ∞ -k ′ ρ ρ(θ).
(2.28)

The product k 0 ρ µ tot defines the kinetic constant k ρ used in the main text.

2.C Details about the specific expression of the CMT depolymerization probability

Eq. (2.10) derives from Arrhenius law: a chemical reaction rate (k A hereafter) is a direct function of its activation energy, i.e. the energetic barrier (E A ) the system has to overcome to go from one chemical state to another:

k A = k 0 e -E A RT .
(2.29)

The product RT corresponds to the thermal energy of the considered molecules are embedded in (R and T respectively stand for the Arrhenius constant and the absolute temperature).

This general experimental law has been amended latter on in the context of biochemistry and cell adhesion in the following way [START_REF] Bell | Models for the Specific Adhesion of Cells to Cells[END_REF]. The energy barrier molecules have to over come can be lowered or raised when mechanical forces are applied to them. The resulting expression reads:

k B = k 0 e -E A ±αf RT (2.30)
where f depicts the force applied to the considered molecules and ±α stands as a coupling parameter. Clearly, the sign of this coupling parameter depends on the action of the force: if the latter eases the reaction, then the coupling parameter must be negative, so the activation barrier is lowered and the transition rate is increased.

Our working hypothesis is that mechanical forces applied to microtubules reduce their depolymerization rate; we therefore postulated a positive coupling parameters. Simplifying the notation as follows:

   k ′ φ def. = k 0 e -E A RT γ def. = α RT (2.31)
and considering that E A is proportional to force, leads to Eq. (2.10).

2.D Microtubule principal axis and anisotropy at steady regime

Generally, the CMT distribution φ is summarized through the angle of its mean orientation (θ φ ) and through its anisotropy (α φ ), respectively given by:

θ φ = - arg φ1 2 and α φ = | φ1 | φ0 (2.32)
where 'arg' is the complex argument, and { φn } n∈N depict the Fourier coefficients of φ.

We define θ ρ and α ρ in the same way, using the CMF distribution ρ.

The steady-state Fourier sequence of CMTs is derived from Eq. (2.11) and read:

∀n ∈ N : φ∞ n = c 0 π I n 2γ| f1 | e -γf 0 πK φ + I 0 2γ| f1 | e -2inθ S (2.33)
where θ S depicts the angle of main stress, and where functions {I n } n∈N depict the modified Bessel functions of the first kind [START_REF] Abramowitz | Handbook of mathematical functions[END_REF]:

∀n ∈ N, ∀x ∈ R : I n (x) = π dθ π e x cos θ cos nθ.
(2.34)

By expressing the complex argument of φ1 from the previous expression, one sees that

θ ∞ φ = θ S .
This means that CMTs align with the main direction of stress. This alignment is moreover symmetric, according to this direction (with the parametrization θ ← θ-θ S , the imaginary part of the Fourier coefficients vanishes).

The anisotropy ratio α ∞ φ (Fig. 2.3(b)) reads:

α ∞ φ = I 1 2γ| f1 | I 0 2γ| f1 | . (2.35)
which is a growing function of | f1 | (if γ > 0). According to the previous expression, the steady-state anisotropy of CMTs depends on stress only, and not on the chemical kinetics.

2.E Details on the simulation pipeline

We detail the numerical procedure employed in Section 2.5.3. All simulations have been performed in the framework presented in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], based on the open-source softwares Sofa [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF] and OpenAlea [START_REF] Pradal | Ope-nAlea: a visual programming and component-based software platform for plant modelling[END_REF]. The numerical procedure employed in this paper is adapted from [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF].

2.E.1 Mechanical equilibrium

2.E.1.1 Space discretization

Finite element procedure We discretize the continuous model through triangular finite elements. Current material positions x are expressed according to the mesh nodes q (n) based on a Lagrange P 1 barycentric interpolation:

x = n∈N N (n) (X) q (n) (2.36)
where N depicts the set of mesh nodes; N (n) (X) depicts the barycentric shape function associated with node n and evaluated at material position X. Spatial differentiation yields one 3 × 2 barycentric deformation gradient tensors per triangle τ :

F (τ ) = n∈N q (n) ⊗ ∂N (n) ∂X τ (2.37)
The 2D Green-Lagrangian strain tensor and the second Piola-Kirchhoff stress are derived from F (τ ) according to:

E (τ ) = 1 2 F (τ ) T • F (τ ) -I and S (τ ) = C (τ )
w : E (τ ) , (2.38) yielding the internal energy of the system (that is only due to deformation):

U = τ ∈T ǫA (τ ) E (τ ) : S (τ ) = τ ∈T ǫA (τ ) 2 E (τ ) : C (τ )
w : E (τ ) , (2.39) where ǫ and A (τ ) respectively stand for the wall thickness (that is uniform) and the surface area (reference configuration) of triangle τ . Nodal inner forces are obtained from the gradient of the strain energy with respect to node positions:

f (n) int = - ∂U ∂q (n) .
(2.40)

The external loads are derived from pressure and expressed at each node n according to its adjacent triangles T n :

f (n) ext = P 3 τ ∈Tn a (τ ) n (τ ) , (2.41)
where n (τ ) stands for the outer normal of triangle τ ; a (τ ) is the surface area of triangle τ in the current configuration.

Meshes All meshes have been generated through the open source software Blender (blender.org). At the exception of the cylindrical stem (for which the original mesh was satisfactory), they have been enhanced via a surface smoothing procedure implemented in the non-commercial meshing tool Graphite (alice.loria.fr/software/graphite).

2.E.1.2 Time discretization

Computing the equilibrium (Eq. (2.3)) consists in minimizing the total energy U (see [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], for further details). This is performed by solving the equation of motion until node forces balance. We use the backward Euler method and the conjugate gradient to perform this step.

To improve stability and to avoid excessive stiffness inhomogeneity, the tensor S employed in Eq. 2.11 is modified on each element and integrates a weighted average of the neighboring stress tensors.

2.E.2 Dynamical resolution of Eq. 2.7

The dynamics of the CMF distribution is computed by an iterative update of the first three Fourier coefficients of ρ through the following forward Euler time scheme with constant time step ∆t (Eq. 2.7 and Fig. 2.8):

∀n : ρt+∆t n = 1 -k ′ ρ ∆t ρt n + 2k ρ ∆t π φ∞ n S t φ ∞ 0 S t .
(2.42)

The stiffness tensor C w is subsequently updated according to Eqs. 2.4 to 2.6 and 2.33.

At t = t 0 , the system is set to its null-stress equilibrium, derived from from Eqs. 2.11 and 2.14:

φ (θ, t 0 ) = c 0 K φ 1 + K φ π and ρ (θ, t 0 ) = K ρ π .
(2.43)

Figure 2.8: Simulation pipeline. The initial state corresponds to the null stress state (no residual stress and no loading). The elastic equlibrium is computed after loading the structure with pressure P > 0. At each elastic equilibrium, the mechanical properties are updated, which brings the systems out of equilibrium.

2.E.3 Visualization of the CMF distributions in Figs. 2.5 to 2.7

We adapted the concept of nematic tensor from the physics of liquid crystals (De Gennes and Prost, 1995) by defining:

Q φ = 4 π dθ π φ (θ) Θ = φ 0 1 + φ 1 φ 0 φ1 φ 0 φ1 φ 0 1 -φ 1 φ 0 .
(2.44)

Notice that, as stated in 2.5.1, the previous definition builds a formal equivalence be-

Q φ,1 -Q φ,2 Q φ,1 + Q φ,2 = 2| φ1 | φ 0 = | φ1 | φ0 (2.45)
Q φ is equivalent to an ellipse with main axis oriented by θ φ = -1 /2 arg φ1 (Eq. 2.32), and with aspect ratio and semi-major axis respectively given by (1α φ ) / (1 + α φ ) (with

α φ = | φ1 |/ φ0
) and φ 0 (1 + α φ ) (see Fig. 2.9).

Figure 2.9: Ellipse-based visualization of the nematic tensor associated with the distribution φ.

Chapter 3

Coupling growth and stress feedback

Introduction

In the previous chapter, we developed a model of the stress feedback valid in the elastic regime (namely neglecting growth). In this chapter, we detail a method to couple the previous model to growth, using the growth model developed in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF].

By doing so, we model the complete feedback loop that was introduced in Chapter 1 and illustrated in Fig. 1.11). For this coupling to be physically consistent, we model the advection-like deformation undergone by the fiber arrays as the infinitesimal volumes of the tissue deform.

Here, we will study the dynamical properties of the system (stress feedback + growth), in particular, in the context of symmetry 1 breaking and symmetry maintenance, by using ideally-shaped pressurized shells. Indeed, the fact that the feedback exploits stress anisotropy to generate anisotropic growth suggests that it has the potential to break an initial symmetry of shape by responding to slight geometrical asymmetries. To study this effect, we here consider the highest level of symmetry, that is the spherical shape, and we assess its stability under various kinds of feedback parameters and shape perturbations. We explore the properties of the model focusing on two main key parameters, i.e. the sensitivity and rapidity of the feedback (detailed hereafter), that need to be modulated in order to promote efficient asymmetry amplification. From a biological standpoint, this kind of process may be involved in the early stages of organ development, for example allowing for small roundish bumps to grow towards axisymmetric, highly elongated organs like branches.

Model

Modeling growth and growth-induced fiber advection

Growth is modeled as in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF]. We consider a morphoelastic material that undergoes rapid elastic relaxation and slow plastic expansion (see Section 1.2.2), both combined through the multiplicative decomposition of the deformation described by Eq. (1.6). Inasmuch as the elastic and plastic regimes occur with very different time scales, one can simulate the incremental evolution of the system's geometry through a two-timing procedure, consisting of i/ computing the elastic equilibrium of the system and ii/ modifying the reference of the material.

The elastic regime of the material is already described in Chapter 2, in which we have detailed a constitutive relation for the material (see Eqs. (2.4) to (2.7) and (2.11)) taking into account stress-induced CMT reorganization and CMF turnover. As for the plastic regime, we use the strain-based growth expressed by Eq. (1.12) [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], and discretized by the following forward Euler procedure:

F g (t + ∆t) = I + ∆tL g (t) • F g (t) (3.1)
where L g (t) depends on the Green-Lagrangian strain E (t)

def.

= 1 2 F T e (t) • F e (t) -I at step t as expressed by Eq. (1.12); ∆t is a constant time step size. The factor f g (t)

def.

= I + ∆tL g (t) depicts the incremental growth that is applied at each growth time step.

In principle, growth induces a change in the fiber density and orientation, which can be assimilated to an advective process. In fact, by virtue of mass conservation, the volumetric density of CMFs decreases as the wall expands if we assume that no cellulose intake is involved. Moreover, anisotropic growth modifies the orientation each fiber, which rotates towards the direction of main deformation.

To model these effects, we assume that the deformation of the fibers is affine, i.e. that these deform according to the macroscopic medium (they neither rotate nor slide relative to the material). The growth procedure detailed previously (Eq. (3.1)) naturally provides an incremental way to simulate the advection of the CMFs. At each growth step, the medium is updated with a small incremental deformation f n g that affects the fiber density function (N.B.: to simplify notations, we drop the exponent n in the rest of this paragraph, since all the next developments are made at time step n). Let R (Θ) and ρ (θ) be the angular density of CMFs (in the reference configuration) respectively before and after applying the incremental deformation f g ). Given a fiber oriented in direction θ after growth, Θ will depict the angle of this fiber before growth (see Fig. 3.1).

Under the affine deformation assumption, the conservation of the cellulose mass that is initially contained within a control volume dV (dv after deformation) and infinitesimal angular sector [Θ, Θ + dΘ], is given by:

ρ (θ) dθdv = R (Θ) dΘdV (3.2)
which amounts to:

ρ (θ) = J -1 R (Θ) dΘ dθ (3.3)
where J def.

= dv dV = det f g ≥ 1 measures the variation of volume. By detailing dΘ dθ (see Appendix 3.A), one obtains:

ρ (θ) = λ 2 1 + λ 2 2 2 - λ 2 2 -λ 2 2 2 cos 2 θ -1 R (Θ) (3.4)
where λ 1 and λ 2 are the eigenvalues of f g (with λ 1 ≥ λ 2 ); and angle θ is the phase between θ and the principal angle of f g . The behavior predicted by Eq. (3.4) is further analyzed in Section 3.3.1.

For the sake of simplicity, we do not consider any growth-induced loading of the fibers, that would imply that the strain of each fiber would not be solely a function of its orientation, but also of its age and mechanical history. For example, two fibers deposited (unloaded) at two different times and in the same orientation could in principle be loaded in a different manner, depending on how much the wall has grown between the two fiber insertions. Here the strain energy associated with each fiber reflects the macroscopic elastic strain of the material, only. This is equivalent to considering that fiber strain energy relaxes at a longer time scale, in the same fashion as the macroscopic medium. Using this assumption, one can use Eq. (2.6) to relate the fiber density to stiffness.

Simulation

Triangular membrane finite elements equipped with P 1 -Lagrangian shape functions are used to model the walls. To compute equilibrium, the equation of motion on nodes is dynamically integrated until the pressure forces f ext and reaction forces f int balance.

Equilibrium is found when the maximal node force falls below a threshold. This step Sofa [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF].

The reaction forces are evaluated from the expression of the total elastic energy of the system:

U = τ ∈T ǫA (τ ) Ψ (τ ) (3.5)
where T is the set of triangle elements; ǫ depicts the wall thickness (assumed to be uniform); A (τ ) is the surface area of triangle τ in the stress-free configuration. The strain energy density Ψ (τ ) at triangle τ is defined as:

Ψ (τ ) def. = Ψ E (τ ) = 1 2 E (τ ) : C (τ )
w : E (τ ) .

(3.6)

By derivation, one obtains the reaction force f

(n)

int at node n associated with node position q (n) :

f (n) int = - ∂U ∂q (n) = - τ ∈T ǫA (τ ) ∂Ψ (τ )
∂q (n) .

(3.7)

The nodal external loads (pressure forces) are computed according to:

f (n) ext = P 3 τ ∈Tn a (τ ) n (τ ) (3.8)
where n (τ ) ( n (τ ) = 1) and a (τ ) are respectively the outward normal and surface area of triangle τ (in the current configuration); T n is the set of the triangle elements that contain node n; P is the excess of pressure relative to the outside.

After equilibrium is reached, we compute the growth step. Two computational solutions can be used for this operation. The first one consists of considering each finite element in its reference configuration, isolated from its neighbors. After growing each element separately, a compatible geometry can be computed by reconnecting the elements. This solution was used in [START_REF] Bassel | Mechanical constraints imposed by 3d cellular geometry and arrangement modulate growth patterns in the Arabidopsis embryo[END_REF] for example. Alternatively, we can modify the elastic deformation gradient

F (τ )
e in order to compensate for a change in the reference shape of the elements, which, from a computational standpoint, is a simpler operation. This solution is employed here.

The state of the system at mechanical balance and at time step t -noted M (t) -is defined through the list of node positions, and element deformation gradients:

M (t) def. = q (n) (t) n∈N , F (τ )
e (t)

τ ∈T

(3.9)

where N depicts the set of nodes. Our general problem is to compute the evolution of the system, i.e. computing M (t + ∆t) from M (t). For this, we re-express Eq. (3.1) according to F e (t), using Eq. (1.6), as detailed hereafter. Algorithm 1 and Fig. 3.2(a)

summarize the procedure employed to update the system.

Starting from a given equilibrium M (t), we compute the rate of plastic growth L (τ )

g (t) from F (τ )
e (t) for all triangles τ (using Eq. (1.12)), which allows to compute F

F (τ ) (t) = F (τ ) e (t) • F (τ ) g (t) = F (τ ) e (t + ∆t) • F (τ ) g (t + ∆t) (3.10)
where F (τ ) e (t + ∆t) is given by: Tensor F (τ ) e (t + ∆t) defines the new elastic deformation of the system after application of growth. The tilde notation used in Eq. (3.11) emphasizes that, a priori, the state

F (τ ) e (t +
M (t + ∆t) def. = q (n) (t) n∈N , F (τ ) e (t + ∆t) τ ∈T
is not at mechanical equilibrium.

Instead, it is used as an initial condition to compute the next equilibrium M (t + ∆t) (Algorithm 1).

After L (τ ) g (t) is computed for each triangle, the growth-induced fiber transport is applied provided the incremental deformation f (τ ) g (t) at step t, and by using Eq. (3.4). Since Eq. (3.4) is expressed in the angular domain (rather than the frequency domain) we first compute the Inverse Fast Fourier Transform of the CMF distribution R (τ ) , then the Fast Fourier Transform to update the Fourier coefficients. In practice, the angular distributions ρ (τ ) and φ (τ ) are approximated by their first 30 Fourier coefficients (Eqs. (2.22) and (2.33)). Although this operation is the costliest operation of the growth step, the computation of the mechanical equilibrium is by far the most important bottleneck of the full algorithm.

The last step consists of computing CMF deposition (according to Eq. (2.42)) and update the stiffness tensor (Eqs. (2.4) to (2.6)). This is done by the classical Runge-Kutta (RK4) procedure (assuming constant stress during the integration step).

For better spatial regularity of in the CMF density, the stress computed at each element is smoothed over its neighbors in a simple manner. The smoothed stress is used to evaluate Eq. (2.11). Since the PK2 stress is expressed in the 2D material coordinate system, it cannot be directly used to perform the smoothing. This is done for each triangle τ , by computing the stress σ (τ ) in the current configuration (Cauchy stress) from the PK2 stress S (τ ) :

σ (τ ) def. = 1 det F (τ ) e F (τ ) e • S (τ ) • F (τ ) e T (3.12)
Then, a mean stress is computed over the element neighborhood (weighted by triangle area). This stress, is generally of rank 3, as the neighborhood may be slightly non planar. It is then re-projected onto the mesh triangles and pulled back to the reference configuration, providing a smoothed PK2 stress field.

Algorithm 1 Growth increment

Require: M = (x, F e ) at equilibrium Ensure:

M ← Grow (M ) is the new equilibrium configuration after growth

1: E ← 1 /2 F T e • F e -I ⊲ Elastic strain from the current equilibrium configuration 2: L g ← Φ E -E y + ⊲ Strain-based growth 3: F e ← F e • I + ∆tL g -1
⊲ Redefining elastic deformation. Leaving equilibrium 4: M ← Equilibrium (x, F e ) ⊲ (FEM)

Meshes

Results

Verification of the model behavior

Convergence of the growth method In this section, we investigate the convergence of the growth method with respect to mesh resolution. In fact, if the solution excessively depends on the discretization, the reliability of the algorithm in accurately resolving the continuous model is not guaranteed. For this, we test the algorithm using spherical meshes with increasing resolution, from 40 to 5,000 nodes (Fig. 3.3(a)), in order to test to what extent an increase in mesh quality affects the various simulated properties. As we augment the mesh resolution, convergence is achieved as soon as the solution is no longer significantly affected by an improvement of the mesh quality.

For each version of the sphere, we simulated growth until the mesh surface area reached twice its initial value. For a given quantity p, we define (for each mesh m) the following measure of the deviation of p m relative to a reference, at time t:

ε m p (t) = |p m (t) -p ref (t)| p ref (t) (3.13)
where p ref (t) is the value associated with the reference mesh (the finest mesh of our set) at time t. We monitored the mesh radius 2 (ε R ), mean strain (ε E ), average CMF concentration (ε CMF ) and isotropy 3 (ε isotropy ).

The results are given in Fig. 3.3(b). For extremely coarse meshes (e.g. the two first meshes on Fig. 3.3(a)), the error (for all the considered quantities) is confined to less than 10%. For finer meshes, the error falls below 1%, which we can consider as a satisfactory error for the present study. Figure 3.3(b), shows that the error tends to zero as a decreasing function of the mesh quality. This indicates that the growth algorithm (considered without stress feedback) in the case of a spherical shell, is convergent with respect to the mesh resolution. In the sequel, we employ meshes with an apparent resolution roughly comparable to that of the third mesh on Fig. 3.3(a).

2 Computed as R def.

= 3 Volume Area . 3 With isotropy def.

= 1 -|ρ1|/ρ0 triangles , see Eq. (2.22).

Fiber advection

We see from Eq. (3.4), that the CMFs tend to accumulate in the direction of main deformation, and are depleted in the perpendicular direction. This captures a tendency of CMFs to passively align with the direction of main growth. In addition, the model captures a dilution of the CMFs: these tend to locally vanish due to the increase in wall volume (see also Fig. 3.5). This can be evidenced from Eq. (3.4) in the particular case where R is an isotropic distribution (say equal to R 0 ). In this case, one may verify that the dilution factor is the variation of volume J:

π ρ = R 0 π λ 2 1 + λ 2 2 2 - λ 2 1 -λ 2 2 2 cos 2 θ -1 dθ = πR 0 J (3.14) (see Appendix 3.B).

Growth of stress-adaptive shells

General idea In Chapter 2, we suggested that the stress feedback may be a mechanism to amplify shape asymmetry. For this reason, it may have the potential to promote symmetry breaking, by amplifying small perturbations in the initial shape. To test this behavior, we simulated the growth of simple abstract surfaces, with various degrees of stress feedback. We considered key parameters of the stress feedback, that are the sensitivity coefficient γ (that governs the sharpness of CMT alignment, see Eq. (2.35)) and the kinetic coefficients k ρ and k ′ ρ that govern cellulose turnover (for the sake of simplicity we impose K ρ = k ρ /k ′ ρ = 1 mol.µm -3 , which allows to vary only one parameter). Parameter k ′ ρ governs the rapidity with which the tissue modifies its stiffness in response to a stress, and characterizes the typical time of CMF reorganization. Intuitively, a sensitive CMT system coupled with a very slow CMF dynamics will fail in rapidly responding to stress by modifying stiffness. Conversely, we may expect that a rapid cellulose renewal will be of little interest if CMTs insufficiently coalign along stress main axis.

As in Chapter 2, we considered pressurized membrane meshes (no inner tissues). Initially, the cellulose system is set in the trivial null-stress steady concentration, which mostly depends on K ρ (see Eq. (2.43)). Hence, all simulations are initialized with the same isotropic elasticity. Simulations are conducted until the mesh surface area reaches twice its initial value.

Spherical mesh

The most symmetric shape is the sphere. Evidently, by a trivial symmetry argument, the sphere is a fixed point of the system. We here study its stability.

We simulated stress feedback and growth using a spherical mesh composed of about 4,000 nodes and 8,000 triangular elements. The results of the simulations on the sphere are shown in Fig. 3.6. Confirming the intuition given above, we verify that for a slow and/or insensitive feedback, growth is uniform and spherical, which means that the sphere is stable for this range of parameters. By contrast, by using a fast and sensitive feedback (upper-right corner of Fig. 3.6), we observed strong growth heterogeneities due to rapid anisotropization of the tissue, suggesting that the spherical shape will be unstable for this range of parameters. Yet, in absence of a remeshing procedure, this kind of heterogeneities occurred with a rapid deterioration of the mesh quality, as quantified by the histograms in Fig. 3.6, showing the distributions of triangle skewness and mesh roughness (these quantities are defined in Appendix 3.C). For extremely fast and sensitive feedback, maximum skewness and roughness in fact exceeded 50%.

To assess the nature of this apparent "patterning" in the growth rate, we tested the sensitivity of this result with respect to mesh quality, using various spherical meshes with a node number ranging between 3,000 and 20,000 (see Fig. 3.7). At the end of the simulation, we quantified for each mesh the discrete correlation length (defined in Appendix 3.D) which provides an estimate of the typical length of the heterogeneity.

For all simulations, this length was comparable to the typical size of the triangular elements, although the number of elements varied significantly between all meshes (see Fig. 3.7). This means that, for this range of extreme (potentially non-physiological) parameters, our FEM pipeline fails in converging with respect to mesh size. Note also that the smoothing operation presented in Section 3.2.2 is performed over the immediate neighbors of a given triangle. Therefore the smoothing radius also varied with mesh size. In future work, we should investigate a more sophisticated, potentially mesh-independent stress-smoothing procedure. Procedures to improve the quality of the computed stress should be also investigated.

Quasi-spherical spheroids

The previous example suggests that the spherical shape is globally stable, although there exists a range of parameters of the feedback that will probably lead to instabilities (confirming this result would require further numerical investigations). The only perturbations considered in that case were those introduced by the meshing, which approximates a sphere by a polyhedron. To study the effect of larger perturbations, we considered the case of almost-spherical spheroids: a prolate spheroid (elongated sphere) and oblate spheroid (flattened sphere). In both cases, we considered a 5% perturbation in diameter.

The stress field in both cases was computed (see Fig. 3.8), showing that in the prolate case, hoop stress slightly dominates, while in the oblate case, the dominating stress is mostly axial. This can be explained by the different natures of the curvature in these two shapes, although both are very close to a sphere. A stress feedback has then the potential to yield utterly different behaviors depending on whether the perturbation is prolate or oblate.

In the case of the slightly prolate spheroid (Fig. 3.9), the system remained roughly spherical provided a weak feedback, as the material stiffness remained more or less isotropic. In the case of an extremely sharp and rapid feedback, growth became highly heterogeneous and noisy as previously observed. A slightly weaker feedback (e.g. γ = 10 MPa -1 , k ρ = 1 mol.µm -3 .s -1 ) allowed for a response to the slight initial stress anisotropy (as CMTs aligned circumferentially) and for the amplification of the initial shape asymmetry. In fact, the shape elongated as a cylinder. Here the mesh skewness increased while, comparatively, mesh roughness was small, which is the mark of a cohesive anisotropic growth.

In the case of the slightly oblate spheroid (Fig. 3.10) the system again developed towards a sphere in the case of a slow or insensitive feedback. Spheroid flatness was slightly amplified for more sensitive and faster feedback, yet this often occurred together with a slight elongation of the shape in the axial plane, that is the plane orthogonal to the initial symmetry axis (not shown). This suggests that the feedback tends to destabilize flat shapes by breaking the initial axial symmetry, and promotes elongation. In these simulations, the stress feedback tended to elongate triangles in a very heterogeneous manner (as shown by the histograms in Fig. 3.10). In particular, singular growth spots appeared at the poles of the structure.

Growth of highly asymmetric shapes

In the previous paragraph, we saw that a slightly prolate sphere may develop axially. We next confirmed this result on a more elongated shape (a cylindrical "pill", shown in Fig. 3.2(c)). In this case, axial cylindrical expansion was promoted by a sharp feedback (γ > 15 MPa -1 ), see Fig. 3.11(a), as we observed an exclusively-axial growth, resulting in a cylindrical expansion with quasiconstant radius. By contrast, lower feedback sensitivity resulted in a radial expansion.

A similar set of simulations, starting from a shape displaying a concave zone (here the shape shown in Fig. 3.2(d)) also confirmed the suggestion made in Chapter 2, that a stress feedback is a potential way to maintain a concave frontier (see Fig. 3.11(b)). In fact, for sufficient values of γ, we observed that the frontier radius barely increased, while the rest of the tissue dilated substantially. This may be a way to amplify the growth heterogeneities between two developmental zones.

Simulation of highly oblate shapes (flat spheroids) generally did not properly converge, even for extremely fine meshes and extremely small time step size. This may partly be due to the emergence of compressive stresses and buckling at the ridge of the spheroid.

Membrane finite elements (that do not take bending energy into consideration) are not appropriate to represent this kind of effect. We further discuss this special case -in particular its biological relevance -in the discussion of this chapter.

Discussion

Amplification of tissue asymmetry by a stress feedback In this chapter, we have studied to what extent a stress feedback can amplify initial perturbations in shape. This is possibly required in the establishment of highly asymmetric organs such as stems or leaves. In vivo, these organs start as small, more or less roundish bumps, which result from a local modification of the epidermal cell wall constitution (see [START_REF] Sassi | An Auxin-Mediated Shift toward Growth Isotropy Promotes Organ Formation at the Shoot Meristem in Arabidopsis[END_REF].

These bumps may further develop in a highly anisotropic manner.

To study how this transition in growth mode occurs, we used simple abstract shapes, i.e. perturbed spheres, as initial shapes. In particular, we employed the classic pressurized vessel assumption, that consists of neglecting the mechanical role of inner tissues and of considering only the outermost wall as an elastic membrane.

This theoretical study suggested that the stress feedback has the potential to promote axial elongation starting from a slightly prolate spheroid (Fig. 3.9), provided a reasonably rapid and sensitive feedback is active, and, secondly, to maintain axial growth while keeping a nearly constant radius (Fig. 3.11(a)). This result is consistent with what was found in [START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF], [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF], that are based on more ad hoc stress feedback models. This amplification was less clear in the case of an oblate spheroid and flat expansion, the simulation of which was more difficult to conduct, due for example to fast and localized mesh deterioration (Fig. 3.10). Further numerical investigations -including the development of a proper remeshing algorithm -will be needed to validate this result. To date, the implementation of a triangle subdivision method is ongoing.

One can still expect, based on mere intuition, that a hollow membrane model will fail in producing very thin shapes. In principle, broadening should result in an increase of the stress in the flattest parts of the structure (namely at the poles), the curvature of which should tend to zero. In the most extreme case, the maintenance of an almost null curvature will require polar tissue tension to approach infinity, in order to balance turgor forces. This is of course very unlikely to happen in plants. In reality, part of this stress should be transferred to inner tissues, that are neglected in the pressure-vessel hypothesis. In such a scenario, one can expect that inner walls may have a crucial mechanical role in maintaining and amplifying flatness. Whether a stress feedback can account for this kind of shape is treated in the next chapter, in the case of Arabidopsis thaliana and tomato leaves and sepals that (as in many other species) are extremely asymmetric organs.

On the maintenance of a homeostatic cellulose density

To achieve efficient shape amplification during the transient organ initiation, fast and sharp feedback may be necessary, as abstracted through the idealized case of a perturbed sphere. In our model, these features are represented through parameters γ and k ′ ρ . We hereafter discuss the role and interpretation of k ′ ρ from mathematical and biological standpoints.

The equation governing the dynamics of cellulose density (Eq. (2.7)) is composed of two terms in its right-hand side: a term representing the CMT-guided polymerization of cellulose, which models a rather well identified biological process (i.e. the CMTguided deposition of cellulose by CSCs, see [START_REF] Paredez | Visualization of Cellulose Synthase Demonstrates Functional Association with Microtubules[END_REF], and an isotropic linear decay. The latter was necessary in Chapter 2 to guarantee the existence of a steady cellulose concentration. However, its mathematical necessity is here less obvious, in so far as the growth-induced advection also contributes to reducing CMF concentration c ρ = π ρ. This can be expressed in the form of a classic advection-reaction equation (that directly derives from the continuity equation):

d t c ρ = k ρ - k ′ ρ c ρ Cellulose decay -c ρ ∇ • v Advection (3.15)
(where d t depicts the material derivative, see [START_REF] Holzapfel | Nonlinear solid mechanics[END_REF]. This shows that growth has the potential to induce a decrease of cellulose concentration, yet in a passive manner (since v is an emergent property of our system, in contrast to k ′ ρ ).

While previously, the steady regime depended upon the independent values of k ρ and k ′ ρ , advection depends upon the growth rate, and therefore is regulated through stiffness and CMF density itself, probably in a less direct manner than a kinetic coefficient like k ′ ρ . One immediate question is then: can advection only (i.e. taking k ′ ρ = 0) overcome cellulose deposition, so that the cellulose concentration can be maintained in a homeostatic, steady value?

Intuitively, the answer should be negative. Indeed, we may expect that such a system (advection + polymerization with constant rate) would be in principle hard to stabilize.

A large value of k ρ would most plausibly result in growth arrest, since stiffness would increase and, likewise, advection should vanish. Conversely, a low value of k ρ would imply soft walls and fast growth, and therefore fast depletion of the fibers, potentially at the expense of tissue integrity. In addition, it is not certain that advection would allow a rapid turnover of the fibers allowing for a rapid response of the system to a change in stress.

To rule out this issue, we can picture two different parsimonious mechanisms. First, we may keep the depolymerization term k ′ ρ ρ (θ) used in Eq. (2.7), which allows an arbitrarily rapid saturation of cellulose concentration and provides flexibility to the system. This was the solution chosen in this chapter and the previous one. This term can reflect a depolymerization process, possibly controlled by the cell through some enzymatic activity. To date however, no such depolymerization mechanism has yet been identified to our knowledge. In the alternative scenario, cellulose polymerization kinetics could depend itself upon growth rate and/or current cellulose concentration.

In other terms, the cell may dynamically regulate the quantity of cellulose it synthesizes by probing its growth and/or structural state and controlling k ρ accordingly. This might also be a plausible scenario, supported by the recent identification of several potential wall integrity and growth sensors like THESEUS or FERONIA involved in the control of cellulose secretion [START_REF] Cheung | THESEUS 1, FERONIA and relatives: a family of cell wall-sensing receptor kinases?[END_REF][START_REF] Hématy | Novel receptor kinases involved in growth regulation[END_REF], Verger and Hamant, 2018[START_REF] Wolf | A receptor-like protein mediates the response to pectin modification by activating brassinosteroid signaling[END_REF]. However, in this scenario, it is not clear how a cell could rapidly renew its wall composition and efficiently perform a change in its growth behavior. Possibly, both scenarios may be actually valid.

Note that we assumed that growth rate depends on fiber concentration mostly. This is a coarse approximation, as, for example, growth is also mediated through the hemicellulose tethers that connect the fibers together (Cosgrove, 2005, Geitmann and[START_REF] Geitmann | Mechanics and modeling of plant cell growth[END_REF] as modeled in [START_REF] Dyson | A model of crosslink kinetics in the expanding plant cell wall: Yield stress and enzyme action[END_REF]. These were here abstracted through the isotropic elasticity tensor C g . To fully understand the loop between wall remodeling and growth, and its integration in global morphogenesis, finer models of the walls and its growth, compatible with a global simulation of mechanics, are required.

Appendices

3.A Proof of Eq. (3.4)

I here detail the establishment of Eq. (3.4) that describe the deformation of the CMF distribution due to the deformation increment f g = I + ∆tL n g .

In particular we determine an expression for dΘ dθ involved in the mass conservation equation (Eq. (2.24)). By definition:

Θ = arctan 2 (j • N , i • N ) = arctan 2 (j * • n, i * • n) (3.16) where i * = f -T g • i and j * = f -T g • j. Recalling that ∇ arctan 2 (y, x) = -y
x 2 +y 2 , x x 2 +y 2 , one obtains, by differentiating Eq. (3.16):

dΘ dθ = dn dθ • dΘ dn = n π 2 • -(j * • n (θ)) i * + (i * • n (θ)) j * i * • n (θ) 2 + j * • n (θ) 2 with n π 2 = n (θ + π /2) = n π 2 • (i * ⊗ j * ) T -i * ⊗ j * f -1 g • n 2 • n = i * ⊗ j * -(i * ⊗ j * ) T f -1 g • n 2 • n π 2 • n = i * ⊗ j * -(i * ⊗ j * ) T f -1 g • n 2 : n ⊗ n π 2 .
From mass conservation (Eq. (2.24)), one derives:

ρ (θ) = R (Θ) λ (θ) 2 J i * ⊗ j * -(i * ⊗ j * ) T : n ⊗ n π 2 . (3.17)
where Θ is obtained from Eq. (3.16) and where λ (θ) = 1 / f -1 g n is the stretch in direction θ.

We can assume, without loss of generality that f g is diagonal in B (f g is a symmetric tensor):

f g B = λ 1 0 0 λ 2 (3.18)
Hence:

i * ⊗ j * = f -1 g • i ⊗ f -1 g • j = 1 J i ⊗ j (3.19)
giving:

i * ⊗ j * -(i * ⊗ j * ) T : n ⊗ n π 2 = 1 J . (3.20)
In addition:

λ (θ) 2 = 1 f -2 g : n ⊗ n = 1 cos θ λ 1 2 + sin θ λ 2 2 = 2λ 2 1 λ 2 2 λ 2 1 +λ 2 2 1 - λ 2 1 -λ 2 2 λ 2 1 +λ 2 2 cos (2θ)
.

(3.21)

Combining the two previous expressions with Eq. (3.17) finally gives:

ρ (θ) = λ (θ) J 2 R (Θ) = 2 λ 2 1 +λ 2 2 1 - λ 2 1 -λ 2 2 λ 2 1 +λ 2 2 cos (2θ) R (Θ) . (3.22)
For the previous expression to be valid in any basis, the argument of the cosine must take into account the phase shift that corresponds to the angle θ f g of main deformation. This is done by replacing θ by θ = θθ g in the previous equation, where θ g is the angle of main growth.

3.B Proof of Eq. (3.14)

Let's note A = λ 2 1 +λ 2 2 2 and α = λ 2 1 -λ 2 2 λ 2 1 +λ 2 2
within the scope of this appendix. We compute:

π dθ A 1 -α cos 2 θ = 1 2A 2π dθ 1 -α cos 2θ = 1 2A 2π dθ 1 + α -2α cos 2 θ = 2 A π 2 0 sec 2 θdθ (1 + α) sec 2 θ -2α sec θ = 1 cos θ = 2 A (1 + α) π 2 0 sec 2 θdθ 1-α 1+α + tan 2 θ sec 2 θ = 1 + tan 2 θ = 2 A (1 + α) ∞ 0 dx 1-α 1+α + x 2 (subst.: x = tan θ) = π A √ 1 -α 2 = π J .

3.C Quantifiers of mesh quality

For each triangle i of a mesh, we define the triangle skewness S i as:

S i = 1 - l min i l max i (3.23)
where l min i and l max i respectively depict the smallest and largest edge length of the triangle. Optimality (S i = 0) is achieved only for equilateral triangles. The skewness is in fact a local quantifier of local triangle shape.

For each mesh node n we define the roughness R n as:

R n = 1 - A min A max (3.24)
where A max and A min respectively depict the surface areas of the biggest and smallest triangles containing node n.

3.D Measure of discrete correlation length

We consider a scalar quantity s defined at mesh nodes. We are interested in quantifying the spatial variation of s. For that we compute the discrete correlation function of s.

For each pair of node (i, j) ∈ N 2 (N is the set of all nodes), we define the topological distance d (i, j) as the minimal number of triangle edges forming a path between nodes i and j (Fig. 3.12(a)). This distance allows to define the topological circle of center i and radius L ∈ N as:

C (i, L) = {j ∈ N |d (i, j) = L} .
(3.25) (Fig. 3.12(a)).

We define the discrete autocorrelation function of s as:

a s (L) = 1 #N i∈N j∈C(i,L) (s i -s) (s j -s) #C (i, L) with s = 1 #N i∈N s i . (3.26)
This extends the classical definition of the autocorrelation function to a mesh. We define the discrete correlation length of s by performing an exponential fit on the discrete values a s (L), which provides a characteristic decay distance (Fig. 3.12(b)).

Abstract

Many plant species have thin leaf blades, which is an important adaptation that optimizes the exchanges with the environment. Here, we provide evidence that their three-dimensional geometry is governed by microtubule alignment along mechanical stress patterns in internal walls. Depending on the primary shape of the primordium, this process has the potential to amplify an initial degree of flatness, or promote the formation of nearly axisymmetric, mostly elongating organs, such as stems and roots.

This mechanism may explain leaf evolution from branches, which is alternative to Zimmermann's influential, but widely questioned, telome theory.

Main text

The formation of thin leaf lamina in plants is an important adaptation that optimizes vital processes, including photosynthesis, transpiration and respiration (Maugarny-Calès and Laufs, 2018). While the regulatory genetic network controlling leaf polarity has been well characterized [START_REF] Kuhlemeier | The Sussex signal: insights into leaf dorsiventrality[END_REF], comparatively little is known on how such a thin structure mechanically arises and maintains itself during development. We addressed this issue by combining computational modeling and a three-dimensional (3D) experimental analysis of leaf morphogenesis in two species (Arabidopsis thaliana and tomato, Solanum lycopersicum). Various leaf types (rosette and cauline leaves, cotyledons and sepals) were analyzed.

Primordia of leaves and leaf-like organs initiate from apical meristems, as rounded, slightly asymmetric bulges (Fig. 4.9c and d). Starting from a ratio of blade width (in the mediolateral axis) to thickness (in the dorsoventral axis) between 1.5 and 2, the leaf and sepal primordia mainly expand in two dimensions, forming a thin lamina with ratios of 10-12 in sepals and even higher in leaves (Poethig and Sussex, 1985) (Fig. 4.1a, d and Fig. 4.9). Growth directions largely rely on the orientation of the cellulose microfibrils in the cell walls [START_REF] Ali | Physical Models of Plant Development[END_REF], Sampathkumar et al., 2014a), which depends on the organization of the cortical microtubule (CMT) arrays guiding the cellulose synthase complexes [START_REF] Paredez | Visualization of Cellulose Synthase Demonstrates Functional Association with Microtubules[END_REF].

To investigate the role of CMTs in leaf development, we first characterized CMT arrangements using immunostaining and in vivo confocal imaging (Fig. 4.1b and e to j). CMT behavior along the inner and outer periclinal walls was highly dynamic. In very young growing sepals at stage 3/4, these CMTs transiently showed some degree of anisotropy (Fig. 4.1e, f), which decreased significantly (from 30% to 10%) early in development (Fig. 4.1h-j; see also [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF]. Similarly, the very young leaf also transiently showed aligned CMTs along the outer periclinal walls, which became more disorganized afterwards (Fig. 4.10a to d).

A very different behavior was found along most of the anticlinal walls. Here CMTs were mainly oriented perpendicular to the surface in Arabidopsis cotyledons, leaves and sepals, as well as in tomato leaves Fig. 4.10e and Fig. 4.11a,b). Staining of cellulose confirmed that this coincided with the main microfibril orientation in stress were dominantly mediolateral (i.e. along the second axis of the ellipsoid). Accordingly, the virtual structure evolved towards a spherical shape if no mechanical feedback was introduced (simulation 1: Fig. 4.2b, d), which echoes the results obtained in vivo using oryzalin treatment. Conversely, when the stress feedback was active throughout the entire 3D tissue, the structure grew longer (and slightly wider), showing that a stress feedback has the potential to promote anisotropic expansion (simulation 2: Fig. 4.2b, d, outer and inner feedback active). In line with the experimental evidence, these simulations showed CMT alignments along the dorsoventral axis on anticlinal walls (Fig. 4.2c, left panel), restricting growth in this direction.

However, simulations where the feedback was active throughout the entire tissue, systematically predicted a mediolateral alignment of CMT arrays along the stiffer outer walls (Fig. 4.2c, left panel), which is not consistently seen in vivo. In addition, we observed, in silico, an emergent loss of cohesion in CMT alignment between layers.

Indeed, CMTs along the outer and inner periclinal walls were oriented perpendicularly (Fig. 4.15). This peculiar effect, also not seen in vivo, probably results from the apical-basal growth of the outer wall (itself prescribed by the mediolateral cellulose orientation), which in turn is actively resisted by the stress responsive inner tissue.

To rule out this effect, we next performed simulations where mechanical feedback was A number of observations on katanin (ktn) mutants provided proof for such an uncoupling between outer and inner walls. KTN is involved in CMT alignment and its mutation leads to the formation of isotropic CMT arrays [START_REF] Bouquin | The Arabidopsis lue1 mutant defines a katanin p60 ortholog involved in hormonal control of microtubule orientation during cell growth[END_REF][START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF]. Different from wild type leaves, the CMTs on the outer periclinal membranes were systematically more isotropic in the ktn mutants bot and lue1 (Fig. 4.3a-c).

However, the degree of anisotropy of anticlinal CMTs was not affected by the mutation (Fig. 4.3d,e). This indicates that in the mutants the CMTs on the outer wall never align with the predicted stress patterns, in contrast to the inner, anticlinal walls. Consistent with our simulations, the mutant leaf and sepal blade were relatively wider, while maintaining a thickness at wild type levels (Fig. 4.3g, h). In summary, the results so far suggest a scenario, where CMTs systematically align along predicted stress patterns in internal walls. While this alignment guarantees leaf flatness, the degree of feedback on the outer cell wall is variable, and accounts for leaf blade width.

The previous simulations were initialized with a relatively flat ellipsoid. We also investigated the response of the system, starting from prolate spheroids only slightly perturbed in their degree of axisymmetry. This showed that the degree of flattening not only depended on the activation of the feedback itself, but also on the initial degree of shape asymmetry. Indeed, in silico, nearly spheroidal structures mainly grew in the apical-basal direction, resulting in the formation of elongated, finger-like shapes (simulation 5: Fig. 4. 2b-d). An axisymmetric structure would maintain itself as such, as the feedback mechanism on inner walls, is on its own not sufficient to break axisymmetry (Fig. 4.2b).

We next tested experimentally the predicted link between initial primordium flatness and final flatness in vivo, using the sepal, which is easily accessible for observation.

Leaf margin genes, such as WOX1 and PRS/WOX3, are expressed in the lateral and middle domains (Fig. 4.16a, b), and are essential for setting up initial asymmetry [START_REF] Guan | Spatial auxin signaling controls leaf flattening in arabidopsis[END_REF][START_REF] Nakata | Roles of the middle domain-specific wuschel-related homeobox genes in early development of leaves in arabidopsis[END_REF][START_REF] Nardmann | Symplesiomorphies in the wuschel clade suggest that the last common ancestor of seed plants contained at least four independent stem cell niches[END_REF]. Accordingly, plants with a double knock-out in both genes showed slightly narrower leaves and a clear reduction in sepal width (Fig. 4.3b and Fig. 4.16c,d). As the primordia of these mutants are still somewhat flattened, we would predict that reducing the initial mechanical feedback on the outer periclinal walls would partially rescue the narrow sepal phenotype. To this end, we introduced bot into a wox1 prs background, and found the width to thickness ratio of the sepals dramatically increased .

We tested a further reduction in primordium width/thickness ratio by combining wox1 prs with a third mutation, asymmetric leaves2 (as2), which has elongated and (nearly) axisymmetric sepals and leaves Nakata et al. (2012) (Fig. 4.4e, f). Although still set up as a slightly flattened structure, the primordia soon became almost axisymmetric.

They then mainly grew in the apical-basal direction (Fig. 4.3e and Fig. 4.17). The GFP-MBD marker revealed highly anisotropic CMT arrays on the anticlinal walls of sepals (not shown) and leaves (Fig. 4.4l). Importantly, along periclinal walls in wox1 prs as2, however, different arrangements were found, i.e. more isotropic in leaves and highly anisotropic in sepals (Fig. 4.4i and k). In other words, there is no specific CMT arrangement on the outer walls that correlates with the elongated shape of these organs. Therefore, axisymmetric or nearly axisymmetric shapes only seem to require CMT alignment along anticlinal walls. Consistent with theory, bot was not able to restore axisymmetric shapes in this background: sepals were shorter and thicker, but remained close to axisymmetric in the wox1 prs as2 bot quadruple mutant (Fig. 4.4g).

The results in Arabidopsis were further confirmed using microsurgery in tomato. Isolation of an incipient leaf primordium from the meristem resulted in compromised WOX expression, leaf margin formation, and flattening [START_REF] Shi | Model for the role of auxin polar transport in patterning of the leaf adaxial-abaxial axis[END_REF]. Both CMT arrays and cellulose microfibrils showed anticlinal arrangement in these axisymmetric leaves .

In conclusion, we have identified a conserved mechanism involving the coordinated behavior of the cytoskeleton in response to mechanical stress in internal, anticlinal walls. Although the precise mechanism behind this behavior remains elusive (Landrein and Hamant, 2013b), we suggest that a stress feedback has both the potential to amplify bilateral asymmetry during leaf development and to promote the elongation of stem-like organs, such as roots. This is a robust property, which is reproduced in our model with a minimum of hypotheses. In addition to the reinforcement of anticlinal walls along stress patterns, leaf flatness could in principle be further enhanced through cell division plane alignment in the same direction. It is known that cells often divide in a plane parallel to the microtubule interphase array (Ehrhardt andShaw, 2006, Wasteneys, 2002). Accordingly, we observed that division planes were mostly perpendicular to the plane of the leaf blade (Fig. 4.1k and m,Fig. 4.12 and Fig. 4.13a and b).

Such anticlinal walls should in principle further increase the resistance of the tissue to thickening, and thus both cellulose deposition and the orientation of new cross walls would contribute synergistically to the final leaf shape.

The amplification of asymmetry potentially provides a parsimonious explanation for leaf evolution. The widely accepted Zimmermann's telome theory proposes that a stem (telome) evolved into a thin leaf through series of shape transformations, which lack plausible molecular evidence [START_REF] Beerling | Zimmermann's telome theory of megaphyll leaf evolution: a molecular and cellular critique[END_REF]. According to our model, once asymmetry is established in a primordium, the CMT-mediated mechanical feedback would amplify the asymmetry to form a plenary leaf blade. Initial symmetry breaking can result from the asymmetric gene expression patterns at the shoot meristem, and likely involves asymmetric patterns of cell wall stiffness and expansion during early stages of development, in particular at the leaf margins [START_REF] Guan | Spatial auxin signaling controls leaf flattening in arabidopsis[END_REF][START_REF] Nakata | Roles of the middle domain-specific wuschel-related homeobox genes in early development of leaves in arabidopsis[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF]. (in contrast to F ). Compatibility is ensured by means of the elastic tensor F e , at the price of introducing mechanical stresses.

We take advantage of the fact that growth (that relies on cell physiological processes) occurs much slower than elastic relaxation. This allows us to treat the plastic and elastic regimes separately. Each step of the algorithm can be decomposed in two suboperations (detailed next), that consist of (i) computing the static equilibrium (Sec- Cartoon of the multiplicative decomposition of the total deformation. The current configuration is obtained by composing the growth tensor F g that defines the stress-free state of the system (that is generally incompatible) and the elastic tensor F e that introduces mechanical stresses.

4.A.2 Elastic regime

4.A.2.1 Elastic model

We locally model the wall as a thin membrane under plane stress, described through the following strain-energy function:

Ψ (E) = 1 2 E : C w : E (4.2)
where operator ':' is the tensor double-dot product; C w depicts the fourth order wall elasticity tensor; and

E def. = 1 2 F T e • F e -I (4.3)
is the Green-Lagrangian strain tensor associated with the elastic deformation (I denotes the second order identity tensor).

To model the fibrous structure of the wall, the elasticity tensor C w is expressed as a function of the angular density of cellulose microfibrils, as detailed in [START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF]. We describe the CMF distribution on each element using the local (π-periodic) density function ρ (in practice, we only store a truncated Fourier decomposition of this distribution). The wall elasticity tensor is detailed using a Cox-like fibrous model [START_REF] Cox | The elasticity and strength of paper and other fibrous materials[END_REF], expressing C w as a function of the low-frequency Fourier coefficients of the CMF angular density. These coefficients are here noted {ρ m } m∈N and {ρ m } m∈N , respectively representing the real and imaginary coefficients, and are given by:

∀m ∈ N : (ρ m , ρm ) def. = 2 (ℜ (ρ m ) , -ℑ (ρ m )) with ρm def. = π dθ π ρ (θ) e -2imθ . (4.4)
where θ is and angle parameter in the orthonormal material basis {i, j} associated with the triangle in its stress-free configuration (in which the stress, strain and elasticity tensors are expressed) The fiber-related stiffness is proportional to some fiber rigidity constant Y f . As originally done in [START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF], we also model the isotropic hydrogel embedding the fibers through an additional isotropic tensor, defined by its

Poisson ratio ν and a reduced Young modulus Y (i.e. the actual Young's modulus divided by 1ν 2 ). The wall elasticity tensor is expressed in Voigt notation as:

[C w ] = πY f ρ 0 16     3 + ρ 2 +4ρ 1 ρ 0 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 1 -ρ 2 ρ 0 3 + ρ 2 -4ρ 1 ρ 0 2ρ 1 -ρ 2 ρ 0 2ρ 1 +ρ 2 ρ 0 2ρ 1 -ρ 2 ρ 0 1 -ρ 2 ρ 0     Microfibril elasticity matrix + Y     1 ν 0 ν 1 0 0 0 1-ν 2     Hydrogel elasticity matrix . (4.5)
The values of Y and ν had little qualitative effect on the simulations (not shown). In the simulations presented here we used ν = 0.2 -which is comparable to the values used in [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF], [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] -and Y f ρ 0 ∼ 10 2 Y , meaning that the mechanical contribution of the soft hydrogel is negligible in our study.

4.A.2.2 Numerical resolution of the static equilibrium

To compute the equilibrium in practice, we use triangular membrane finite elements equipped with P 1 -Lagrangian shape functions [START_REF] Zienkiewicz | The finite element method: solid mechanics[END_REF]. Let N and T be the respective sets of the mesh nodes and triangular elements. We iteratively integrate the equation of motion on nodes until the pressure forces f ext and reaction forces f int balance. This step is performed using the backward Euler method and the conjugate gradient to solve the intermediate linear systems.

At each sub-step, the nodal reaction force f

(n)

int (for each node n ∈ N ) is evaluated by differentiating the total potential strain energy U of the system with respect to node position q (n) :

f (n) int = - ∂U ∂q (n) = - ∂ ∂q (n) τ ∈T ǫ (τ ) S (τ ) Ψ (τ ) = - τ ∈T ǫ (τ ) S (τ ) ∂Ψ (τ ) ∂q (n) (4.6)
where S (τ ) is the surface area of triangle τ ∈ T (reference configuration); ǫ (τ ) is the wall thickness at triangle τ (assumed to be constant in time, and either equal to ǫ in for inner walls or ǫ out = 3ǫ in for outer periclinal walls). The element-wise strain energy density (τ ) is obtained from Eq. (4.2):

Ψ
Ψ (τ ) def. = Ψ E (τ ) = 1 2 E (τ ) : C (τ )
w : E (τ ) . (4.7)

In practice, Eq. (4.6) is evaluated by splitting ∂Ψ (τ ) /∂q (n) in virtue of the chain rule:

∂Ψ (τ ) ∂q (n) = ∂Ψ (τ )
∂E (τ ) :

∂E (τ )

∂F

(τ ) e

: ∂F

(τ ) e ∂q (n) , (4.8) then using Eqs. 4.3 and 4.2 and the following expression:

F (τ ) e = n∈N q (n) ⊗ ∂N (n) ∂X τ (4.9)
where X depicts the 2D material coordinates in the grown configuration, and N (n) denotes the shape function associated with node n.

The external nodal loads (pressure forces) are given by:

f (n) ext = 1 3 τ ∈Tn s (τ ) P (τ ) n (τ ) (4.10)
where T n is the set of the finite elements that contain node n; s

(τ ) = S (τ ) det F (τ )
e is the surface area of triangle τ in the current configuration; n (τ ) is the normal to triangle τ ( n (τ ) = 1); P (τ ) is the difference in pressure between both sides of element τ (the sign depends on the orientation of n (τ ) ). We assume steady and uniform pressure P in within the tissue. Hence, P (τ ) vanishes if triangle τ does not belong to an outer periclinal wall, and is equal to P = P in -P 0 otherwise (P 0 being the outer atmospheric pressure).

4.A.3 Growth

4.A.3.1 Growth kinetics

To model growth we use the strain-based model developed in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], expressing the rate of growth tensor L g as a piece-wise linear function of E at equilibrium:

L g def. = ∂F g ∂t • F -1 g = Φ E -E y I + (4.11)
The previous expression provides a multidimensional extension of Lockhart's onedimensional model [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF]. Parameters Φ and E y respectively depict the wall extensibility and yield threshold; ( • ) + depicts the tensor ramp function. The latter is defined for any second order symmetric tensor T as:

(T ) + = k max (0, T k ) T k ⊗ T k (4.12)
where T k and T k are respectively the eingenvalues and corresponding normed eigenvectors of T .

4.A.3.2 Numerical resolution of growth

Starting from a system at static equilibrium at t, Eq. (4.11) is integrated, for each triangle τ , using the forward Euler method with constant time step ∆t:

F (τ ) g (t + ∆t) = I + ∆tL (τ ) g (t) • F (τ ) g (t) . (4.13)
Since immediately after this operation, the node positions have not moved yet, Eq. (4.1) can be rewritten as:

F (τ ) (t) = F (τ ) e (t + ∆t) • F (τ ) g (t + ∆t) , (4.14)
where F (τ ) e (t + ∆t) defines the new elastic deformation, that compensates for the new growth tensor. In virtue of Eqs. 4.1,4.13 and 4.14,F (τ ) e (t + ∆t) can be computed from the previous elastic deformation F g (t) is nonsingular). Eq. (4.15) expresses the idea that augmenting the plastic component amounts to relaxing some part of the elastic deformation by the same amount. This allows to simulate growth without explicitly storing F g , in so far as only the dilation factor I + ∆tL (τ ) g (t) is needed to evaluate Eq. (4.15).

Note the presence of both times t and t + ∆t in Eq. (4.15). This reflects the fact that, a priori, F (τ ) e (t + ∆t) does not satisfy mechanical equilibrium, hence the tilde notation. The deformation gradients F (τ ) e (t + ∆t) and the current positions give the initial conditions to compute the next mechanical equilibrium (Section 4.A.2.2).

4.A.3.3 Growth-induced advection of the microfibrils

In principle, growth induces a change in the fiber density and orientation, which can be assimilated to an advective process. In fact, by virtue of mass conservation, the volumetric density of CMFs may decrease as the wall expands (if we disregard cellulose intake). Moreover, anisotropic growth modifies the orientation of each fiber, which rotates towards the direction of main deformation.

To model these effects, we assume that the deformation of the fibers is affine, namely that these deform according to the macroscopic medium. We update the microfibril angular density incrementally, for each triangle τ and at each growth step t (N.B.: for the sake of clarity, we drop the t and τ notations from now on). This is done by replacing the old fiber distribution ρ with the new one (noted ρ), which verifies:

ρ (θ) ∝ (1 -α cos (2 (θ -ϑ))) -1 ρ (Θ) and π ρ = J -1 π ρ (4.16)
where J def.

= det I + ∆tL g depicts the incremental variation of volume; θ and Θ respectively parameterize fibers in the deformed (after incremental growth is applied) and initial (before incremental growth is applied) configurations; α and ϑ respectively measure the anisotropy and angle of the incremental growth.

In practice, this step is computed in the angle-domain, by first using the Inverse Fast Fourier Transform computed on the Fourier coefficients of ρ, which allows to compute Eq. (4.16); and then the Fast Fourier Transform, which provides the new Fourier coefficients of the microfibril distribution.

4.A.4 Stress feedback

To model the stress feedback, we express the dependency of cortical microtubules organization upon stress, and the CMT-guided cellulose deposition as detailed in Oliv- = cos θ i + sin θ j.

Eq. (4.17) captures a CMT accumulation in the first direction of stress (that maximizes the argument of the exponential). Moreover, the sharpness of this alignment is a growing function of both stress anisotropy and stress amplitude [START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF].

We model the polymerization/depolymerization of cellulose fibers through (i) a deposition of cellulose in direction θ proportional to φ (θ) (characterized by a kinetic constant k ρ ), coupled with (ii) a linear decay of cellulose (characterized by a kinetic constant k ′ ρ ):

.

ρ (θ) = k ρ φ (θ) -k ′ ρ ρ (θ) . ( 4
(θ) t=t 0 = k ρ πk ′ ρ , (4.19)
see [START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF]. After estimation of the first mechanical equilibrium, the CMF distribution is equated with the constant-stress steady solution of Eq. (4.18), that is proportional to φ. This allows for a coarse initialization of the CMF system at the very beginning of the simulation, which provides a more realistic initial condition than isotropic stiffness.

For better spatial regularity of the CMT density, we smooth the stress computed at each element over each cell interface, and use the smoothed expression of the stress in Eq. (4.17). Since the PK2 stress is given in the material coordinate system of each triangle, i.e. in the incompatible stress-free configuration, it cannot be directly used to perform a smoothing. In practice, this is performed by computing the stress in the current configuration of the material (Cauchy stress) from the PK2 stress S. Then, we compute a mean stress over each cell interface (weighted by triangle area). The outcome of this operation is generally of rank 3, as the interfaces may be slightly non planar. Therefore, the result is then re-projected onto the mesh triangles and pulled back to the reference configuration, providing a smoothed PK2 stress field.

4.A.5 Meshes

Each cell wall is composed of ∼10-20 triangular finite elements, and is common to two cells (unless it is an outer periclinal wall). Meshes are generated following the procedure detailed in [START_REF] Ali | Simulations and analysis of turgorinduced stress patterns in multi-layered plant tissues[END_REF], itself based on the open source algorithm DRACO-STEM [START_REF] Cerutti | DRACO-STEM: An Automatic Tool to Generate High-Quality 3d Meshes of Shoot Apical Meristem Tissue at Cell Resolution[END_REF], which allows to build meshes dedicated to finite element modeling, from incoming cell segmentations, which are 3D images. of each region coincide. It is here generated using Lloyd's algorithm [START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF].

The main benefit of a CVT, with respect to a more general Voronoi diagram, is that it allows more uniform cell volumes and shapes. This is a convenient property for us, as it allows to rule out the possible effects of cell shape and size variability on global morphogenesis. The CVT also builds cell layers that have a visually uniform thickness, which qualitatively resemble the actual cell layers seen in leaves.

4.A.6 Main parameters

In the simulations shown in the main text, we used the following parameters. For turgor pressure, we chose P = 0.25 MPa, which is of the order of magnitude of the values reported in [START_REF] Beauzamy | Mechanically, the Shoot Apical Meristem of Arabidopsis Behaves like a Shell Inflated by a Pressure of About 1 MPa[END_REF]. Meshes where scaled so that the average cell volume was equal to V cell ≈ 5 × 5 × 5 µm 3 , which is close to what was measured in our experiments (each mesh contains 800 cells, total volume is 800V cell ∼ 10 5 µm 3 ). The mean wall thickness was parameterized as (ǫ out + ǫ in ) /2 = 0.25 µm (with ǫ out = 3ǫ in ), which is close to the thickness measured experimentally in this work. The fiber rigidity is characterized by the Young modulus Y f k ρ /k ′ ρ ∼ 10 2 MPa [START_REF] Chanliaud | Mechanical properties of primary plant cell wall analogues[END_REF]. The growth characteristic time was set as Φ -1 ≈ 20 min. For a strain of 5%, this roughly corresponds to a expansion speed of about 15% per hour (cell length doubles every 5h), which is close to the parameter used in [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF].

Other parameters, which are less easily measurable, where set in a qualitative manner, in order to obtain a coherent simulation behavior: growth threshold E y = 1%, feedback

4.A.7 Softwares

In-house Python code was implemented in order to simulate growth and stress feedback. To compute the static equilibrium, we employed the finite element method as implemented in the open source software Sofa [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF] The meshing algorithm (Section 4.A.5) was initialized using meshes generated with the open source software Blender (www.blender.org). The DRACO-STEM tool [START_REF] Cerutti | DRACO-STEM: An Automatic Tool to Generate High-Quality 3d Meshes of Shoot Apical Meristem Tissue at Cell Resolution[END_REF]) is available at gitlab.inria.fr/mosaic/draco stem.

4.A.8 Sensitivity to initial position and coefficient γ

4.B Shape diagram

4.B.1 Construction

In order to monitor the change of a given growing shape, we track its principal dimensions over time, namely its length (λ 1 ), width (λ 2 ) and thickness (λ 3 ). In this work, we mostly focus on the various aspect ratios of the shape, and we disregard the volume of the structure, which, by contrast, characterizes the size of the organ. Hence, we can represent the shape at each time point through a normalized vector Λ defined as: In this model, we assumed that CMTs were stabilized in the direction of main stress, which in average results in a directional selection, promoting CMT coalignment in the axis of main stress (Eq. (2.9)). To represent the dynamic modification of the wall elasticity, we expressed the deposition of CMFs by transmembrane cellulose synthases as a function of CMTs distribution. One of the main assets of our approach, compared to other published works, is that the dynamic evolution of the wall elasticity explicitly depends on the kinetic constants of the chemical processes involved in wall turnover.

Λ i = λ i j λ 2 j . ( 4
In addition, coupling chemical kinetics with a growth-induced advection of the fiber raises a number of questions (Chapter 3). In particular, we discussed in Chapter 3 some minimal conditions on cellulose kinetic processes, allowing for a cell to maintain a homeostatic stiffness, while keeping the ability to perform a rapid change in growth direction if required.

As a limit of our approach, the affine deformation hypothesis (which is a prerequisite for the constitutive behavior developed in Chapters 2 and 3) will by construction not capture the possible heterogeneities of strain that may appear within the wall at a small scale. These could emerge when the dynamic evolution of the wall structure occurs in a time scale comparable to that of growth, resulting in differences in strain between the wall constituents. These may confer interesting properties to the mechanics of the wall, as for example modeled by Dyson and colleagues in the context of hemicellulose synthesis and breakage [START_REF] Dyson | A model of crosslink kinetics in the expanding plant cell wall: Yield stress and enzyme action[END_REF]. One may also remark that simple, linear strainstress constitutive relations (like that used in this thesis), may also lack realism. In fact, more nonlinear elastic effects are frequent in biological mechanical systems [START_REF] Chagnon | Hyperelastic energy densities for soft biological tissues: a review[END_REF]. These effects are not fully characterized in the context of the meristem, which is not easily accessible to mechanical testing, although Kierzkowski et al. have recently proposed, on an experimental ground, that nonlinear effects may occur at the SAM, in the form of a strain hardening effect [START_REF] Kierzkowski | Elastic domains regulate growth and organogenesis in the plant shoot apical meristem[END_REF].

This thesis is an effort to describe the link between wall structure and wall elasticity. Yet, the plastic regime was described in a phenomenological manner, inspired by Lockhart's approach [START_REF] Lockhart | An analysis of irreversible plant cell elongation[END_REF], which is itself an empirical rheological formulation. Finer growth laws were proposed to reflect the dynamics of wall turnover and irreversible creep. For example, Dyson and colleagues established a mechanistic explanation of Lockhart's law, where the viscous wall creep is the result of an imbalance between fast wall expansion (causing hemicellulose breakage) and slow wall turnover [START_REF] Dyson | A model of crosslink kinetics in the expanding plant cell wall: Yield stress and enzyme action[END_REF]. This model was later adapted to 2D growth in silico in [START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF].

Stress feedback in multi-layered tissues

Often enough, theoretical and experimental studies on meristem mechanics have focused either on the outermost cell layer or on the outer cell wall. In fact, according to the prevailing hypothesis, these bear most part of the tissue stresses (at least in the SAM, see Beauzamy et al., 2015, and references therein). This allows to assimilate the whole tissue to a pressurized shell, neglecting the inner walls of the tissue (as in Chapters 2 and 3 and [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Bozorg | Models of Mechanics and Growth in Developmental Biology: A Computational Morphodynamics approach[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF][START_REF] Kierzkowski | Elastic domains regulate growth and organogenesis in the plant shoot apical meristem[END_REF] for example).

The actual mechanical role of these inner walls is yet not straightforward. For instance, it was suggested, by simulation, that the presence of inner walls can yield stress patterns different from those predicted through a pressurized-shell hypothesis [START_REF] Ali | Simulations and analysis of turgorinduced stress patterns in multi-layered plant tissues[END_REF]. We also stressed that most of the shell-based models, coupling growth and stress feedback, in the literature [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] had been used in the context of prolate axisymmetric organs (typically stems). In modeling flat structures like leaves (which have a different type of symmetry), it appeared that a shell-based approach is probably not appropriate to model planar expansion (Chapters 3 and 4).

To represent the latter, we needed a finer representation of the tissue, taking the inner walls into account (Chapter 4). In fact, we suggest that oriented CMFs at inner walls are crucial for this type of growth regime. This would mean that a shell-based approximation is not only coarse, but fundamentally insufficient to represent planar expansion.

To model the expansion of this kind of extremely asymmetric structures, growth and stress feedback were simulated using ellipsoidal meshes, containing both outer periclinal wall and inner tissues, and composed of 800 cells. To generate these meshes, a dedicated method (presented in Chapter 4) was implemented, based on [START_REF] Cerutti | DRACO-STEM: An Automatic Tool to Generate High-Quality 3d Meshes of Shoot Apical Meristem Tissue at Cell Resolution[END_REF].

This full 3D approach allowed us to predict fine properties emerging from the coupling of growth and stress feedback in 3D alveolar tissues. We first observed that, in silico, a complete "activation" of the feedback in all membranes of the tissues, over-constrained the system, which resulted in a loss of cohesion in orientation between inner and outer CMTs/CMFs (see Chapter 4 and Fig. 4.15). This can be explained by the anisotropy of the outer wall that promotes apico-basal growth, which in turn stretches the inner cells.

These dynamically respond to this traction by orienting their CMTs accordingly. The anisotropy of the surface walls was not consistently observed in vivo, as in leaves and sepals, the outer periclinal microtubules rapidly became disorganized. If our modeling assumptions are correct, in particular if a stress feedback does exist in leaves, this could indicate that a fine regulation of the system's sensitivity is required, preventing 123 the effect predicted by simulation to occur in reality. Disabling the stress feedback on the outer wall of the tissues in our simulations (which is consistent with in vivo observations) in fact further allowed to generate more stable, flat growing shapes. This potentially provides a mechanistic explanation of anisotropic leaf growth.

Cell division and meshes

Another crucial aspect of development is cell division. This aspect was not addressed in this thesis. While simulating division in discrete 2D approaches (e.g. spring-based models) is trivial, no satisfactory computational solution has so far been proposed to perform this operation in cell-scale 3D models like those used in Chapter 4. In fact, simulating a cell division means adding a new wall within a cell and arbitrarily cutting the already existing walls in two. Performing this cut directly on the finite elements forming a given cell wall would probably result in a dramatic decrease of mesh quality.

Due to this limitation, the 3D simulations of growth presented in Chapter 4 do not integrate division and, for this reason, can only capture short developmental times.

Correlatively, anisotropic growth results in distortions of meshes that cause the FEM to become increasingly unreliable in the long run. A major future challenge, which can be a subject of research per se, is to design algorithms to dynamically update meshes and simulate cell division in an efficient manner. The aim here would be to sustain sufficient mesh quality in time, while ensuring the continuity of the various properties associated with the mesh triangles (in particular the strain). Remeshing/division methods tailored for planar [START_REF] Fozard | Vertex-element models for anisotropic growth of elongated plant organs[END_REF] or 2D curved [START_REF] Bozorg | A continuous growth model for plant tissue[END_REF] cells can be of interest while addressing those various issues with 3D polyhedral cells.

Do plant cells "measure" stress?

More than a decade ago, Hamant and coworkers reexamined the seminal idea (originally proposed by Green and King, 1966, in Nitella) of a causal relationship between mechanics and tissue physiology, mediated via microtubules [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF]. The appealing idea that forces may provide a positional and directional cue to cells has become increasingly popular in the plant biology literature and has motivated various experimental (e.g. [START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF], Hervieux et al., 2016[START_REF] Hamant | Mechano-devo[END_REF], Sampathkumar et al., 2014a[START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF][START_REF] Uyttewaal | Mechanical Stress Acts via Katanin to Amplify Differences in Growth Rate between Adjacent Cells in Arabidopsis[END_REF] or theoretical [START_REF] Alim | Regulatory role of cell division rules on tissue growth heterogeneity[END_REF][START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF][START_REF] Sapala | Why plants make puzzle cells, and how their shape emerges[END_REF] investigations. To date however, the hypothesis of a mechanical control of plant morphogenesis still lacks a clear mechanistic understanding (for reviews, see [START_REF] Hamant | Life behind the wall: sensing mechanical cues in plants[END_REF][START_REF] Hamant | Are microtubules tension sensors? Under review[END_REF], Landrein and Hamant, 2013a).

Despite the increasing experimental evidence underpinning the stress feedback hypothesis [START_REF] Hamant | Are microtubules tension sensors? Under review[END_REF], one may object that the original terminology employed in [START_REF] Hamant | Developmental Patterning by Mechanical Signals in Arabidopsis[END_REF] inter alia, proposing that "stresses orient microtubules" remains unsatisfactory from a physical standpoint. In fact, defining a force as a physical object on its own, other than the physical cause of acceleration, is not self-evident (see Feynman et al., 1965, chap. 12). Likewise, it is known that a force cannot be measured by other means than its effect (motion). That cells actively react to forces (and stresses) is hence a rather abstract statement, that does not easily allow to picture a reasonable physical and chemical mechanism. Indeed, evidence wants that the "perception" of a force should in any way rely upon the perception of a strain.

Could that strain be the tangential elastic strain of the walls? In anisotropic materials like the plant cell wall, the respective main axes of stress and strain may not coincide.

Instead they are perpendicular in a number of cases, as explored in Chapter 2. Thereby, it is equally reasonable to suggest that microtubules align perpendicular to main strain, rather than along main stress [START_REF] Schopfer | Mechanical signals in arabidopsis: Stress or strain[END_REF]. Nonetheless, this scenario fails in the rare cases where main stress and strain coincide, like in the concave crease around organs, where CMTs, main strain and main stress are supposed to be parallel [START_REF] Ali | Simulations and analysis of turgorinduced stress patterns in multi-layered plant tissues[END_REF][START_REF] Burian | A correlative microscopy approach relates microtubule behaviour, local organ geometry, and cell growth at the Arabidopsis shoot apical meristem[END_REF][START_REF] Dumais | Analysis of surface growth in shoot apices[END_REF]. Theoretical modeling also tends to disprove such a possibility, by suggesting that an alignment of CMTs/CMFs perpendicular to strain would be unstable [START_REF] Bozorg | Stress and Strain Provide Positional and Directional Cues in Development[END_REF].

One could also imagine that CMTs could be anchored to the wall and put under tension by wall deformation. As suggested by [START_REF] Inoue | Sensing surface mechanical deformation using active probes driven by motor proteins[END_REF] for example, CMTs under tension would in turn stabilize (as also discussed in [START_REF] Hamant | Are microtubules tension sensors? Under review[END_REF]. However, according to this scenario, CMTs would always align with the wall main strain. This would not be compatible with cases where main stress and strain are perpendicular and where CMTs are aligned along the main stress direction. Alternatively, CMT organization could be controlled by means of other intermediate molecular actors. CLASP proteins for instance, that accumulate at cell edges, regulate CMT alignment, by promoting CMT stability [START_REF] Ambrose | A CLASPmodulated cell edge barrier mechanism drives cell-wide cortical microtubule organization in Arabidopsis[END_REF]. A correlation between predicted stress and CLASP -or any other chemical actor -distribution, which has not yet been explored to my knowledge, could potentially provide insight into the problem.

The scenario of an indirect chemical regulation of CMT alignment would require that some information on wall stress would be chemically transferred to CMTs. This could be based on the tension-induced stabilization of some molecular actor embedded within the wall, acting as a stress sensor. However, how such a mechanism could provide a directional information is not self-evident. For example, this sensor could have an asymmetrice.g. filamentous -shape, resulting, through stabilization, in an average coalignment, and therefore a directional information. Yet, again, this tension should a priori be caused by wall strain (perpendicular to the wall stress): we find ourselves in the same predicament as before.

In this thesis, I have mostly tried to represent the emergent effects of a stress feedback, without questioning its actual physical foundations. It yet appeared to me that the various explanatory scenarios proposed in the literature are not completely satisfactory and are often abstract. In this conclusion, I would like to sketch a possible general mechanism, based on parsimonious physical principles and inspired by the mechanoperception mechanisms better identified in animals [START_REF] Charras | Tensile forces and mechanotransduction at cell-cell junctions[END_REF]. Briefly, the idea developed here is that the strain energy stored nonuniformly by the various adhesion zones around the cell, has the potential to provide an anisotropic, cell-scale chemical information on the external loading.

In the majority of plant cellular biomechanical models, it is more or less implicitly assumed that stresses are borne by the cell walls only, as cells are forced to be in perfect contact. This was for example an assumption of our model, where we assumed that two cells shared the same wall. In reality, each cell is embedded within its own wall exoskeleton, and contiguous cells are connected to each other through a sticky, pectinrich middle lamella [START_REF] Daher | How to let go: pectin and plant cell adhesion[END_REF][START_REF] Jarvis | Intercellular adhesion and cell separation in plants[END_REF][START_REF] Knox | Cell adhesion, cell separation and plant morphogenesis[END_REF]. This middle lamella potentially plays an important role in transferring mechanical forces from cell to cell, as experimentally investigated by Verger and coworkers, on Arabidopsis hypocotyl and cotyledon (Verger et al., 2018).

Consider a cell subject to an anisotropic -say tensile -external loading. A priori, the strain energy stored in the middle lamella should be heterogeneously distributed around this cell. In fact, zones of highest strain energy density should be those where cells tend to be split apart with the strongest pull, i.e. where the lamella is the most mechanically solicited. In a multicellular context the main axis of tissue-scale stress could be normal to these zones. This was suggested by Verger et al. who evidenced that in the adhesion-deficient mutant quasimodo1 [START_REF] Bouton | Quasimodo1 encodes a putative membranebound glycosyltransferase required for normal pectin synthesis and cell adhesion in arabidopsis[END_REF], cells detach along the lines of maximal predicted tension (Verger et al., 2018).

The heterogeneous distribution of adhesion strain energy provides a cell-scale directional information on the stress field. Mechanistically, this information could be transduced through a cascade of chemical reactions (that would remain to be characterized), from an adhesion strain sensor, all the way to CMT polymerization/depolymerization processes. CMT nucleation or stabilization could be for example promoted by specific chemical actors nonuniformly distributed around a cell (Fig. 5.1). In addition, CMT behavior is known to depend on the mutual interactions between CMTs [START_REF] Dixit | Encounters between Dynamic Cortical Microtubules Promote Ordering of the Cortical Array through Angle-Dependent Modifications of Microtubule Behavior[END_REF], which may amplify their coalignment.

This scenario, summarized in Fig. 5.1 could provide a conceptual basis for a mechanistic explanation of the stress feedback. In animals, adhesion-based mechanosensing -mediated via cadherins for instance -was better characterized (see [START_REF] Charras | Tensile forces and mechanotransduction at cell-cell junctions[END_REF]. In plants, the above-detailed mechanism could potentially account for the alignment of periclinal CMTs. Yet, it remains unclear whether it will be the case for other membranes. For instance, anticlinal membranes at the SAM often display CMTs pointing towards the outside, irrespective of how anisotropically distributed periclinal CMTs are. Further knowledge on the mechanical state of the middle lamella between the epidermis and the second cell layer is required. It also appears that this mechanism could provide a cell-scale information. Whether it can account for subcellular CMT patterns like those reported in Sampathkumar et al. (2014a) is unclear.

To my knowledge, this general scenario has not yet been explored, neither experimentally nor theoretically. Modeling could provide a tool to assess the theoretical plausibility of such a mechanism, from a macroscopic standpoint, for example by first allowing for the quantitative estimation of the adhesion-related strain energy distribution. A cell subject to an anisotropic loading would store adhesion-related strain energy in a differential manner, depending on the orientation of the membranes with respect to loading axis. In turn, the various levels of strain energy stored in the middle lamella may provide a directional information allowing to regulate CMT alignment. (a) schematic cylindrical organ (stress is mainly circumferential), close-up on the middle lamella containing mechanical sensors. (b) Cartoon of a cell aligning its microtubules along the axis of main tensile loading. The strain energy stored in the middle lamella is non uniformly distributed around the cell. This in turn modulates some chemical process (for example stabilization of some molecular actor) that provides the directional information.

Glossary

  Figure 1.1:(a) Cartoon of Arabidopsis thaliana's inflorescence (Image: www.pflanzenforschung.de) and close-up on the inflorescence meristem (scanning electron microscopy, scale bar = 100 µm, Image courtesy of J. Traas). (b) General structure of the shoot apical meristem (top view). Organs form at the periphery of the meristem at a regular pace, as the meristem tip progresses. (c) General structure of the animal and plant cells. Unlike animal cells, plant cells are embedded in a stiff wall. This cell wall prevents relative cell displacements, as two cells cannot detach from each other. By contrast, animal cells, that have their membranes in contact, may slide against one another more freely. Redrawn from[START_REF] Coen | The art of genes: How organisms make themselves[END_REF].

Figure 1

 1 Figure1.2: Within a multicellular context, the wall is subject to both cell-autonomous pressure-induced forces and cell-cell contact forces.

  Figure 1.3: (a) Instantaneous velocity at position x. (b)Growth of a 1D rod. Between t and t + dt, the infinitesimal portion of length dx as moved and increased in length. The new length is obtained by substracting the position of its right extremity x + dx + v (x + dx) dt and that of its left extremity x + v (x) dt. The line dilation rate of the element is equal to the ratio of the new length to the previous length, or 1 + L (x) dt. (c) Schematic representations of the various kind of deformations applied to an infinitesimal square volume. Pure strain, pure rotation, and general deformation.

  Figure1.6: (a) Strain-independent and (b) strain-dependent growth of a beam. Red arrows indicate the force that is required to maintain the extremity of the beam at a constant position. Right plots indicate associated strain energy profiles (qualitative). In (a), the growth may favor the accumulation of compressive stress, that eventually may result in mechanical instabilities. The beam may buckle under compression, reaching a more energetically-favorable state. This does not occur in case (b), as growth stops in absence of strain.

  Figure 1.7: (a) Piece-wise linear regime corresponding to Eq. (1.11): growth rate vs stress. (b) Ellipse-based visualization of the growth multidimensional threshold in three cases. From left to right: all the eigenvalues of E are below the threshold (no growth occurs); one eigenvalue only exceeds the threshold (uniaxial growth); all eigenvalues exceeds the threshold (two-dimensional growth).

  Figure 1.9: Control of growth anisotropy through wall stiffness (a-b) Atomic force microscopy scans of (a) disorganized and (b) organized cellulose microfibrils of onion cell wall. Adapted from Zhang et al. (2014) by permission from Springer. (c) The cell wall as a regulator of growth direction. Schematic of two pressurized cells with different wall properties and same initial shape. The first one has an isotropic wall and displays no privileged axis of growth. The other one displays hoop-shaped organization of microfibrils and expands axially. Reprinted from Oliveri et al. (2019).

  cortical microtubules (CMTs), as CSCs bind with CMTs (see Fig. 1.10(b) and Bringmann

  Figure 1.10: (a) Cartoon of the microtubule general structure. Redrawn from T. Splettstößer, Wikimedia Commons. (b) Cartoon showing the microtubule-guided transport of cellulose synthase complexes CSCs. The CSCs move along cortical microtubules and across the plasma membrane, and synthesize microfibrils in their wake. Reprinted from Oliveri et al. (2019).

Figure 1

 1 Figure 1.11: Hypothetical feedback loop governing plant morphogenesis: stress feeds back on CMTs, and regulates growth.

Foreword

  This is a post-peer-review, pre-copyedit version of an article published in Journal of Mathematical Biology (here referred as[START_REF] Oliveri | Regulation of plant cell wall stiffness by mechanical stress: a mesoscale physical model[END_REF]. The final authenticated version is available online at doi.org/10.1007/s00285-018-1286-y. Several minor corrections have been added to this chapter with respect to the published version.

  Figure 2.1: (a) The cell wall as a regulator of growth direction. Schematic of two pressurized cells with different wall properties and same initial shape. The first one has an isotropic wall and display no privileged axis of growth. The other one displays hoopshaped organization of CMFs and expands axially. (b) Cartoon showing the CMTguided transport of cellulose synthase complexes (CSCs). CSCs move along cortical microtubules and across the plasma membrane, and synthesize microfibrils in their wake. Redrawn after Landrein and Hamant (2013a).

  Figure 2.2: Mechanical model of the epidermis. (a) The tissue is represented by a pressure-filled thin membrane. (b) Sketch of the wall structural model: CMFs are immersed in an elastic solid isotropic hydrogel.

  Figure 2.3: Steady state of the CMT distribution. Microtubules respond nonlinearly to stress anisotropy and amplitude. Microtubule anisotropy depends equally on both stress anisotropy and amplitude. (a) Visualization of the CMT distribution φ according to the stress anisotropy α f (f 0 = 1 MPa) and the coupling coefficient γ. The red dashed line represents the polar distribution of forces f . (b) Anisotropy α ∞ φ of the CMT distribution as a function of f 0 (expressed in MPa) and stress anisotropy α f with stress feedback. (c) Relative concentration of polymerized tubulin c ∞ p /c 0 = π /2φ ∞ 0 /c 0 as a function of f 0 (expressed in MPa) and stress anisotropy α f with stress feedback. Parameters for (b-c): c 0 = 1 mol.µm -3 , γ = 10MPa -1 , K φ = 1.

  Figure2.4: The effect of stress feedback on elastic strain. Stress feedback confers nonlinear elastic properties to the material. Two noticeable effects are the effect of CMF reorganization on areal expansion (that may differ relative to the isotropic linear case) and the strain and stress respective directions (that may be orthogonal). (a) Areal dilation (Tr E ∞ ) ∆ OFF without (γ = 0) and ∆ ON with (γ = 1, unit: 10 -6 Pa -1 ) feedback. In absence of feedback, dilation is a linear function of f 0 , which is not the case with feedback. (b) Relative variation (∆ ON -∆ OFF ) /∆ OFF (in %) of areal dilation between a stress-responsive material (γ = 1) and a non-responsive material (γ = 0), as a function of f 0 (expressed in MPa) and stress anisotropy α f . Blue bold contour indicates level 0 (∆ OFF = ∆ ON ). (c) Ratio E ∞ ⊥ /E ∞ of strain (in the stress eigenbasis, where the basis vectors are sorted according to the descending order of eigenvalues) as a function of f 0 (expressed in MPa) and α f , with stress feedback. Main stress and strain are orthogonal iff E ∞ ⊥ /E ∞ > 1; domains labeled ⊥ and respectively indicate those where main stress and strain are orthogonal or collinear. Blue and red bold contours respectively indicate levels 0 (null orthogonal strain E ∞ ⊥ ) and 1 (isotropic strain E ∞ ⊥ = E ∞ ). Parameters: c 0 (arbitrary), γ = 10 MPa -1 , K φ = 1 (arbitrary), K ρ = 1 mol.µm -3 , Y f = 2 GPa, ν = 0.2, Y = 0.2 GPa.

  Figure 2.5: Virtual stem. Circumferential hoop stresses promote circumferential coalignement of CMTs resulting in axial elongation. (a-b) Stress field. Heat map: mean stress f 0 /2 (in MPa), elliptic representation of the tensor. (a) In absence of feedback, (b) with feedback. (c-d) Microtubule distribution. Heat map: CMT anisotropy α ∞ φ , ellipses: alignment tensor (see Appendix 2.E). (c) In absence of feedback, (d) with feedback. (e-f) Strain field. Heat map: dilation Tr E. (e) In absence of feedback, (f) with feedback. Parameters: stem dimensions 300 µm× 200 µm, turgor pressure P = 0.5 MPa, c 0 (arbitrary), γ = 2 × 10 -7 Pa -1 , K φ = 1, K ρ = 1 mol.µm -3 , Y f = 20 GPa, ν = 0.2, Y = 1 GPa.

  Figure 2.6: Virtual outgrowing organ. Highly anisotropic stresses at the concave frontier contribute to maintain the latter that expands mainly axially. (a) Original mesh and visualized cut (in blue). (b-c) Stress field. Heat map: mean stress f 0 /2 (in MPa), elliptic representation of the tensor. (b) In absence of feedback, (c) with feedback. (d-e) Microtubule distribution. Heat map: CMT anisotropy α ∞ φ , ellipses: alignment tensor (see Appendix 2.E). (d) In absence of feedback, (e) with feedback. (f-g) Strain field. Heat map: dilation Tr E. (f) In absence of feedback, (g) with feedback. Parameters: mesh characteristic dimension 100 µm (order of magnitude), turgor pressure P = 0.5 MPa, c 0 (arbitrary), γ = 10 MPa -1 , K φ = 1, K ρ = 1 mol.µm -3 , Y f = 20 GPa, ν = 0.2, Y = 0.3 GPa.

  Figure 2.7: Holed cap. Circular reorganization of CMTs allows to limit the effect of the breach by promoting radial expansion around the hole. (a) Original mesh and visualized cut (in blue). (b-c) Stress field. Heat map: mean stress f 0 /2 (in MPa), elliptic representation of the tensor. (b) In absence of feedback, (c) with feedback. (d-e) Microtubule distribution. Heat map: CMT anisotropy α ∞ φ , ellipses: alignment tensor (see Appendix 2.E). (d) In absence of feedback, (e) with feedback. (f-g) Strain field. Heat map: dilation Tr E. (f) In absence of feedback, (g) with feedback. Parameters: mesh characteristic dimension 100 µm (order of magnitude), turgor pressure P = 0.5 MPa, γ = 2 × 10 -7 Pa -1 , c 0 (arbitrary), K φ = 1, K ρ = 1 mol.µm -3 , Y f = 20 GPa, ν = 0.2, Y = 0.5 GPa.

  tween the definition of anisotropy based on Fourier coefficients (Eqs. (2.13) and (2.32)) and the nematic-based ratio employed in Boudaoud et al. (2014) (Eq. (2.12)). In fact, the matrix expression of Q φ becomes diagonal in the basis associated with the rotation that zeroes the odd Fourier coefficient φ1 (in this basis | φ1 | = φ1 = φ 1 /2 ∈ R + ). Since | φ1 | and φ 0 are base-invariant, the two eigenvalues of Q φ are equal to φ 0 ± 2| φ1 |, which gives:

Figure 3 . 1 :

 31 Figure 3.1: Deformation of an arbitrary angular distribution by an incremental deformation f g = I + ∆tL g .

  ∆t) by Eq. (3.1). The total deformation gradient can now be expressed according to the new irreversible deformation F

  Figure 3.2: (a) Schematic representation of the growth procedure applied to a given triangle element. (b) General simulation pipeline. (c, d) Surface meshes: (c) Cylinder closed with hemispherical caps. (d, e) Two partially merged spheres.

  Figure 3.3: Convergence of the growth method (a) Various meshes of a sphere with increasing resolution. (b) Error in radius (ε R ), strain (ε E ), CMF concentration (ε CMF ) and isotropy ε isotropy presented in the main text (Eq. (3.13)). Simulation parameters: turgor pressure P = 0.1 MPa, mesh typical length = 40 µm 3 , wall extensibility Φ = 1 s -1 , growth threshold E y = 2%, γ = 0, wall thickness ǫ = 0.4 µm, cellulose deposition constant k ρ = 1 mol.µm -3 .s -1 , cellulose decay constant k ′ ρ = 1 s -1 , surface fiber Young's modulus ǫY f = 40 MPa.µm and reduced surface Young's modulus of the gel ǫY ≈ 1 MPa.µm, matrix Poisson's ratio ν = 0.2.

Figure 3 .

 3 Figure 3.4 shows minimal examples of these behaviors, considering a given distribution subject to an arbitrary anisotropic and steady growth field. In a more general case, i.e. where growth depends itself on fiber distributions (Fig. 3.5), we use FEM simulations. Figure 3.5 shows simulations of growing pressurized membranes (Figs. 3.2(d) and 3.2(d)) with initial uniform and isotropic CMF distribution, where stress feedback, as well as cellulose turnover are inactivated (e.g. k ρ = k ′ ρ = γ = 0). Under this assumption, we verify that the total mass of cellulose is a conserved quantity (as shown by the dotted line in Fig. 3.5(b)). Moreover, we observe that the decrease of cellulose concentration is more rapid in fast growing zones (Figs. 3.5(c) and 3.5(d)), and that the CMFs tend to passively orient in the main direction of growth in the concave part of the mesh, where growth is anisotropic (Figs. 3.5(c) and 3.5(d)).

  Figure 3.4: Advection of a fiber distribution by a steady unidirectional growth. The arrows indicate the successive deformations of an initial orthonormal basis. Plain and dashed lines show the current and initial polar plots of the distribution. From top to bottom: initially isotropic distribution, anisotropic distribution orthogonal to the growth axis, general distribution.

  Figure 3.5: Coupling growth and CMF transport in a FEM simulation of growth. Growth simulation for mesh M3 (Fig. 3.2(e)). The simulation is initialized with uniform isotropic CMF density. No CMF polymerization/depolymerization are considered (k ρ = k ′ ρ = 0). (a) CMF concentration at various time steps. (b) Variation (value(t)/value(0) -1) of wall area, of average CMF concentration, of CMF total mass and of CMF average anisotropy. (c, d) Simulation result after several growth steps. (c) CMF density (as a % of the initial density) and anisotropy. Rods indicate direction of main alignment. (d) Growth rate and anisotropy. Rods indicate main growth axis. This simulation confirms the depletion of cellulose in fast growing zones, and the tendency of CMFs to rotate in the direction of main growth. Simulation parameters: turgor pressure P = 0.1 MPa, mesh typical length = 100 µm, wall extensibility Φ = 1 s -1 , growth threshold E y = 1%, γ = 0, reduced surface fiber Young's modulus ǫY f = 100 MPa.µm and reduced surface Young's modulus of the gel ǫY ≈ 1 MPa.µm, matrix Poisson's ratio ν = 0.2.

Figure 3 . 6 :

 36 Figure 3.6: Growth of a stress-adaptive sphere. Growth is computed until surface area is twice the initial one, for various values of the sensitivity coefficient γ and kinetic constant k ρ (taking K ρ = 1 mol.µm -3 ). Heatmap shows the areal dilation of triangles (namely the ratio between final area vs initial area). Histograms quantify mesh quality at the end of the simulation: skewness and roughness (Appendix 3.C) are respectively represented in orange and blue (vertical lines show median and maximum values). Simulation parameters: turgor pressure P = 0.1 MPa, mesh typical length = 40 µm, wall extensibility Φ = 1 s -1 , growth threshold E y = 1%, γ = 0, wall thickness ǫ = 0.4 µm, reduced surface fiber Young's modulus ǫY f = 40 MPa.µm and reduced surface Young's modulus of the gel ǫY ≈ 1 MPa.µm, matrix Poisson's ratio ν = 0.2. 79

  Figure 3.7: Growing sphere equipped with fast and sharp feedback: γ = 20 MPa -1 and k ρ = 5 mol.µm -3 , k ′ ρ = 5 (all other parameters are equal to those used in Fig. 3.6). Areal dilation (triangle current area divided by initial area, interpolated at nodes) shown in the initial configuration and discrete correlation length (see Appendix 3.D) for increasing mesh resolution.

Figure 3 . 8 :Figure 3 Figure 3

 3833 Figure 3.8: Stress at equilibrium in the oblate and prolate spheroids. Rods and colormap respectively indicate the main axis and trace of the Cauchy stress. Arrows indicate the axis of the shape elongation/flattening.

  Figure 3.11: Growth profiles of (a) the closed cylinder (Fig. 3.2(c)) and (b) merged spheres (Fig. 3.2(d)), for various values of γ. Each dot represents the distance from symmetry axis and height of the mesh nodes. Bold lines represent a polynomial fit. Red dots show the profile taken after the first equilibrium is reached. The other profiles correspond to the configuration when the mesh surface area has reached twice its original surface area. Simulation parameters: turgor pressure P = 0.1 MPa, mesh typical length = 40 µm, wall extensibility Φ = 1 s -1 , growth threshold E y = 2%, cellulose deposition constant k ρ = 1 mol.µm -3 .s -1 , cellulose decay constant k ′ ρ = 1 s -1 , reduced surface fiber Young's modulus ǫY f = 40 MPa.µm and reduced surface Young's modulus of the gel ǫY ≈ 1 MPa.µm, matrix Poisson's ratio ν = 0.2.

Figure 3

 3 Figure 3.12: (a) Topological distance and circles. b) Discrete autocorrelation function and exponential fit (example).

  this time active on inner walls only. This scenario follows our in vivo observations showing disorganized CMTs on the outer periclinal walls of growing leaves and sepals (Fig. 4.1h-j, Fig. 4.10b-d). Here, simulations led to further amplification of flatness (simulation 4: Fig. 4.2b-d), while the predicted arrangements of CMTs along both outer and inner walls were qualitatively in line with the in vivo observations (Fig. 4.2c, right panel). This scenario, involving an uncoupling between CMTs in inner and outer walls, seems therefore more plausible than the previous one. Note that, by contrast, activating mechanical feedback on outer walls only, resulted in reduced asymmetry as the virtual organ developed towards an axisymmetric elongated shape (simulation 3: Fig. 4.2b and d).

Figure 4 . 1 :

 41 Figure 4.1: Shape, CMT organization, and cell division in developing leaves and sepals. (a) Cross sections of Arabidopsis leaf primordia showing highly anisotropic growth. (b) CMT organization by immunostaining (green) with nuclei stained by DAPI (red) (c) CMFs stained by Direct Red 23 (white) in Arabidopsis rosette leaf primordia. (d) Overview of the same flower bud and cross section through the abaxial sepal at stages 3/4 and 5/6. (e-i) flowers sepals expressing GFP-MBD. (e) Overview of sepal at stage 3/4, inset indicates detail given in (f) showing anisotropic CMTs at the outer wall. (g) Same image stack as (f) but tilted to show highly anisotropic anticlinal CMTs (arrows). (h) Overview of sepal at stage 5/6, inset indicates detail given in (i) to show isotropic CMTs. (j) Quantification of CMTs on periclinal and anticlinal walls in sepals using Fibriltool (Boudaoud et al., 2014), showing differences in the degree of anisotropy along anticlinal (N=32 walls from 4 stage 3/4 sepals; N=52 walls from 5 stage 5/6 sepals) and periclinal (outer) (N=100 cells from 4 stage 3/4 sepals; N=207 cells from 5 stage 5/6 sepals) walls during sepal development. (k-m) Cell division pattern by mPS-PI staining in optical cross sections of tomato P3 treated with DMSO (k) or oryzalin (n) White, divisions perpendicular to the epidermis; blue, divisions parallel (angle < 30 • ) to mediolateral axis in inner cells or to the epidermis; green, other divisions (30 • ≤ angle ≤ 90 • ). Statistics are provided in (m). (n-p) Effect of oryzalin treatment on Arabidopsis sepal development after 24h (n) and 48h (o). (p) Quantification of width/thickness ratios. Treated sepals do not flatten (N=3 biological repeats). Scale bars, 20 µm in (a-c), 20 µm in (d),10 µm (e-i).

Figure 4 . 2 :

 42 Figure 4.2: Computational modeling. (a) Overview of the 3D model. Virtual tissues with outer and inner cells are composed of triangular FEs. At each time step, the element stiffness is updated in function of the stress tensor (Chapter 3). (b) Outcome of five different scenarios (result after 20 time steps). Starting from a flattened ellipsoid (simulation 1-4) different levels of flatness amplification can be achieved depending on whether the feedback is active (green dot) or inactive (red cross) on outer and/or inner walls. Without feedback (simulation 1) the structure becomes thicker, with feedback everywhere (simulation 2) the structure becomes slightly flatter, with feedback on the outside only (simulation 3) the structure will become thicker and move towards a cylindrical shape. Maximal flattening is obtained with feedback on inner walls only (simulation 4). A prolate spheroid (simulation 5) will remain axisymmetric with the same feedback on inner walls only. (c) Predicted CMT orientations in simulation 2 (upper panel) and 4 (lower panel). Both simulations predict anisotropic CMTs/CMFs along anticlinal walls as observed in vivo (see white line segments on cross-sections). Simulation 2 leads to highly anisotropic CMTs/CMFs on outer walls (white line segments on outer wall), which is not always observed in vivo. (d) Top diagram: ellipsoid shape changes. These can be represented as respectively points and trajectories on a 2D diagram (Appendix 4.B). Feedback in the inner tissues causes flattening (trajectories below the dotted line). Bottom diagram: in (near) axisymmetric structures elongation dominates (trajectories in green area). 99

Figure 4 . 3 :

 43 Figure 4.3: CMT arrangements in ktn mutant.(a) Overview of stage 3 bot flower bud expressing GFP-MBD showing more or less isotropic arrays. (b) Overview of stage 5 bot flower bud with isotropic CMTs. (c) Detail showing random orientation of periclinal CMTs in sepal in (b). (d) Tilted detail of (b). Anticlinal CMTs remain highly anisotropic in bot (indicated by arrows). (e) Quantification of microtubule orientations using Fibriltool. Wild type sepals at stage 3/4 have a higher degree of anisotropy on their outer periclinal walls than at stage 5/6, while in bot mutants anisotropy is low from early stages onwards. By contrast CMTs on anticlinal walls remain highly anisotropic throughout development. N=100 cells from 4 stage 3/4 sepals and N=255 cells from 7 stage 5/6 sepals for periclinal analysis; N=31 walls from 3 stage 3/4 sepals and N=36 walls from 4 stage 5/6 sepals for anticlinal CMT analysis. (f) Quantification of width/thickness ratios in Col-0 and bot sepals. N=10 WT and 7 bot sepals at stage 3/4. N=9 WT and 11 bot sepals at stage 5/6. N=11 WT and 8 bot sepals at stage 14/15. (g,h) Cross sections of mature leaves of Col-0 wild type (g) and lue1 (h). Scale bars, 10 µm in (a-d); 50 µm in (g, h).

Figure 4 . 4 :

 44 Figure 4.4: CMT and CMF arrangements in polarity mutants and effect of ktn. (a-c) Phenotypes of wild type (a), wox1 prs (b) and wox1 prs bot (c). (d) Quantification of width/thickness ratio, showing that bot increases width in the double mutant. Note that there is still some margin identity left. N=10 WT, 7 wox1 prs and 8 wox1 prs bot sepals at stage 3/4. N=9 WT, 9 wox1 prs and 8 wox1 prs bot sepals at stage 5/6. N=11 WT, 10 wox1 prs and 8 wox1 prs bot sepals at stage 14/15. (e) Radialized sepal primordia showing increased boundary domains. (f) Radialized organs (overview) in flower bud of a triple wox1 prs as2 mutant. (g) Quadruple wox1 prs as2 bot mutant organs remain close to axisymmetric. (h-l) Overview of phenotype and CMT alignment in a finger-like sepal (h-i) and leaf (j-k) of wox1 prs as2. (i, k) Details of (h) and (j) showing anisotropic (j) and random (k) CMT arrangements on periclinal walls, respectively. (l) tilted detail of (j) showing that CMTs in anticlinal directions. (m-p) Isolation of an incipient leaf primordium in tomato from the meristem results in the formation of radialized leaves (m, n). CMTs (o) and CMFs (p) are mostly oriented in anticlinal directions. Scale bars: 100µm in (a-c), 50 µm in (e), 100 µm in (f, g), 20 µm in (h), 10 µm in (i), 20 µm in (j-l), 100 µm in (m, n), 20 µm in (o, p).

Figure 4

 4 Figure4.5: Morphoelastic description of growth kinematics. Cartoon of the multiplicative decomposition of the total deformation. The current configuration is obtained by composing the growth tensor F g that defines the stress-free state of the system (that is generally incompatible) and the elastic tensor F e that introduces mechanical stresses.

  figuration of each triangle element, expressed in the 2D material coordinates); and θ def.

  .18) Eq. (4.18) is solved by the classical Runge-Kutta procedure (RK4) with constant time step size (the time step size ∆t for growth and cellulose deposition are the same), assuming constant stress during the time step. At the beginning of the simulation, CMT and CMF distributions are isotropic and uniform, taking the trivial null stress steady solution of Eq. (4.18) as an initial condition for ρ: ρ

  Fig. 4.6(a)). A CVT is by definition a Voronoi tessellation wherein the centroid and seed

  available at www.sofaframework.org. Data structures (cellularized meshes) are implemented in the Python library CellComplex (gitlab.inria.fr/mosaic/cellcomplex). Visualization of 3D simulations was based on the open source software TissueLab (github.com/VirtualPlants/tissuelab).

Fig. 4

 4 Fig.4.7 provides a sensitivity analysis of the model with respect to initial shape (i.e., initial position in the shape diagram) and feedback sensitivity γ, in the case where outer stress feedback is inhibited. This was performed by generating 7 different ellipsoidal meshes with equal volume and cell number, following the procedure detailed in Section 4.A.5. In this way, only the various aspect ratios varied between all the meshes. For null γ we robustly obtained trajectories heading towards a sphere (represented by blue trajectories in Fig.4.7), which corresponds to the intuitive scenario where all walls are isotropic in stiffness. Higher γ generally resulted in flattening shapes (red and green trajectories in Fig.4.7), for all the tested initial positions.

  Figure 4.8: Morphospace. (a) 3D visualization of a shape vector Λ. (b) 2D spherical projection of the triangle morphospace.

Figure 4 .

 4 Figure 4.10: CMT organization in developing Arabidopsis leaves. (a) Young leaf primordium expressing GFP-MBD, top view, showing anisotropic CMTs . (b-d) Slightly older leaf showing different CMT arrangements. Boxed areas detailed in (c) (showing cells with anisotropic CMT arrays) and (d) (cells with isotropic CMT arrays). (e) Tilted version of (c), showing anticlinally oriented CMTs (arrow). Scale bars, 10 µm.

Figure 4 . 11 :

 411 Figure 4.11: Organization of anticlinal CMTs and CMFs in tomato and Arabidopsis. (a-c) The organization of CMTs in the cross section of a tomato P3 by immunostaining with tubulin antibody. (a) overview, (b) shows the magnification of a part of (a), Arrow heads show CMTs. (c) The orientation of CMFs in the cross section of a tomato P3 stained by Direct Red 23 dye. (d-i) Live imaging of 3xYFP-labelled Cellulose Synthase Interacting 1 (CSI1/ POM2, green) in mCherry-TUA5 (pink) background sepals. (d-f) average projections to reveal CSI trajectories along microtubules on periclinal surface membrane. (g-i) show trajectories along anticlinal walls. Note anticlinal trajectories of CSI (arrows). Imaging by z-stacking (0.3 µm steps) and time-lapsing (1 min intervals) for 20 min. Scale bars, 20 µm in (a-c), 10 µm in (d-i).

Figure 4 . 12 :

 412 Figure 4.12: Division orientations in inner and epidermal cells of leaf primordia. (ac) Cell linage tracing analysis in Arabidopsis leaf primordia using a heat-shock induced Cre-loxP system. No endoplasmic reticulum localized GFP signal (GFPer) is available without heat shock (w/o HS) (a). In contrast, in leaf primordium 72 hours after heat shock (w/ HS), GFPer is observed in continuous cell files (white arrowheads) along the mediolateral direction (b). (c) shows the magnification of a part of (b). (d-e) The distribution of inner cell division angle which is against the mediolateral axis labeled as blue (d) in a collection of optical cross sections of one entire tomato P3 from the tip to the base (e). (f, g) The distribution of epidermal cell division angel which is against the corresponding tangent labeled as blue (f) in a collection of optical cross sections of one entire tomato P3 from the tip to the base (g). Scale bars, 20 µm in (a) and (b).

Figure 4 . 13 :

 413 Figure 4.13: Importance of CMTs for orientated cell division and planar leaf formation. (a, b) Cell division pattern by mPS-PI staining in optical cross sections of tomato P2 (a) and P3 (b) show that the inner cells tend to divide perpendicularly to the mediolateral axis (blue line), while the epidermal cells exclusively perform periclinal divisions with the division plane perpendicular to the outer surface. The newly formed cell walls are labeled green for inner cells and white for epidermal cells. The depth of the optical cross sections from the tip of corresponding leaf primordium is shown. (c) Cell division pattern by mPS-PI staining in optical cross sections of tomato P3 24h after the treatment of DMSO (Mock) (left column) and oryzalin (right column). The newly formed cell walls are labeled green for inner cells and white for epidermal cells. The depth of the optical cross sections from the tip of corresponding leaf primordium is shown. (d, e) The morphology of tomato P3 48h after the treatment of DMSO (d) or oryzalin (e). Arrows indicate the primordia treated with chemicals. Leaf primordium treated with DMSO can perform normal anisotropic growth and generate lateral leaflet primordia (asterisk). The anisotropic growth and planar leaf form is compromised in oryzalin treated samples. Scale bars, 20 µm in (a-c), 100 µm in (d, e).

Figure 4 . 14 :

 414 Figure 4.14: Cell wall thickness in tomato leaf primodia. (a) TEM micrograph of the epidermal cell lay in tomato P3 leaf primordium cross section showing the out cell wall (between blue arrowheads) is thicker than inner cell walls (between red arrowheads). (b) Quantification of (a). Scale bar, 2 µm in (a).

Figure 4 . 15 :

 415 Figure 4.15: Mechanical conflict generated by feedback on all walls. When all wallsshow mechanical feedback, the outer layer will mainly orient its CMTs perpendicular to the longest axis of the ellipsoid and cause the structure to lengthen. This will cause resistance of the inner cells, which will orient their CMT along that axis. This will generate distortions in the long run. Note that in the axial plane, all CMTs orient anticlinally as also observed in vivo. Green arrows indicate main CMT orientations in different planes and layers.

Figure 4 .

 4 Figure 4.16: Expression pattern of WOX genes and phenotype of wox1 prs in Arabidopsis. (a, b) Expression patterns of pWOX1:NLS-GFP3 (a) and pPRS:NLS-GFP3 (b) in optical cross sections of Arabidopsis leaf primordia. Both reporters show GFP expression in the middle domain, while pWOX1 also shows additional activity in several cell layers in leaf margins. (c, d) Mature leaf phenotypes of 17 DAG old Col (c) and wox1 prs (d). Scale bar, 20 µm.

Figure 4 . 17 :

 417 Figure 4.17: Boundary formation in wild-type and wox1 prs as2 mutant flowers. (a) Segmented 3D reconstruction of wild type (35S: GFP-Lti6b) at 0h, showing degree of Gaussian curvature. Blue color indicates negative curvature which is a marker for the boundary. (b) Confocal, 3D reconstruction showing the same flower bud at two time points. The site of negative curvature is marked blue. The width of this boundary does not change and remains about 10 µm wide. (c) Segmented 3D reconstruction of wox1 prs as2 (35S: GFP-Lti6b) at 0h. (d) Development of boundary in the triple mutant at three time points. The initial zone of negative curvature is slightly broader than the wild type and then gradually increase in size (from about 15 µm to 45 µm, marked with color from blue to yellow and red). Scale bars, 20 µm.

Figure 5 . 1 :

 51 Figure5.1: Sketch of a conceptual adhesion-based mechanism for the perception of tissue-scale stresses in plants. A cell subject to an anisotropic loading would store adhesion-related strain energy in a differential manner, depending on the orientation of the membranes with respect to loading axis. In turn, the various levels of strain energy stored in the middle lamella may provide a directional information allowing to regulate CMT alignment. (a) schematic cylindrical organ (stress is mainly circumferential), close-up on the middle lamella containing mechanical sensors. (b) Cartoon of a cell aligning its microtubules along the axis of main tensile loading. The strain energy stored in the middle lamella is non uniformly distributed around the cell. This in turn modulates some chemical process (for example stabilization of some molecular actor) that provides the directional information.

apoptosis

  Programmed cell death. . . . . . . . . . . . . . . . . . . 7 Cauchy stress Physical quantity that expresses the internal forces acting within an infinitesimal volume of matter in the current configuration. . . . 39 cell division Process by which a parent cell divides into two daughter cells. . . . . . . . . . . . . . . . . . . . . . . . 7 cell growth Irreversible expansion of the cells due to turgor-induced forces. . . . . . . . . . . . . . . . . . . . . 6 cell wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 CMFs cellulose microfibrils . . . . . . . . . . . . . . . . . . . . . 24 CMTs cortical microtubules . . . . . . . . . . . . . . . . . . . . 25 CSCs cellulose synthase complexes . . . . . . . . . . . . 25 deformation gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 elasticity tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23, 36 Eulerian strain rate tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 FEM finite element method . . . . . . . . . . . . . . . . . . . 15 Green-Lagrangian strain Geometric measure of how much a displacement locally differs from a rigid body displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22, 23 growth tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 meristem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 morphoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 morphogenesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 multiplicative decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 organogenesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 osmosis Spontaneous movement of water through the (semipermeable) cell membranes due to a difference in solute concentration between both sides of the membrane. This allows cells to absorb water. . . . . . . . . . . . . . . . . . . . . . . . . . . 6 rate of growth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 SAM shoot apical meristem . . . . . . . . . . . . . . . . . . . . . 5 second Piola-Kirchhoff stress Physical quantity that expresses the internal forces within an infinitesimal volume of matter in the reference configuration. . . . . . . . . . . . 23 spin tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 strain energy density Volumetric density of potential strain energy (in J.m -3 ) which is stored by a mechanical system undergoing a deformation. . . . . . . . 10 stress Physical quantity that expresses the internal force within an infinitesimal volume of matter. It can be mathematically expressed in various manners, in particular depending on whether we express the forces in the current or reference configuration of the material. 10 tissue A contiguous set of cells.. . . . . . . . . . . . . . . . . . 2 turgor pressure Cell internal hydrostatic pressure resulting from osmotic water intake, which is resisted by the cell wall. . . . . . . . . . . . . . . . . . . . . . . . . 2, 6 velocity gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

There are other types of meristems (root meristems, axillary meristems, cambium, etc.).

By contrast, in animals, cells are freer to slip against one another. Many animal developmental processes involve collective cell migration[START_REF] Friedl | Collective cell migration in morphogenesis, regeneration and cancer[END_REF] or apoptosis (e.g. formation of digits,[START_REF] Jacobson | Programmed cell death in animal development[END_REF].

In a certain sense, this amounts to looking for a balance of force, yet based on an explicit formulation of the total potential.

See for example en.wikipedia.org/wiki/List of finite element software packages.

This term was coined by Alan Turing in "The chemical basis of morphogenesis"[START_REF] Turing | The chemical basis of morphogenesis[END_REF].

Due to the symmetries C ijkl = C klij and C ijkl = C ijlk = C jikl = C jilk , a given elasticity tensor C can be expressed as a 3 × 3 matrix (Voigt notation).

The actual Young's modulus divided by 1 -ν 2 .

In reality, fibers are connected through the hemicellulose network, that is yet assumed to be significantly softer than the fibers and then negligible.

 4 For given positions x and angle θ (angle with ex in the tangential plane), ρ (θ) dθ is the number of fibers at the vicinity of x whose orientation angle lies in the infinitesimal interval [θ, θ + dθ

].5 A fiber can be modeled by a straight line, equivalently defined by the unit vectors e θ = cos θex + sin θey or -e θ .

Note that this definition makes sense only for tensile stress (S1 ≥ S2 ≥ 0).

This can be mathematically intuited by additionally assuming that c0 ≫ π φ, which reduces Eq. (2.9) to a linear ordinary differential equation. Its solution displays an initial slope (|dtφ (θ, t0) |) proportional to k ′ φ (S, θ) (that is a decreasing function of f0).

In the transient regime however, orthotropy cannot be assumed, since the choice of the reference axis does not eliminate both ρ1 and ρ2 in general.

The notion of symmetry is here to be understood in the sense used in physics, that is the existence of a transformation that leaves some property -here the shape -of an object unchanged.

Structurally symmetric around the axes of anisotropy.

Acknowledgements viii

Acknowledgements

This PhD work was performed in the French National Institute for Computer Science and Automation (Inria). Funding was provided by the Inria Project Lab Morphogenetics. Foremost, I want to express my immense gratitude to my supervisors Christophe Godin and Jan Traas, for their constant support, patience and advice. Christophe, I am grateful for your infinite scientific enthusiasm and for the moral and professional support that you gave me. I hope to become a scientist as curious, pedagogical and passionate as you are. Jan, I am thankful for all the stimulating scientific discussions, your wise advice, and your sense of humour! Working with you was a very pleasant and enriching experience, professionally, intellectually and humanly. I want to warmly thank the other members of the jury, Hervé Delingette and Yoël Forterre -for kindly accepting to review this manuscript -and Andrea Parmeggiani. I truly enjoyed the discussions we had during my defense, which was a very pleasant moment of science. I thank Olivier Ali, who closely followed my work and was an indispensable actor of this thesis. Thank you for your support and friendship, and, naturally, for all the hoppy scientific discussions! Many thanks also to Romain who needs to practice his chess openings, to my sicilian mate Bruno, to Guillaume, Bea and Ibrahim for the memorable musical sundays in Montpellier, Jean-Phi for your kindness and for your help, Jo, Eugenio, Emmanuel, Fred, Julia, Stéphaneverger, Thomas, Florian & Florian, Antoine, Julie, Feng, Leia, Alice, Hanna, Corentin, Svenja, Long, Shogo, Matthieu, for your friendship and support. To my dear friends Stan & Pauline, Thomas, Davide, Sezos ¸& Gizos ¸, Astrid, Lo, Pierrick, Sibel, and to my copains du foyer -Carlito, Lise, Léa, Laia, Jeff, Jérémy, Mathilde, Benoît, Arthur, Douglas, Raquel, Chie, Michiel, Fab, Jack, Aline, Thomas, Anaïs, Caro, Nico, Laurine, Léo, Yoan, and all others -who have become my second family in Lyon. Sincerely, thank you to all.

I would also like to thank the Inria team assistants, Laurence and Sylvie, for their great work and help, as well as all the rest of the RDP staff. I am grateful to my family -my mother, my father, my sister, my brothers, and mia zia Antonietta -for always being proud of me. Last but not least, to Morgane, for the infallible moral support and affection that you showed me, in the best and in the worst times, thank you so much.

Finally, our model raises further questions about the dynamics of microfibrils. The latter is described by Eq. (2.7) that contains a degradation term. This term is required mathematically for the existence of a steady regime. Without it, the cell wall would undergo ceaseless accumulation of cellulose. Whether such a limiting term has an enzymatic origin or could be explained by a growth-related dilution/advection process remains to be established. It is in fact not self-evident that growth-induced advection/diffusion could suffice for rapid changes in growth direction (as observed during organ initiation). Indeed, in principle, newly deposited microfibrils would be in mechanical competition with older layers with different orientation. This competition would further slow the mechanical responsiveness of the tissue down. In depth theoretical modeling should provide further insights into this problem.

Acknowledgments

The authors would like to thank Guillaume Cerutti for assistance with the visualization tool TissueLab (github.com/VirtualPlants/tissuelab). Funding was provided by Inria Project Lab Morphogenetics and European Research Council (Grant No. 294397). these walls (Fig. 4.1c and Fig. 4.11c), while the cellulose synthase-associated proteins followed anticlinal paths along the CMTs (Fig. 4.11d-i).

To further evaluate the role of the CMTs in leaf development, we treated primordia with the CMT-depolymerizing drug oryzalin, at concentrations where they continued to grow. After treatment, the width to thickness ratio did not increase, in contrast to the untreated controls . Outgrowing leaves and sepals were thicker, while lateral expansion was compromised. When cells continued to divide, division plane alignment became randomized. This shows that CMTs arrays are crucial for asymmetric leaf expansion.

How do these heterogeneous and dynamic CMT arrangements on anticlinal and periclinal membranes emerge? In plants, turgor pressure and differential growth both generate mechanical stresses within the cell walls [START_REF] Peters | The history of tissue tension[END_REF]. It has been proposed that these stresses serve as a regulatory cue for cellular growth [START_REF] Ali | Physical Models of Plant Development[END_REF], Sampathkumar et al., 2014a). Indeed, CMTs often align along the axis of maximal tension (Hamant et al., 2008, Landrein andHamant, 2013b). This in turn would lead to the CMT-guided deposition of cellulose microfibrils (Landrein and Hamant, 2013b) and wall reinforcement restricting growth in this orientation of maximal tension (Fig. 4.11di).

We first investigated if this so-called stress feedback mechanism could provide, on theoretical grounds, a plausible scenario for leaf morphogenesis. To this end, we developed a computational modeling approach. The models of leaf development proposed in the literature -(e.g. [START_REF] Hervieux | A mechanical feedback restricts sepal growth and shape in arabidopsis[END_REF][START_REF] Kuchen | Generation of leaf shape through early patterns of growth and tissue polarity[END_REF][START_REF] Qi | Mechanical regulation of organ asymmetry in leaves[END_REF] -are in 2 dimensions. Therefore, out-of-plane walls are not taken into account, although they could significantly impact the mechanics of the system. To alleviate this limitation, we developed a 3D finite element (FE), multicellular model (adapted from [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], Oliveri et al., 2019, see also Appendix 4.A), to analyze the effect of mechanical feedback at the level of the entire growing organ. The effect of cell division, not taken into account in these simulations, was considered to be negligible as the simulations were only carried out over short time periods.

Incipient leaves were represented as ellipsoidal alveolar structures (composed of 800 cells) under steady and uniform pressure (see Appendix 4.A). We used ellipsoids of initial aspect ratios comparable to that of young primordia (Fig. 4.2). To account for the differences in thickness between outer and inner walls observed in vivo (Fig. 4.14), the outer walls were made 3 times stiffer than the inner walls in the model. Inside these structures the dominant wall forces were in the dorsoventral direction (i.e. along shortest axis of the ellipsoid). By contrast, at the outer surface, wall strain and

Acknowledgments

The authors would like to thank Olivier Hamant for reviewing the manuscript, and Thomas Laux, Ari Pekka Mäh önen, Ben Scheres, ABRC and NASC for providing seeds. This work was funded by the National Natural Science Foundation of China grants Appendices

4.A 3D growth model Author's note

This section is of the supplementary materials associated with the present preprint.

It here mostly serves the self-containment of this chapter and contains a number of redundancies with respect to Chapters 2 and 3. However, several modeling specificities are here introduced, motivated by the use of more complex structures.

4.A.1 General ideas

The tissue is modeled as a multicellular alveolar structure, each cell being described as a set of connected walls, which together form a 2D continuum, constantly loaded with a steady and uniform pressure. To model both the elastic and plastic effects occurring in the virtual organ, we adapted the model detailed in [START_REF] Boudon | A Computational Framework for 3D Mechanical Modeling of Plant Morphogenesis with Cellular Resolution[END_REF], which represents the growing tissue as a morphoelastic system (Goriely andBen Amar, 2007, Rodriguez et al., 1994). This consists in postulating that the total deformation gradient F (from the initial to the current configuration) is the product of an irreversible (plastic) component F g and a reversible (elastic) component F e (the subscripts 'g' and 'e' respectively stand for growth and elastic):

(see Fig. 4.5). Constitutively, in the morphoelastic approach, we assume that the strainenergy density Ψ is a function of F e only. By contrast, the growth tensor F g defines a stress-free configuration. In so far as each individual region of the domain may grow independently from its neighboring, the growth tensor F g in general does not define a compatible configuration, i.e., it is not the gradient of a continuous displacement field sensitivity γ = 30 MPa -1 .µm -1 (Section 4.A.8 provides a sensitivity analysis), cellulose

.7: Sensitivity analysis with respect to initial shape and γ (in MPa -1 .µm -1 ), see Section 4.B. We used here k ′ ρ = 1 h -1 and Φ -1 ≈ 10 min.

In the sequel, we will note Λ = (L, W, T) the Λ-position of a given 3D shape. As Λ is of norm 1, it is equivalent to a point of the unit sphere (Fig. 4.8(a)). Hence, it can be equivalently represented using a longitude-latitude spherical parameterization (θ, ϕ). Intermediately, edges ac and bc capture all possible shapes having two equal dimensions (prolate for ac or oblate for bc). Last edge ab refers to all the degenerated shapes with null thickness.

4.C Supplementary figures