Tomáš Ligurský
email: tomas.ligursky@ugn.cas.cz

Yves Renard

A Method of Piecewise-Smooth Numerical Branching

Keywords: Numerical branching, piecewise smooth, steady-state problem, contact problem, Coulomb friction. MSC (2010) 65P30

A method of numerical branching is proposed for piecewise-smooth steady-state problems when any analytical expressions are not known for the regions of smoothness of the function involved. The performance of the method is shown for model bifurcations in discretised contact problems with Coulomb friction. First, different parameter settings in its input are tested so that the optimum one can be proposed. Then, the method is used to investigate the behaviour of the bifurcations for different meshes. All solution branches seem to be discovered reliably if the parameters are set up properly.

Introduction

The steady-state bifurcation problem:

Find y ∈ U such that H(y) = 0, (P)
where U ⊂ R N +1 and H : U → R N , has been the subject of large number of studies in the last decades. A variety of numerical bifurcation methods has been constructed if H is smooth, say continuously differentiable (see, e.g., [START_REF] Allgower | Numerical path following[END_REF][START_REF] Georg | Matrix-free numerical continuation and bifurcation[END_REF][START_REF] Hughes | A bisection-like algorithm for branch switching at a simple branch point[END_REF] and references therein). On the other hand, there are many equilibrium problems in economics and diverse engineering fields whose models lead to a system of non-smooth equations [START_REF] Alexander | Piecewise smooth homotopies[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]. Let us mention frictional and frictionless contact problems in solid mechanics, which are of our specific interest, for instance. However, methods of numerical branching of solutions of non-smooth problems depending on a parameter are still very little explored to our knowledge: Branching of static equilibrium curves of discretised frictionless contact problems were treated in [START_REF] Björkman | Path following and critical points for contact problems[END_REF][START_REF] Schulz | Equilibrium paths of mechanical systems with unilateral constraints I[END_REF], where the tangential directions of curves emanating from points of non-smoothness were characterised as the solutions of a certain linear complementarity problem, and a method based on resolution of this problem was suggested for branch switching during numerical continuation (see also the theoretical survey [START_REF] Klarbring | Stability and critical points in large displacement frictionless contact problems[END_REF]). Furthermore, the possible evolutions in quasi-static plane frictional contact problems were determined similarly in [START_REF] Da | The evolution and rate problems and the computation of all possible evolutions in quasistatic frictional contact[END_REF].

In our recent paper [START_REF] Ligurský | A continuation problem for computing solutions of discretised evolution problems with application to plane quasi-static contact problems with friction[END_REF], we developed a restarted predictor-corrector method for numerical continuation of solution curves of Problem (P) provided that H is piecewise C 1 (P C 1 ). The method consists in continuing smooth solution branches by a standard predictor-corrector method and joining the smooth branches continuously.

More precisely, our method computes a sequence of points {y k } lying approximately on a selected solution branch together with a sequence of the corresponding tangent vectors {t k } with a unit weighted norm t k w . Starting from a couple (y k , t k ), one predictor-corrector step yields a new couple (y k+1 , t k+1 ) in the following way: In the predictor step, an initial approximation (ỹ, t) of the new couple is generated in the direction of t k as ỹ := y k + ht k , t := t k , where h is a step size. Then, the corrector steps, which are iterative steps of Newton's type, are run with the initial approximation (ỹ, t). If they succeed in bringing the predicted point back to the currently approximated branch, the resulting couple is accepted, the step size h is adapted for the next predictor-corrector step, and the current step is done. Otherwise, h is reduced and the predictor and the corrector steps are repeated. If these predictor-corrector steps fail in computing a new couple, a special procedure is carried out for locating a new smooth branch (the so-called simple tangent H = H (1) H = H (2) H = H (3) H = H (4) yk-1 tk-1 yk tk ȳ Fig. 1 An example of the considered situation. switch). Afterwards, the predictor-corrector steps are restarted for continuation of the branch found. Let us emphasise that the whole method does not obey any analytical expressions of the regions of smoothness of H, and it is implementable generically. Similar continuation techniques constructed specially for static plane frictional contact problems can be found in [START_REF] Haslinger | Path-following the static contact problem with Coulomb friction[END_REF][START_REF] Haslinger | Qualitative analysis of solutions to discrete static contact problems with Coulomb friction[END_REF]. Nevertheless, all these methods were proposed for finding solely one continuing branch when an end point of the most recently traced smooth branch is encountered. Besides, a bifurcation criterion was proposed merely in [START_REF] Ligurský | Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction[END_REF].

The present paper deals with purely non-smooth numerical branching of solutions of (P) when H is P C 1 . We develop an approach that does not rely on any analytical expressions of the regions of smoothness of H in Section 2, and we test its performance on model examples of contact problems in Section 3.

Throughout the paper, we employ the following definition of a P C 1 -function [START_REF] Scholtes | Introduction to piecewise differentiable equations[END_REF]:

Definition 1.1 A function H : U → R N defined on an open set U ⊂ R M is P C 1 if it is continuous and for every ȳ ∈ U , there exist an open neighbourhood O ⊂ U of ȳ and a finite family of C 1 -functions H (i) : O → R N , i ∈ I(ȳ), such that ∀y ∈ O : H(y) ∈ {H (i) (y); i ∈ I(ȳ)}.
The functions H (i) are termed selections of H at ȳ.

Numerical Branching

Consider that an approximation of a smooth solution branch of (P) has been obtained by the restarted predictor-corrector method from [START_REF] Ligurský | A continuation problem for computing solutions of discretised evolution problems with application to plane quasi-static contact problems with friction[END_REF], and either the method has failed in locating a new smooth branch or it has located one but we want to find other ones because a bifurcation is expected by virtue of the criterion from [START_REF] Ligurský | Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction[END_REF], for example. Restricting ourselves to purely non-smooth branching (Assumptions (I)-(III) below), we shall describe a method designed for discovering potentially all continuing branches.

Let ȳ be an end point of the recovered branch. We shall assume the following:

(I) The whole solution set of (P) in a vicinity of ȳ is formed by one-sided smooth solution branches emanating from ȳ into mutually distinct regions of smoothness {y ∈ O; H(y) = H (i) (y)} for some i ∈ I(ȳ). (II) The gradients

∇H (i) (ȳ), i ∈ I(ȳ), have a full rank. (III) Ker ∇H (i) (ȳ) ∩ Ker ∇H (j) (ȳ) = {0}, ∀i, j ∈ I(ȳ), i = j.
Note that these assumptions guarantee particularly that tangent vectors at ȳ to any two different branches are linearly independent (see Fig. 1 for an illustration of branching under our consideration).

Since it seems to be hardly possible to encounter a point of non-smoothness in practical computations, we shall suppose in addition that H is smooth in any point considered in the procedures hereinafter. Then, a unit tangent vector t at such a point y lying in a vicinity of ȳ and belonging to a solution branch from the region {z ∈ U ; H(z) = H (i) (z)} is uniquely determined up to a direction by the conditions:

∇H (i) (y)t = 0, t w = 1,
where . w denotes a suitable weighted norm [11, pp. 86 and 87].

Our branching method consists of two steps: an approximation of the end point ȳ and a subsequent location of new branches from a neighbourhood of the approximate end point.

Approximation of the end point

The end point can be approximated by a bisection-like procedure that is based on the predictor-corrector method used in the process of numerical continuation of smooth branches. The proposed procedure relies on the fact that the individual branches emanating from the end point have numerically distinguishable tangents, which is justified by Assumptions (I)-(III). The algorithm can be sketched as follows:

Algorithm 2.1 Input: c diff ∈ (0, 1), h > h min > 0, y 0 , t 0 ∈ R N +1 such that H(y 0 ) ≈ 0, ∇H(y 0 )t 0 = 0, t 0 w = 1,
and the corrector steps for the prediction ỹ := y 0 + ht 0 , t := t 0 do not converge or they converge to a couple (y, t) with t t 0 /( t t 0 ) < c diff .

Step 1: Set h := h/2.

Step 2: Do the prediction ỹ := y 0 + ht 0 , t := t 0 .

Step 3: Run the correction with (ỹ, t).

Step 4: If the correction converges to a new couple (y, t) with t t 0 /( t t 0 ) ≥ c diff , then set y 0 := y, t 0 := t.

Step 5: If h < h min , break. Otherwise, go to Step 1.

Output:

The point y 0 with the corresponding unit tangent t 0 .

The couple (y 0 , t 0 ) in the input is supposed to be the last couple obtained by numerical continuation of the most recent smooth branch (the couple (y k , t k ) in Fig. 1), and h equals typically the minimal step size from the predictor-corrector method employed during the continuation. We denote the Euclidean norm by . and c diff is the minimal value of the cosine of the angle of two tangent vectors considered to correspond to the same branch. Hence, if the conditions in Step 4 are satisfied, the newly found point y is supposed to be from the same branch as the one stored in y 0 at that time, and it is accepted as a better approximation of ȳ. The minimal step size h min determines the precision required for the resulting approximation of ȳ.

Location of new branches

Suppose that y 0 calculated by Algorithm 2.1 is a good approximation of ȳ, it belongs to some smooth branch, and t 0 is the corresponding unit tangent pointing out from the region of smoothness with y 0 (see Fig. 2). Inspired by the methods for branch switching at a smooth simple bifurcation point [START_REF] Georg | Matrix-free numerical continuation and bifurcation[END_REF][START_REF] Hughes | A bisection-like algorithm for branch switching at a simple branch point[END_REF], we shall propose a heuristic technique composed of a sequence of predictor-corrector steps and yielding points on other smooth branches.

Our heuristic stems from the idea used in the simple tangent switch [START_REF] Ligurský | A continuation problem for computing solutions of discretised evolution problems with application to plane quasi-static contact problems with friction[END_REF]: Let t be a vector satisfying

∇H(ỹ) t = 0, t w = 1, (1)
where ỹ is a point close to ȳ and lies in the region of smoothness containing a new smooth branch (Fig. 3). Then, either t or -t should be a suitable prediction direction for a predictor-corrector step giving a point on this new branch.

Imagine for a while that you have a sequence of points ỹ that pass through all regions of smoothness intersecting at ȳ at your disposal. Under Assumptions (I)-(III), you could discover all solution branches by computing the sequence of vectors t from (1) and trying the corresponding predictor-corrector steps with t and -t successively. Hence, we would like to guess such a sequence of points ỹ from all the regions intersecting at ȳ.

For this purpose, consider the most probable branching scenarios described in [START_REF] Ligurský | Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction[END_REF]: By linearisation and projection, the corresponding regions of smoothness can be represented by cones in a two-dimensional vector space with vertices at 0. This representative space can be spanned by any two linearly independent vectors, say v 1 and v 2 . Thus, a sequence of linear combinations of v 1 and v 2 whose directions are deployed densely enough in the space comprises the directions pointing into the representatives of all regions of smoothness as indicated in Fig. 4. From here, one can conclude that if v 1 and v 2 represent some vectors t1 and t2 , then the corresponding sequence of linear combinations of t1 and t2 contains directions pointing into all regions and it can be used for calculating the sequence of points as desired.

In the general setting with (y 0 , t 0 ) from the output of Algorithm 2.1, let us choose a step length h sufficiently large in comparison with the precision of the approximation of ȳ by y 0 , in particular so large that y 0 + ht 0 leaves the region with y 0 . Next, let us take t1 := -t 0 and calculate a vector t2 satisfying

∇H(y 0 + ht 0 ) t2 = 0, t2 w = 1.
According to the imposed assumptions, t1 and t2 are linearly independent. One can choose a set B ⊂ R 2 so that {β 1 t1 + β 2 t2 } (β1,β2)∈B is a sequence of densely deployed directions. With regard to the argumentation for the most probable branching scenarios, one can expect that the sequence {y 0 + h(β 1 t1 + β 2 t2 )} (β1,β2)∈B might be the desired one containing points from all regions intersecting at ȳ. These considerations lead us to the following algorithm:

Algorithm 2.2 Input: c diff ∈ (0, 1), n dir ∈ N, h, h init > 0, y 0 , t 0 ∈ R N +1 such that H(y 0 ) ≈ 0, ∇H(y 0 )t 0 = 0, t 0 w = 1,
and y 0 + ht 0 leaves the region of smoothness with y 0 .

Step 1: Take t1 := -t 0 and compute t2 such that

∇H(y 0 + ht 0 ) t2 = 0, t2 w = 1.
Set n br := 1, y n br := y 0 , t n br := t1 .

Step 2: Set i := 0.

Step 3: Set

α := 2πi n dir , v := (sin α) t1 + (cos α) t2 , v := v v w ,
and compute t such that ∇H(y 0 + hv) t = 0, t w = 1.

Step 4:

If i = 0 or | t t 0 |/( t t 0 ) < c diff , set t 0 := t.
Otherwise, go to Step 7.

Step 5: Try one predictor-corrector step with the couple (y 0 , t) and the initial step size h init . If it converges to a new couple (y, t), then (i) if t (y -y 0 ) < 0, take t := -t;

(ii) if t t j /( t t j ) < c diff , j = 1, . . . , n br , then set n br := n br + 1, y n br := y, t n br := t.

Step 6: Repeat the procedure from Step 5 with t := -t.

Step 7: Set i := i + 1. If i < n dir , go to Step 3. Otherwise, break.

Output: Points y 1 , . . . , y n br with the corresponding unit tangents t 1 , . . . , t n br .

Here, c diff is the minimal absolute value of the cosine of the angle of two vectors from the kernels of gradients considered to correspond to the same selection similarly as in Algorithm 2.1. Further, n dir is the total number of linear combinations of t1 and t2 for seeking points from distinct regions. The step sizes h and h init should be appropriately large in comparison with the expected precision of the approximation of ȳ by y 0 so that the points determined with them fall into the regions in the respective directions.

The conditions in Step 4 serve for determining whether t has been computed from a gradient from the same region as previously. If it is so, it is superfluous, hence the corresponding predictor-corrector steps are skipped. In Step 5(ii), one verifies whether (y, t) does not belong to any of the branches already found before. This ensures that the points in the output are from different smooth branches. The corresponding tangents are oriented in the directions of branching from ȳ, which is guaranteed by Step 5(i). Afterwards, any of the branches located in this way can be traced by predictor-corrector steps.

Remark 2.3 To increase the probability of recovering all branches, one can restart Algorithm 2.2 several times from Step 2 with other choices of t1 and t2 . We propose the following possibilities:

(i) Take ti from the set {t 1 , . . . , t n br } at your current disposal. (ii) Take ti := t+ i with t+ i satisfying

∇H y 0 + h t- i + 0.1 t t w t+ i = 0, t+ i w = 1,
where the value of h is the same as in Algorithm 2.2, ti is equal to the vector ti employed in the previous run of Algorithm 2.2 and t is chosen randomly.

Numerical Tests

We have tested our technique on two examples of bifurcations in discretised elastostatic plane contact problems with Coulomb friction from [START_REF] Ligurský | Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction[END_REF]. We have chosen the following strategy for restarting Algorithm 2.2 based on Remark 2.3:

(1) Pass with t1 and t2 successively through the combinations from {t 1 , . . . , t n br } at your disposal.

(2) If all combinations which are available so far have already been employed, try to let t1 be the same and choose t2 according to Remark 2.3(ii). In our tests, we denote the total number of selections of t1 and t2 by n span , that is, the total number of the restarts of Algorithm 2.2 equals n span -1. Our calculations have been done with the finite-element library GetFEM++ [START_REF] Renard | Getfem++, An open source generic C++ library for finite element methods[END_REF].

Elastostatic plane contact problems with Coulomb friction

Let us consider static deformation of an elastic body whose reference configuration is the closure of a bounded domain Ω ⊂ R 2 . Let the boundary ∂Ω be Lipschitz-continuous and split into three disjoint relatively open subsets Γ D , Γ N and Γ c . The body is fixed along Γ D whereas an applied surface force of density h is prescribed on Γ N . A flat rigid foundation supports the body along Γ c , and the contact is modelled by unilateral conditions and the Coulomb friction law. We suppose that there is no initial gap between the body and the foundation. We consider that the surface force at a point x ∈ Γ N depends on a real parameter γ and it may depend on x either, that is, h = h(γ, x) in our model, γ being a loading parameter. We seek equilibrium states of the body for the values of γ from an interval I of our interest.

Discretisation of this problem is done by applying a conforming Lagrange finite-element method to a mixed variational formulation of the problem with Lagrange multipliers enforcing the Dirichlet and the contact boundary conditions. The contact conditions are approximated nodally.

Let {φ i } 1≤i≤2nΩ and {ξ i } 1≤i≤2n D be bases of the shape functions of the finite elements for the displacement and the Lagrange multiplier corresponding to the Dirichlet condition, respectively. We introduce the following vectors and matrix:

A(u) = (A i (u)) 1≤i≤2nΩ , A i (u) := Ω σ x, I + 2nΩ j=1 u j ∇φ j (x) : ∇φ i (x) dx, u = (u j ) 1≤j≤2nΩ , L(γ) = (L i (γ)) 1≤i≤2nΩ , L i (γ) := Γ N h(γ, x) • φ i (x) ds, B D = (B D,ij ) 1≤i≤2n D 1≤j≤2nΩ , B D,ij := Γ D ξ i (x) • φ j (x) ds.
Here, σ stands for the first Piola-Kirchhoff stress tensor, which is replaced by the Cauchy stress tensor in the case of smalldeformation elasticity (the second example). Further, let {a i } 1≤i≤nc be the set of the finite-element nodes corresponding to {φ i } 1≤i≤2nΩ and lying on Γ c \ Γ D , and ν and τ denote the unit inward normal and the unit tangent to the foundation, respectively. One can define matrices B ν , B τ ∈ R nc×2nΩ so that for any u h = 2nΩ j=1 u j φ j ,

(B ν u) i = u h ν (a i ), (B τ u) i = u h τ (a i ), i = 1, . . . , n c ,
where u h ν and u h τ are the components of u h in the directions of ν and τ :

u h ν (a i ) = u h (a i ) • ν, u h τ (a i ) = u h (a i ) • τ .
The discretised problem can be written in the following form of (P) [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers[END_REF]:

Find y := (γ, u, λ D , λ ν , λ τ ) ∈ I × R 2(nΩ+n D +nc) such that H(y) = 0,
where

H : I × R 2(nΩ+n D +nc) → R 2(nΩ+n D +nc) is defined by H(y) :=     A(u) -L(γ) -B D λ D -B ν λ ν -B τ λ τ B D u -1 r λ ν,j -((λ ν -rB ν u) j ) -, j = 1, . . . , n c -1 r λ τ,j -P [F ((λν -rBν u)j )-,-F ((λν -rBν u)j )-] ((λ τ -rB τ u) j ) , j = 1, . . . , n c     . (2)
Here, u ∈ R 2nΩ is the vector of nodal displacements, λ D ∈ R 2n D is the Lagrange multiplier corresponding to the Dirichlet condition, and λ ν , λ τ ∈ R nc are the normal and tangential Lagrange multipliers on the contact zone, respectively. The notation (.) -means the negative part (or equivalently, the projection onto (-∞, 0]), and P [a,b] (.) stands for the projection onto an interval [a, b], a ≤ b. The friction coefficient is represented by a non-negative constant F here and r > 0 is an augmentation parameter. Since the projections onto an interval are P C 1 -functions, H is also P C 1 provided that A and L are continuously differentiable.

Let us mention that if h denotes the mesh size for the displacement, the parameter r can be viewed as a coefficient of the proportionality between the Lagrange multipliers divided by h, which gives approximations of the contact stresses, and the displacements divided by h on the contact zone (compare to the physical interpretation of the modulus of elasticity). Nevertheless, the (exact) solutions of the discrete problem do not depend on the value of r because the formulations are mathematically equivalent for all r > 0. Its choice affects only numerical solving of the problem.

Problem with a rectangular body

Consider deformation of a rectangular block that is 40 mm wide and 80 mm high, fixed along its upper side and loaded along its lateral sides (Fig. 5). A plane-strain approximation of the nonlinear Ciarlet-Geymonat constitutive law [START_REF] Ciarlet | I: Three-dimensional elasticity[END_REF]Chapter 4] is used for the material of the block. Namely, the first Piola-Kirchhoff stress tensor σ is given by

σ(x, F ) = ( σ( F )) 1≤i,j≤2 , F = F 0 0 1 , F ∈ R 2×2 , σ( F ) = 2b tr( F F ) I + 2(a -b F F ) F + 2c det( F F ) -d F -, F ∈ R 3×3 , where λ = 4000 N/mm 2 , µ = 120 N/mm 2 , a = 30 N/mm 2 and b = µ 2 -a, c = λ 4 - µ 2 + a, d = λ 2 + µ.
We have prescribed h(γ, x) = γ(-2, 0.12(x 1 -20)) (in N/mm 2 ) on both lateral sides of the block, r = 10 and F = 1.

Discretisation is done by approximating the displacement and the Lagrange multiplier for the Dirichlet condition with continuous piecewise bilinear functions on a uniform mesh with 800 squares and 21 contact nodes. There are six solution branches intersecting at γ = 0, which are illustrated in Fig. 5. Branches 1, 2 and 3 correspond to forcing the block to the right with no contact, contact-stick and contact-slip to the right of the lower right vertex of the block. Branches 4, 5 and 6, which are symmetric with respect to the axis x 1 = 20 of the block, correspond to forcing the block to the left with no contact, contact-stick and contact-slip to the left of the lower left vertex.

In the testing of our method, we have set c diff = 0.99999, h = 5e-7, h min = h/2 14 . = 3.1e-11 and h init = 5e-4 in Algorithms 2.1 and 2.2. To start with, we have chosen various values of n dir , and we have looked into sufficient values of n span for finding all six branches. Our computations are summarised in Table 1, where the entry Input means the branch from which y 0 is chosen in the input of Algorithm 2.1 (only the listings for Branches 1, 2 and 3 are written down because the behaviour for Branches 4, 5 and 6 is the same due to the symmetry of the solutions). The cost of each computation is characterised by n pc , which stands for the total number of the predictor-corrector steps. These constitute the most expensive part of our technique. Note that random generations according to Remark 2.3(ii) are needed for some input branches and some values of n dir . In such cases, we have run our branching method three times, and the behaviour of the performed computational process has varied from one run to another. We denote the total number of the restarts of Algorithm 2.2 following a random generation by n rand here and in what follows. Afterwards, we have prescribed n dir as well as n span , and we have compared the obtained values of n pc to get an idea about the cheapest but quiet reliable strategy. The overview of the performance of our method is presented in Table 2. All six branches have been recovered in all cases. According to the results, we can propose to set n dir to about 20 and n span to about 15 in this particular example. In addition, we can suggest to increase rather n dir than n span to increase the probability of finding all branches at the lowest price in general.

To explore the bifurcation phenomenon in the contact problem thoroughly, we have performed the method of numerical branching for discretisations on different uniform meshes, namely with 200, 800, 3200 and 12800 bilinear squares (M200, 8 720 3

Table 3 The relative error of the displacements for γ = -1 for the rectangular body. The displacements u h 1 , u h 2 , u h 3 and u h 4 have been computed for M200, M800, M3200 and M12800, respectively.

Branch 1 2 3 u h1 -u h4 1,Ω / u h4 1,Ω 6.1e-2 1.7e-1 4.3e-2 u h2 -u h4 1,Ω / u h4 1,Ω 2.3e-2 1.1e-1 1.1e-2 u h3 -u h4 1,Ω / u h4 1,Ω
6.9e-3 5.1e-2 5.4e-3 M800, M3200 and M12800). Having taken n dir = 30 and n span = 15, we have found six branches distinctly separated for all meshes. The branches seem to be stable and to converge with the meshes, see the bifurcation diagrams in Fig. 6 obtained by the method of piecewise-smooth numerical continuation from [START_REF] Ligurský | A continuation problem for computing solutions of discretised evolution problems with application to plane quasi-static contact problems with friction[END_REF], and the behaviour of the relative error for Branches 1, 2 and 3 in Table 3 (which is the same as the one for Branches 4, 5 and 6).

Fig. 6 The normal and tangential displacements uν and uτ at the most left contact node of the rectangular body.

Problem with a triangular body

Next, let us consider deformation of an isosceles right-angled triangle with the legs 1 m long, lying on a leg, loaded along the other leg and fixed along the hypotenuse (Fig. 7). Let us restrict ourselves to the small-deformation framework with Hooke's law with the Lamé constants λ = 100 GN/m 2 and µ = 82 GN/m 2 , and let r = 10 and F = 1.7. Both the displacement and the Lagrange multiplier for the Dirichlet condition are approximated with continuous piecewise linear functions over triangular meshes. As described in [START_REF] Ligurský | Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction[END_REF], there are four solution branches intersecting at γ = 0 for h = h(γ) = γ(-26, -7.5) GN/m 2 and the discretisation with a uniform mesh with 4096 triangles and 64 contact nodes, see Fig. 7. They correspond to a partial contact and slip of the triangle to the right (Branch 2), and to no contact, contact-stick and contact-slip to the left of the lower left vertex of the triangle with pulling the whole triangle to the left (Branches 1, 3 and 4).

To explore the bifurcation in this contact problem, we have taken uniform meshes with 4096, 16384, 65536 and 262144 linear triangles (M4096, M16384, M65536 and M262144), and we have used our method with c diff = 1 -1e-10, h = 5e-8, h min = h/2 14 . = 3.1e-12, h init = 5e-4, n dir = 30 and n span = 10. The bifurcation behaviour has shown to be more complex here. Branches 1 and 4 approach one another for finer meshes, and they disappear both for the finest mesh for the force h(γ) = γ(-26, -7.5), as illustrated in Fig. 8.

To investigate this behaviour more precisely, we have run an auxiliary numerical continuation of the same parametrised problem except the force h(γ) = (-20 -10γ, -7.5). As one can observe in Fig. 9, there are values of γ for which the contact problem has multiple solutions for this parametrisation, but these values differ a bit for each mesh. Namely, the problem has three solutions for γ approximately from (0.56541, 0.66665), (0.5883, 0.70731), (0.59996, 0.73289) and (0.60585, 0.7499) for M4096, M16384, M65536 and M262144, respectively. Moreover, two solutions approach one another when γ approaches a bound of the interval corresponding to a fixed mesh.

Besides, note that A in ( 2) is linear in this case. If h is considered linear with respect to γ, then L is linear, as well, and H given by ( 2) is piecewise linear and especially positively homogeneous. This yields that the solutions obtained by the auxiliary continuation characterise the solution behaviour of the corresponding unparametrised contact problem with forces h = (h 1 , h 2 ) for any h 2 < 0.

In particular, we have discovered three solutions of the unparametrised contact problem with forces h = (h 1 , -7.5) for h 1 approximately from (25.6541, 26.6665), (25.883, 27.0731), (25.9996, 27.3289) and (26.0585, 27.499) for the meshes M4096, M16384, M65536 and M262144, respectively. This and the positive homogeneity give us cones in the h 1 -h 2 plane such that the unparametrised problem has three solutions for the forces from their interiors for each mesh. Furthermore, if a force approaches a half-line from the boundary of the cone for a given mesh, two of the three corresponding solutions approach one another.

Since the line h(γ) = γ(-26, -7.5) is closer and closer to a half-line from the boundaries of the cones for finer and finer meshes and out of the cone for the finest mesh, this explains the behaviour of the bifurcation in the present example. Nevertheless, regarding the branching in the whole h 1 -h 2 plane, one can find it stable and convergent, again. This is illustrated by the bifurcation diagrams in Fig. 10 and Table 4 computed for a bit modified parametrised force h(γ) = γ(-26.4, -7.5), for which there are four branches distinctly separated for all meshes. Fig. 8 The normal and tangential displacements uν and uτ at the most left contact node of the triangular body for h(γ) = γ(-26, -7.5) (in GN/m 2 ).

Fig. 9 The normal and tangential displacements uν and uτ at the most left contact node of the triangular body for h(γ) = (-20 -10γ, -7.5) (in GN/m 2 ).

Fig. 10 The normal and tangential displacements uν and uτ at the most left contact node of the triangular body for h(γ) = γ(-26.4, -7.5) (in GN/m 2 ).

Table 4 The relative error of the displacements for γ = -1 (Branch 2) and γ = 1 (Branches 1, 3 and 4) for the triangular body. The displacements u h 1 , u h 2 , u h 3 , u h 4 have been computed for M4096, M16384, M65536 and M262144, respectively.

Conclusion

We have developed a method of numerical branching for piecewise-smooth problems, and we have suggested how to choose its parameter setting to discover reliably all solution branches at the lowest price. We have used our method for investigating the behaviour of bifurcations in discretised contact problems with Coulomb friction for different meshes. As far as we know, this is the first attempt to devise a numerical method for piecewise-smooth branching that does not rely on any analytical expressions of the regions of smoothness of the function involved.

Fig. 2

 2 Fig. 2 Output of Algorithm 2.1.

Fig. 3 Fig. 4

 34 Fig. 3 Prediction direction for finding a new branch.

Fig. 5

 5 Fig. 5 Deformed bodies corresponding to the solutions with |γ| = 1: (a), (b), (c) γ = -1; (d), (e), (f) γ = 1.

Fig. 7

 7 Fig. 7 Deformed bodies corresponding to the solutions with |γ| = 1: (b) γ = -1; (a), (c), (d) γ = 1.

Table 1

 1 Rectangular body, n dir given a priori. Input n dir n span n pc n rand Input n dir n span n pc n rand

		20	8 168	1		20	4 136	1
		20	10 214	3		20	7 220	3
		20	13 266	6		20	8 240	4
		40	7 292	2		40	4 260	1
	Branch 1	40	8 330	3		40	5 294	2
		40	14 478	7		40	5 316	2
		50	3 162	0		50	4 274	1
		60	5 232	0	Branch 2	50	5 342	2
		100	3 206	0		50	7 490	4
		20	5 128	0		60	5 374	2
		40	3 134	0		60	8 620	5
	Branch 3	50	1	72	0		60	9 654	6
		60	4 202	0		100	4 372	1
		100	1 102	0		100	4 406	1
							100	7 662	4

Table 2

 2 Rectangular body, n dir and nspan given a priori.Inputn dir n span n pc n rand

					Input	n dir n span n pc n rand
		20	15 310	3		20	15 350	2
		20	15 316	5		20	15 360	4
		20	15 322	1		20	15 416	8
		20	16 326	3		20	16 374	2
	Branch 1	20	16 344	5		20	16 386	4
		20	16 346	1		20	16 434	8
		50	10 400	0		50	10 518	1
		60	10 416	0	Branch 2	50	10 532	1
		100	8 474	0		50	10 562	2
		20	15 318	0		60	10 618	2
		20	15 318	0		60	10 642	3
		20	16 334	1		60	10 734	6
	Branch 3	20	16 334	1		100	8 642	1
		50	10 396	0		100	8 684	2
		60	10 414	0		100
		100	8 456	0

  excellent research teams and inter-sectoral mobility at Palacký University in Olomouc II" (CZ.1.07/2.3.00/30.0041), and by The Ministry of Education, Youth and Sports of the Czech Republic from the National Programme of Sustainability (NPU II) project "IT4Innovations excellence in science -LQ1602".

Acknowledgements The basis of this study was worked out at a postdoctoral stay of T. Ligurský at INSA-Lyon funded by Manufacture Franc ¸aise des Pneumatiques Michelin. During further developments, T. Ligurský was supported from the project "Support for building


C Brun

C Crauste-Thibierge

D L'h Ôte

F Ladieu

G Biroli

J.-P Bouchaud

Evidence of growing spatial correlations at the glass transition from nonlinear susceptibility measurements

Keywords: glasses, nonlinear, susceptibility, dielectric

Evidence of growing spatial correlations at the glass transition from nonlinear susceptibility measurements

We have measured the ac nonlinear dielectric response χ 3 (ω, T ) of an archetypical glassformer (glycerol) just above its glass transition temperature T g ∼ 190 K. Our measurements were performed at temperatures between T g + 4 K and T g + 35 K, which corresponds to 5 decades of variation of the relaxation time τ (T ) of the liquid. We find that χ 3 (ω, T ) is peaked for ωτ ∼ 0.2 and that the height of the peak grows as one approaches T g . A powerlaw in frequency beyond the peak of χ 3 is observed and χ 3 (ω, T ) displays scaling as a function of ωτ . These features correspond to the main predictions of some recent theoretical works relating the nonlinear susceptibility of supercooled liquids to the average number N corr of dynamically correlated molecules. We thus interpret the experimental increase of the peak of nonlinear susceptibility when T decreases, as an increase of N corr when T → T g . This strongly reinforces the collective picture of the dynamics of supercooled liquids close to the glass transition. Finally, we discuss some heating effects and show that they do not affect significantly our conclusions.

1 Introduction Upon fast enough cooling, most liquids avoid crystallization and enter in a supercooled state. In this supercooled state, the "viscous slowing down" phenomenon arises [1], i.e. the viscosity η of supercooled liquids increases extremely fast when the temperature T is decreased, this increase being even faster than Arrhenius in the so-called fragile glasses. Below T g the viscosity is so high that the system is, in practice, a solid, the glass. This glass transition is a longstanding issue of condensed matter physics, since no direct structural signature has ever been detected around T g , contrary to what happens at the crystallization transition [1]. Some progresses about the glass transition were made in the two past decades, where some experimental breakthroughs [2,3] as well as numer-ical works [6] established the heterogeneous nature of the dynamics of supercooled liquids. These progresses led to the important concept of "dynamical heterogeneities" [4,5], according to which the relaxation events come from collective motions of groups of typically N corr molecules: some of these groups are faster than others, but none of them are permanent entities, they all evolve in time. It has been proposed that the huge increase of η in supercooled liquids comes from the fact that this number N corr increases as T decreases towards T g . The variation of N corr with temperature thus became an important issue in the field, triggering new theoretical ideas [7][8][9][10]. One of these new ideas came from an analogy [10] with the well known spin glass physics where the nonlinear susceptibility di-verges at T c , reflecting the long range amorphous order which sets in at the spin glass transition temperature T c . Let us underline that no divergence occurs in the linear susceptibility of spin glasses due to the fact that two point correlation functions are "blind" to this amorphous long range order. This peculiar order is captured only by higher order correlations functions, and thus by their associated nonlinear responses (through fluctuation-dissipation relations). Based upon a generalized fluctuation-dissipation relation, it was proposed in Ref. [10], that a similar phenomenon occurs in supercooled liquids, with the key difference that the peak of the nonlinear response χ 3 was predicted to occur at finite frequencies of order 1/τ (T ) where τ (T ) is the relaxation time of the supercooled liquid. For the first time, a susceptibility was directly related to the quantity N corr of interest, allowing to scrutinize its temperature dependence.

2 Experimental setup Motivated by the new prediction of Ref. [10], we have developed an experiment allowing to measure the nonlinear susceptibility of supercooled liquids by the detection of the third harmonics of the dielectric polarisation P (t) [11]. The supercooled liquid fills the space between two thick metallic electrodes separated by a distance e, and is submitted to a field E cos(ωt). As a result, a dielectric response P 1 arises at ω. P 1 is dominated by the linear response P lin , corresponding to the linear susceptibility χ lin (ω) + iχ lin (ω) [11], and hereafter the frequency where χ lin is maximum defines the relaxation frequency f α , the latter being of order 1/τ . Due to the nonlinear behavior of the system, P (t) contains higher harmonics at frequencies nω, where n is an odd integer for symmetry reasons (the even values of n would be allowed if the field would contain a static part). The response P 3 at 3ω is thus the first signal which is directly proportional to the nonlinear properties of the liquid. In pratice, P 3 is dominated by the third harmonics susceptibility χ 3 of the liquid (see [11,12] for a precise definition of χ 3 ). Due to the weakness of P 3 /P 1 , we have developed a special bridge technique, containing two samples of different thicknesses e thin < 20μm and e thick 2e thin . This allows to get rid of all the spurious 3ω signals coming both from the source and the voltage amplifier, and yields a sensitivity for P 3 /P 1 in the range of 10 -7 . Our experiment is described in full details in Refs. [11,12]. The results reported here are for glycerol, an archetypical glass former with T g 190 K.

3 Results and discussion This high sensitivity setup was used to test the predictions of Ref. [10] which, as evoked in the introduction, gives a relation between χ 3 and the T dependent average number of dynamically molecules N corr (T ) through the following scaling form :

χ 3 (ω, T ) ≈ 0 (Δχ 1 ) 2 a 3 k B T N corr (T ) H (ωτ ) , (1)
1 (Color online) For glycerol, the modulus of the dimensionless nonlinear susceptibility X3(ω, T ) = χ3(ω, T ) × kBT /[(Δχ1) 2 a 3 0] is shown at two temperatures in the supercooled regime. For comparison, the corresponding quantity for independent polar molecules, labelled |X3| trivial , is plotted as a solid line (which does not depend on T ). Both the humped shape of χ3 versus frequency, as well as its T dependence, are distinctive features of glassy correlations. Inset: phase of X3 (same symbols, phase of X 3,trivial not shown).

where

Δχ 1 = χ 1 (ω = 0) -χ 1 (ω → ∞)
is the part of the static linear susceptibility corresponding to the slow relaxation process we consider, a 3 the volume occupied by one molecule, and H a certain complex scaling function that goes to zero both for small and large arguments. This humped shape of H comes out due to glassy correlations [10]. To emphasize this point, let us consider the case of independent polar molecules undergoing (non inertial) rotational Brownian rotation. Their nonlinear dielectric susceptibility is given (see [12,13]) by the first prefactor of the right hand side of Eq. ( 1) times a function of ωτ which has two main features: i) it reaches its maximum value at ω = 0, ii) it does not depend on T , once plotted as a function of ωτ . This is why, in the following, we shall divide χ 3 by the prefactor of the right hand side of Eq. (1), i.e. we shall focus onto the dimensionless nonlinear susceptibility 1) is relevant for supercooled liquids, one has X 3 (ω, T ) ≈ N corr (T ) H (ωτ ). On the contrary, for non-interacting molecules, |X 3 | is a steadily decreasing function of ωτ and does not depend on T . To illustrate this point, we have used the analytical results of Déjardin et al. [13], and plotted, on Fig. 1, the values of |X 3 | trivial for independent molecules as a continuous solid line. Here the label "trivial" emphasizes the fact that the solid line on Fig. 1 has nothing to do with glassy correlations.

X 3 (ω, T ) = χ 3 (ω, T ) × k B T /[(Δχ 1 ) 2 a 3 0 ]. If Eq. (
In Fig. 1, two sets of data for glycerol are given, corresponding to T = 210.2 K and to T = 217.8 K, to illustrate the main results of Ref. [14] where much more data are reported. First, as predicted in Ref. [10], one observes on Fig. 1 that |X 3 | is "humped". Following Eq. ( 1), this humped shape comes from the transient nature of the glassy correlations, the lifetime of which is of the order of ∼ α : for f /f α 1, the glassy correlations are not yet established, while for f /f α 1 the liquid "flows", which destroys the glassy correlations, leaving only the "trivial" correlations above mentioned. This is why it was anticipated in [10] that the maximum of |X 3 | should arise for a frequency f of the order of f α ; and Fig. 1 shows indeed that f 0.21f α . Furthermore, |X 3 | decays as a power law when f ≥ f α as predicted in Ref. [15], this point is illustrated in Fig. 2 of Ref. [14]. Last, defining, for each considered T , X 3,max (T ) as the maximum of |X 3 (ω, T )| over frequency, one observes that |X 3 (ω, T )/X 3,max (T )| falls onto a unique master curve, depending only on f /f α , in agreement with the prediction of Eq. ( 1). This scaling property can be guessed from the main part of Fig. 1, see also Fig. 2 of Ref. [14], as well as from the inset where it is shown that the phase of X 3 falls onto a T -independent master curve. Departures from this scaling occur at very low frequencies and/or at the highest temperature of 225 K reported in Ref. [14]. This is not surprising since in these two limits the role of the trivial nonlinearities evoked above is expected to contribute significantly to the measured X 3 .

Having checked that the main predictions of Eq. ( 1) are relevant for the nonlinear susceptibility of glycerol, we now move to the T -dependence of X 3,max (T ). As shown in Fig. 2, X 3,max increases when T decreases. By using Eq. ( 1), we interpret this behavior as giving directly the temperature dependence of N corr (T ) when T decreases towards T g . This is one of the main results of our study, and a strong experimental evidence of a growth of the dynamical correlation length close to T g . To further support this interpretation, we plot on Fig. 2 the T dependence of T χ T defined by T χ T = max ω (T ∂(χ lin (ω)/Δχ 1 )/∂T ). As explained in, e.g. Ref. [8], T χ T should be proportional to N corr (T ) up to an unknown prefactor proportional to χ 0 , where χ 0 is related to the correlation between the lo-Figure 3 (Color online) For glycerol |X3| is compared at T = 204.5 K to the heating contribution arising from the modulation of temperature at the frequency 2ω coming from the dissipation of electrical power in the supercooled liquid. The difference between the overestimated and damped heating contribution is explained in the text and in Ref. [12]. The heating contribution is in any case negligible with respect to the values of |X3| measured in glycerol. To ease the comparison with heating at very large values of f/fα, an extrapolation of the data of glycerol is plotted as a dotted line, following the power law behaviour reported in Ref.

[14] at 210.2 K. cal fluctuation of enthalpy and the orientation of molecules [16]. If the amplitude of this "enthalpy-orientations" correlation is T -independent, then T χ T captures the variation of the range of this correlation with T , and χ 0 is Tindependent. In that case, T χ T is a fair estimator of N corr , up to an unknown constant prefactor [16]. Assuming that χ 0 does not depend on T for Glycerol, one can compare the T -dependence of T χ T and of X 3,max . In Fig. 2, one sees that the T -dependence of both quantities are close (the unknown prefactor was used to make T χ T coincide with X 3,max arbitrarily at 202 K). This shows that the simpler method, introduced in [7] and used extensively in [9], where N corr (T ) is drawn from T χ T is valid at least for glycerol.

Expanding on the arguments given in Refs. [12,14], we now comment on the possible contribution of the heating of the liquid to our third harmonics measurements. This heating comes simply from the fact that the strong applied electric field leads to some dissipated electrical power p with a d.c. component as well as a component p 2 (t) oscillating at 2ω. Before being absorbed by the thermal reservoir setting the base temperature T of the experiment (hereafter represented by the metallic electrodes), this electrical power has to travel across the sample of thickness e. As a result a small temperature increase arises, containing a 2ω component δT 2 ∼ p 2 . By using, as an order of magnitude estimate δP (t) (∂P lin (t)/∂T )δT (t), which involves a product of two terms oscillating at ω and 2ω, one sees that a spurious contribution to the measured P 3 arises from δT 2 (t), i.e. from a part of the heating phenomenon.

In order to estimate this spurious effect, we have to compute δT 2 (t) for a supercooled liquid of total specific heat c, and thermal conductance κ. As shown in Ref. [17], the latter quantity does not depend on frequency, while the former depends on f /f α due to the fact that the slow degrees of freedom relevant for the glass transition cannot equilibrate for frequencies larger than f α . The heat propagation equation for such a liquid surrounded by two thermal reservoirs (i.e. the electrodes) separated by the distance e, is solved by decomposing the heating into spatial modes of wave vector K = mπ/e with m an odd integer. After spatially averaging the heating over the sample volume, one finds, keeping only the dominant m = 1 spatial mode [12]:

δT 2 (t) = 0 χ lin ωE 2 e 2 24κ th cos(2ωt -φ 2 ) 1 + (2ωτ th ) 2 (2)
where φ = -π + 2 arctan(χ 1 /(χ 1χ 1 (∞))), τ th = ce 2 /(κ th π 2 ), φ 2 = φ + arctan(2ωτ th ) and 0 is the dielectric permittivity of vacuum.

The simplest way to estimate the heating contribution P 3,h is to set P 3,h = (∂P 1 (t)/∂T )δT 2 (t), as evoked above. This clearly overestimates P 3,h in the case f ≥ f α because δT 2 (t) oscillates at a frequency larger than the typical relaxation frequency f α of the dipoles. A first estimate of the damping of the effect of δT 2 on the polarisation is given at the end of Ref. [12]. We shall not reproduce the equations here. We shall just illustrate the heating contribution by Fig. 3, where the overestimated and damped heating contributions to |X 3 | are plotted at T = 204.5 K and compared to the measured values of |X 3 | at the same temperature. The heating contribution is clearly negligible at any frequency. A detailed analysis (to be published elsewhere) shows that the heating contribution |X 3 (f /f α )| h globally increases with T -contrarily to |X 3 (f /f α )|-. The fact that heating is negligible at 204 K at any value of f /f α thus guarantees that this is also true for any lower temperature. Last in the vicinity of f 0.21f α where |X 3 | is maximum, the heating contribution is so small that it does not affect the data reported in Fig. 2, even at the highest temperature of 225.3 K.

Let us emphasize that the heating that we have just discussed above must not be confused with that of the "inhomogeneous heating model" [18][19][20]. This "inhomogeneous heating model" intends to give a phenomenological description of the intrinsic nonlinear effects in supercooled liquids, and was shown to account for the nonlinear measurements at ω in, e.g., references [18][19][20]. In this framework, each dynamical heterogeneity has its own fictitious temperature on top of the temperature of the phonon bath. Some more work is needed to give a thorough comparison of this model at 3ω and our data.

Conclusion

To summarize, we have discussed the behavior of the nonlinear susceptibility χ 3 (ω, T ) of glycerol for a temperature interval from T g + 4K to T g + 35K.

The predictions of Eq. ( 1) have been shown to be relevant for the experimental data, and this is why we interpret the T dependence of the height of the peak of χ 3 as reflecting directly that of N corr (T ). We think that any model aiming at describing the nonlinear susceptibility near the glass transition should be able to account for both the modulus and the phase of X 3 reported here. Indeed, these two quantities should contain important information about the dynamical correlations in the considered supercooled liquid. Furthermore accounting for the modulus and the phase of X 3 will allow to obtain the absolute value of N corr (T ), which is not possible at present due to the unknown expression of the function H. This is why, one can hope that a detailed understanding of the kind of data reported here will yield substantial progress in the modelisation of the glass transition, and perhaps to unveil its possible link with critical phenomena.

Figure 2 (

 2 Figure2(Color online) For glycerol, the maximum value X3,max of the modulus of X3(ω) over frequency is plotted as a function of T . Its increase when T decreases is interpreted as reflecting directly that of Ncorr(T ) (see text). The solid line is the T dependence of T χT with an arbitrary normalisation at 202 K (see text).

www.pss-c.com© 2011 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

Acknowledgements We thank R. Tourbot for his outstanding technical help. We acknowledge interesting discussions with C. Alba-Simionesco, A. Lefèvre, R. Richert, M. Tarzia, and support by ANR grant DynHet.


domain: sdv.bv.bot phys nlin

text:
9!

While many molecular players involved in growth control have been identified in the past 10! decades, it is often unknown how they mechanistically act to induce specific shape changes 11! during development. Plant morphogenesis results from the turgor-induced yielding of the 12! extracellular, load-bearing cell wall. Its mechano-chemical equilibrium appears as a fundamental 13! link between molecular growth regulation and the effective shape evolution of the tissue. Here, 14! we focus on force-driven polymerization of the cell wall as a central process in growth control.

15!

We propose that mechanical forces facilitate the insertion of wall components, in particular 16! pectins, a process that can be modulated through genetic regulation. We formalize this idea in a 17! mathematical model, which we subsequently challenge with published experimental results.

Introduction

19!

The link between gene regulation and growth, i.e. the irreversible spatial expansion of a cell or a 20! tissue, is a major issue in plant and animal developmental biology. However, although many 21!

5!

components. Such a coupling between mechanics and biochemistry, termed force-driven 63! polymerization, has also been described for the assembly of actin bundles in animal systems [START_REF] Kruse | Aster, Vortices and Rotating Spirals in Active Gels of Polar 327! Filaments[END_REF].

64!

Starting from known descriptions of cell wall structure (see [START_REF] Wolf | Growth Control and Cell Wall Signaling in Plants[END_REF] and references therein) we 65! expose how force-driven polymerization could be central in scaling up molecular processes into 66! viscoplastic behavior and ultimately growth of multicellular structures.

67!

The cell wall is a fiber-reinforced hydrogel

68!

The cell wall is composed of a network of rigid cellulose microfibrils embedded in a matrix, 69! mostly composed of water, polysaccharides such as pectins and hemicelluloses, proteins and ions 70! [START_REF] Caffall | The structure, function, and biosynthesis of plant cell 331! wall pectic polysaccharides[END_REF]. This matrix corresponds to a biphasic mixture between a porous solid phase and a liquid 71! one [START_REF] Iijima | Swelling behaviour of calcium pectin hydrogels by 333! thermomechanical analysis in water[END_REF]. Therefore, the cell wall as a whole can be considered as a fiber-reinforced hydrogel.

72!

From a functional point of view, the wall performs two contrasting functions, as it has to resist 73! turgor-induced constraints whereas simultaneously it should be able to expand. Because of their 74! mechanical characteristics [START_REF] Vincent | From cellulose to cell[END_REF][START_REF] Zhang | Visualization of the nanoscale pattern of recently-deposited cellulose 336! microfibrils and matrix materials in never-dried primary walls of the onion epidermis[END_REF], the cellulose microfibrils tethered to each other by shorter 75! hemicellulose chains are major actors in both of these functions (for reviews see: [START_REF] Wolf | Growth Control and Cell Wall Signaling in Plants[END_REF][START_REF] Cosgrove | Re-constructing our models of cellulose and primary cell wall 339! assembly[END_REF] [START_REF] Cosgrove | How do plant cell walls extend?[END_REF].

76!

The available experimental evidence suggests that cellulose microfibrils deposition is involved in 77! growth inhibition rather than in growth promotion per se. Indeed, after drug-induced inhibition of 78! cellulose deposition, cells can continue to grow during several tens of hours [START_REF] Heisler | Alignment between PIN1 Polarity and Microtubule Orientation 342! in the Shoot Apical Meristem Reveals a Tight Coupling between Morphogenesis and Auxin 343! Transport[END_REF][START_REF] Pelletier | A role for pectin de-methylesterification in a developmentally 345! regulated growth acceleration in dark-grown Arabidopsis hypocotyls[END_REF]. Similarly, 79! growth also occurs in certain mutant backgrounds where cellulose synthesis is impaired [START_REF] Fagard | PROCUSTE1 Encodes a Cellulose Synthase Required for Normal 348! Cell Elongation Specifically in Roots and Dark-Grown Hypocotyls of Arabidopsis[END_REF][START_REF] Hématy | A receptor-like kinase mediates the response of Arabidopsis cells 351! to the inhibition of cellulose synthesis[END_REF]. 80! These experimental clues suggest that growth of the cell wall relies on the expansion of the 81! embedding matrix. In other words a scenario has emerged, where cell wall expansion is restricted 82! by the synthesis of cellulose microfibrils and promoted by the synthesis of matrix material. In the 83! following paragraphs we will discuss an important component of this matrix, i.e. pectin. As we 84! ! 6! will see, there is increasing evidence that pectin could play an important role in growth control.

85!

Pectin is an important structural component of the wall matrix

86!

As indicated above, hemicelluloses such as xyloglucans have been classically considered as 87! essential structural components, central to concepts of wall loosening and induction of growth. 88! However, the recently produced Arabidopsis thaliana double mutant (xxt1/xxt2), showing a loss 89! of xyloglucan synthesis, does not show major growth deficiencies [START_REF] Park | A revised architecture of primary cell walls based on 353! ! 18! biomechanical changes induced by substrate-specific endoglucanases[END_REF][START_REF] Cavalier | Disrupting two Arabidopsis thaliana xylosyltransferase genes 356! results in plants deficient in xyloglucan, a major primary cell wall component[END_REF], suggesting that the 90! irreversible expansion of the cell wall does not rely, at least totally, on hemicellulose dynamics 91! only. This has revived the interest in the role of pectins in growth control. Pectins form up to 35 92! % of the primary cell wall of dicots [START_REF] Caffall | The structure, function, and biosynthesis of plant cell 331! wall pectic polysaccharides[END_REF] i.e. the wall surrounding young and fast growing cells.

93!

Their dynamics have been associated with growth control in several cell types, in particular tip-94! growing cells such as pollen tubes and root hairs, for review see [START_REF] Wolf | Growth Control and Cell Wall Signaling in Plants[END_REF] and references therein, see 95! also e.g. [START_REF] Rojas | Chemically Mediated Mechanical Expansion of the Pollen Tube 359! Cell Wall[END_REF], [START_REF] Parre | Pectin and the role of the physical properties of the cell 361! wall in pollen tube growth of Solanum chacoense[END_REF]. Moreover, AFM-based assays on Arabidopsis thaliana, demonstrated that 96! pectin demethyl-esterification was required and tightly regulated during cell wall softening prior 97! to aerial organ outgrowth [START_REF] Braybrook | Mechano-Chemical Aspects of Organ Formation 363! in Arabidopsis thaliana: The Relationship between Auxin and Pectin[END_REF], [START_REF] Peaucelle | The Control of Growth Symmetry Breaking in the Arabidopsis 365! Hypocotyl[END_REF]. In addition, recent evidence indicates the existence of close 98! interactions between cellulose and pectins which could function as tethers between the 99! microfibrils [START_REF] Cosgrove | Re-constructing our models of cellulose and primary cell wall 339! assembly[END_REF][START_REF] Wang | Cellulose-Pectin Spatial Contacts Are Inherent to Never-Dried 367! Arabidopsis Primary Cell Walls: Evidence from Solid-State Nuclear Magnetic Resonance[END_REF]. Dwarf phenotypes in mutants where pectin synthesis or delivery is impaired 100! further illustrate the importance of pectin production for plant growth [START_REF] Bouton | QUASIMODO1 encodes a putative membrane-bound 370! glycosyltransferase required for normal pectin synthesis and cell adhesion in Arabidopsis. 371![END_REF][START_REF] Gendre | Trans-Golgi network localized ECHIDNA/Ypt interacting protein 373! complex is required for the secretion of cell wall polysaccharides in Arabidopsis[END_REF]. Noteworthy, 101! pectins have also been identified as important regulators of growth in more primitive plants such 102! as the green algae [START_REF] Proseus | Calcium pectate chemistry controls growth rate of 376! Chara corallina[END_REF][START_REF] Proseus | Pectate chemistry links cell expansion to wall 378! deposition in Chara corallina[END_REF]. roughly account for 65% of all of them [START_REF] Mohnen | Cell wall polysaccharide synthesis[END_REF]. A fundamental feature of homogalacturonans is 108! their ability to form multiple intermolecular chelation bonds in presence of calcium cations.

109!

These cations-dependent bonds are able to glue together pectin molecules, forming structures 110! called "egg-boxes" (Fig. 1). In a cation-rich environment, these adhesive interactions lead to the 111! formation of a hydrogel [START_REF] Wolf | Homogalacturonan Methyl-Esterification and Plant Development. 382![END_REF]. In addition, pectins are able to interact with other cell wall 112! components including cellulose and hemicellulose [START_REF] Popper | Xyloglucan-pectin linkages are formed intra-384! protoplasmically, contribute to wall-assembly, and remain stable in the cell wall[END_REF].

113!

An interesting feature is that newly synthesized homogalacturonans cannot spontaneously form 114! egg-box structures. Indeed, in order to trap divalent cations, the carboxyl groups of the 115! homogalacturonan chains must be negatively charged. When synthesized, carboxyl residues are 116! shielded by methyl groups, which must be removed to induce gelification. Therefore, two types 117! of pectin molecules can be distinguished: (i) "inactive" methyl-esterifed M-pectin and (ii) 118! "active", demethyl-esterified pectin, able to aggregate into macromolecular assemblies.

119!

Demethyl-esterification is carried on within the cell wall by specialized enzymes, the pectin 120! methyl-esterases (PMEs). In turn, PME-inhibitors (PMEIs) can modulate PME activity. This 121! process is essential, for both the degree and the pattern of demethyl-esterification influence the 122! mechanical/rheological properties of the constituted gel [START_REF] Peaucelle | Pectin-Induced Changes in Cell Wall Mechanics Underlie Organ 387! Initiation in Arabidopsis[END_REF]. Accordingly, modifications in the 123! levels of PME and PMEI expression significantly alter growth rates [START_REF] Peaucelle | Arabidopsis phyllotaxis is controlled by the methyl-esterification 389! status of cell-wall pectins[END_REF]. Note that other charged 124! residues such as acetyl-esters [START_REF] Caffall | The structure, function, and biosynthesis of plant cell 331! wall pectic polysaccharides[END_REF] can also be present on homogalacturonans and their regulation 125! could play a similar role. However, for the sake of simplicity we only consider methyl-ester 126! groups here. Noteworthy, 66 PMEs and 69 PMEIs have been reported in the Arabidopsis 127! genome, attesting for the importance of the process [START_REF] Wolf | Homogalacturonan Methyl-Esterification and Plant Development. 382![END_REF][START_REF] Sénéchal | Homogalacturonan-modifying enzymes: structure, expression, and 391! roles in plants[END_REF]. The PME/PMEI system appears as a 128! key regulator of cell wall expansion:

129! ! 8!
• First, by regulating the amount of "active" demethyl-esterified pectins, within the matrix 130! liquid phase, ready for insertion it influences its viscosity: The more active molecules 131! available, the faster the expansion. 132!

• Secondly, by removing more or less methyl groups on pectin molecules already stuck 133! together, it regulates the strength of the bonds within the matrix solid phase. This can be 134! interpreted as the tuning of the plastic yielding threshold of the cell wall.

135!

In conclusion, in the previous paragraphs we have seen that the cell wall can be considered as a 136! fiber reinforced hydrogel, which is put under tension by the osmotic pressure within the cells. To 137! grow, the cells have to yield to this pressure by expanding their walls in an irreversible manner.

138!

We will next discuss the chemical and physical processes that are at the basis of this expansion, 139! before indicating how molecular regulation can act to generate specific shapes of tissues and 140! organs. To be able to go beyond qualitative descriptions of growth regulation, we will present a 141! description in the form of a mathematical model. Finally we test this model using available 142! quantitative data 143!

Chemical equilibrium of the pectin hydrogel

144!

If single molecules in solution can bind to each other, macromolecular assemblies are generated.

145!

Assuming a high enough initial concentration of those molecules, such a mechanism eventually 146! leads to the formation of a hydrogel, i.e. the percolation of a macromolecular solid phase within a 147! solution of single molecules. The pectin matrix corresponds to such a biphasic system [START_REF] Iijima | Swelling behaviour of calcium pectin hydrogels by 333! thermomechanical analysis in water[END_REF]. When 148! the cell is growing, new matrix material has to be added: new components are exocytosed within 149! the liquid phase and components from the liquid phase have to be transferred into the solid one.

150!

Because the wall is a closed shell of constant thickness (i.e. with no border where an hypothetic 151! ! 9! nucleation mechanism could take place), a way to add new components to the solid phase is by 152! inserting them in between existing ones. Such a mechanism is sketched in Fig. 1a and formalized 153! in Eq.1, where ! !! and ! respectively depict one bond between two molecules in the solid phase 154! and a free molecule in the liquid one. 155!

! ! ! + !!! ⇄ !2! !! (1) 156!
This chemical equilibrium between the liquid and solid phase of the gel can be formalized by 157! standard concepts from thermodynamics. Indeed, phase equilibrium is reached once the free 158! energy (we will drop the "free" for now on) difference ∆! (= ! ! !"#$% ! -! !"#$!"% ) between the 159! two sides of Eq.1 vanishes. This energy difference reads

∆! = ! ! ! -! ! ! • !"!( ! /! !! ) and 160!
encompasses two general classic features. First, the stabilization of molecules by the adhesive 161! interactions between them is accounted for by the first term ! ! < 0 (the negative sign meaning 162! that adhesion stabilizes the system, i.e. lowers its energy). Second, the insertion of one molecule 163! from the liquid phase into the solid one comes with the loss of some degree of freedom. This is 164! expressed by the second term -! ! ! • !"!( ! /! !! ) a standard expression of mixing entropy, where 165! ! ! ! represents the thermal energy of the system (! ! and ! respectively stand for the Boltzmann 166! constant and the absolute temperature), ! depicts the monomer concentration in the liquid phase 167! and ! !! the water concentration in that same phase (≃ !55!"# • ! _ 1 ). In view of the above, the 168!

equilibrium condition ∆!! = !0 provides a value for a corresponding critical free molecule 169! concentration (with ! = 1 /! ! !): 170! ! ! * = ! ! ⋅ !"# !! ! ! (2) 171!
In other terms, the actual concentration c of monomers in solution controls the chemistry of the

172! ! 10! gel: 173! • If ! > ! ! * ! ⇔ ! ∆!! < 0, monomers spontaneously insert into the solid phase, growth 174! occurs. 175! • If ! < ! ! * ! ⇔ ! ∆! > 0, no spontaneous insertion, no growth occurs.

176!

The critical concentration ! ! * only depends on one parameter: the adhesive energy between 177! molecules in the solid phase (! !! in Eq.2). In order to modulate this equilibrium the cell only has 178! to tune the adhesive links between molecules. This can be achieved using enzymatic activities, 179! i.e. through molecular regulation. However, another efficient way to do so is by pulling on those 180! links (Fig. 1), as detailed in the following paragraph. 181!

Influence of mechanical forces on the chemical equilibrium

182!

A tensile force (!) applied on the solid phase can ease the breaking of existing bonds, thus 183! facilitating the insertion of soluble molecules. This shifts the chemical equilibrium towards 184! polymerizationi.e. towards the right-hand side of Eq.1. In this case, the energy difference has 185! to be lowered by the work !! = !! • ! generated by the tensile force (!) over a distance !, 186! characteristic of the deformation of loaded bonds:

∆! ! = ∆! -! • ! . Consequently, the 187! equilibrium condition ∆! ! = 0 yields, this time, a force-dependent expression for the equilibrium 188! concentration ! * : 189! ! * (!) ! = ! ! ! * • ! !!•!•! (3) 190!
This mechanism assures that tensile forces applied on a gel at equilibrium decrease its critical 191! concentration. As a consequence molecules in solution are recruited into the solid phase, (Fig. 1

& 192! ! 11!
2). This transduction mechanism is called force-induced polymerization and has been shown to 193! play an important role in the dynamics of major load-bearing structures such as the actin 194! cytoskeleton [START_REF] Kozlov | Processive capping by formin suggests a force-393! driven mechanism of actin polymerization[END_REF][START_REF] Courtemanche | Tension modulates actin filament polymerization mediated 395! by formin and profilin[END_REF][START_REF] Shemesh | Actin Polymerization upon Processive Capping by 397! Formin: A Model for Slowing and Acceleration[END_REF][START_REF] Shemesh | normalize concentration of 502! free active pectin molecules in the liquid phase. The orange curve is a visual representation of 503! Eq.3 and therefore represents the chemical equilibrium between the solid and the liquid phase of 504! the gel[END_REF]. Note that this mechanism is by essence irreversible, for when an increased 195! force is applied on a given configuration (black arrow on Fig. 2a) and then removed (dashed black 196! arrow on Fig. 2a) the system does not go back to its original position.

! ! ! = !! • !" • ! ! !•!!! ! (4) 205!
where the relative expansion rate (! ! ) is expressed as a function of the tensile force (!), 206! formalizing the out-of-equilibrium insertion reaction. ! > 0 is a proportionality constant and !" 207! represents the excess of free monomer in solution compared to the equilibrium concentration.

208!

Note that in an exhaustive approach the coefficient ! should depend on the elastic energy stored 209! in the solid phase. But in a turgid growing tissue this stored energy should be roughly constant 210! and can be ignored for the sake of simplicity.

211!

Reinforcing the gel: the cellulose microfibrils network steps in

212! ! 12!
Eq.4 depicts the expansion of a hydrogel under mechanical load, applied to the specific case of 213! pectin within the cell wall. A corollary question concerns the integration of the cellulose 214! microfibril network in this description.

215!

The entanglement between the pectin-based matrix and the cellulose fibrils leads to a complex 216! distribution of mechanical constraints at the molecular scale. An accurate mechanical analysis 217! would require a precise structural description of the considered wall and intensive numerical 218! simulations, similar to the ones produced by Kha et al [42] and Puri et al [43] in the case of cell 219! wall models containing only cellulose and hemi-cellulose fibers. Such a study lies outside of the 220! range of this article and we will only provide a qualitative analysis here.

221!

In a wall containing both pectin and cellulose, the force (!) felt by a pectin strand in a given 222! direction should in principle only be a fraction of the total force (! !"! ) exerted on a unit volume of 223! the wall. This fraction would depend, at the very least, on the relative concentration of cellulose 224! fibers compared to pectin strands ( ! !!! ) and the rigidity of those fibers ( ! ! ):

225! ! = !(! !"! , ! !!! , ! ! ) ≤ ! ! !"! .

226!

Controlling both the angular distribution of cellulose fibers and their rigidity appears as a way to 227! regulate the force exerted on pectin strands and consequently the expansion rate (! ! ) as 228! expressed in Eq.4. Following this idea, two stereotypical asymptotic cases can be evoked: 229!

• All the cellulose fibers are oriented perpendicular to the expanding pectin strand (c.f. Fig.

230!

1b-2). In this case, the cellulose fibers and the pectin strand can be considered in series 231! and the whole force exerted on the wall is transmitted to the strand: ! = ! !"! .

232!

• All the cellulose fibers are oriented in the same direction as the expanding pectin strand 233! ! 13! (c.f. Fig 1b-3). In this case, the cellulose fibers and the pectin strand lie in parallel and the 234! force felt by the pectin strand is smaller than the total one: ! ≈ ! !"! 1 + ! !!! ⋅ ! ! where 235! ! !!! and ! ! respectively account for the ratio of relative concentrations of pectin and 236! cellulose and the rigidity ratio between a pectin strand and a cellulose fiber.

237!

The qualitative point we want to make here is that enhancing the deposition of cellulose in a 238! given direction should directly down-regulate the expansion rate of the matrix in this specific 239! direction.

240!

Finally, it is worth mentioning that a detailed thermodynamical description of such a fiber-241! reinforced hydrogel should not only include dissipation between matrix molecules but also 242! between the fibers themselves and between the fibers and the matrix.

243!

Experiments on Chara fit the force-driven polymerization model

244!

We were not aware of quantitative experimental studies on higher plants focused on the 245! relationship between molecular composition of the cell wall, turgor-generated forces and the 246! morpho-dynamics of the whole tissue. Fortunately, such studies have been conducted on the 247! green algae Chara corallina [44,45] and we confronted the results with our model ( Fig. 3).

248!

In the case of isolated Chara cells, the mechanical force (!) exerted on the cell wall is 249! proportional to the pressure differential (∆!) between the inside and the outside of the wall. 250! Combined with Eq.4, this yields an exponential relationship between the pressure differential 251! exerted on a cell wall and its expansion rate. All other parameters kept constant, the ratio between 252! two relative growth rates ( ! ! /! ! ! ) as a function of the corresponding ratio of pressure 253!

differentials (∆!/∆! ! ) reads: 254! ! 14! ! ! ! ! ! = !"# C ! ⋅ ∆! ! ⋅ ∆! ∆! ! -1
(5) 255! see Fig. 3a for a comparison between experimental points and the fitting function.

256!

If we assume that the insertion of demethyl-esterified pectins is the molecular mechanism behind 257! cell wall expansion, the growing rate is an affine function of the concentration of demethyl-258! esterified pectin molecules in the solute phase. Fig. 3b depicts this behavior and qualitatively 259! matches experimental evidence of this proportionality exposed on Fig. 3

Discussion, perspectives & conclusion

269!

Our ambition was to demonstrate that a fundamental molecular mechanism (force-driven 270! polymerization) applied to a key actor (pectin) can be at the heart of a complex integrated 271! behavior such as morphogenesis.

272!

Although force-driven insertion appears as an appealing mechanism to explain cell wall 273! expansion at the molecular scale, it is worth noting that it might not be the only one possible. For 274! ! 15! instance, Rojas et al [START_REF] Rojas | Chemically Mediated Mechanical Expansion of the Pollen Tube 359! Cell Wall[END_REF] proposed a detailed study of pollen tube expansion using an apposition 275! mechanism: Newly demethyl-esterified pectin molecules are not inserted into the solid phase but 276! rather directly polymerized as a new layer next to the old one. This new layer modifies the 277! mechanical properties of the load-bearing cell wall and induces its viscous relaxation. Such a 278! mechanism does not exclude the model we propose here but is rather complementary: Rojas this mechanism provides a gateway to develop models that integrate quantitative geometrical 295! variables and biochemical equations (Fig. 4). In other words, it becomes possible to express 296! molecular activities in terms of geometrical outputs at the level of entire tissues and organs in a 297! ! 16! fully mechanistic manner. For this purpose, Eq.4 can be coupled to other differential equations 298! describing the dynamics of enzymatically-regulated variables. For instance one could assume that 299! the activation energy ! ! in Eq.4 depends on the demethyl-esterification degree of pectin 300! molecules and that pectin methyl esterases regulate its value. The influence of PME activity on 301! growth by coupling Eq.4 could then be explored using an equation describing PME dynamics.

302!

To assess the validity of the model a wide range of experiments can be considered. The first step 303! would be to extend the quantitative analysis produced by Proseus and Boyer to cells and tissues 304! from higher plants: quantitative tracking to follow cells and tissues submitted to various levels of 305! turgor pressure could be a start. In complement, rescue assays could be tested to restore a normal 306! growth phenotype to mutant or drug-treated plants by modulation of mechanical forces. Finally 307! and more speculatively, the use and analysis of biomimetic systems such as artificial walls could 308! also be very useful.

Glossary box:!

435!

Rheology: Study of the deformation of matter at a mesoscopic or macroscopic scale. 436! Viscosity: Rheological/mechanical property that relates the deformation rate of a material to its 437! mechanical load. For a given loading force, the more viscous a material is the slower it will 438! deform. Note that its inverse is called extensibility and is often used in growth modeling.

439!

! Thermal energy: The portion of the energy of a system due to the random relative microscopic 440! movements of its components. Absolute temperature of the system is, by definition, a direct 441! measure of this internal energy. For living systems which temperature is around 300 K the usual 442! value of their thermal energy is: The conservation equation of the number of bonds in the solid phase reads:

ℰ !! ∝ ! ! ! = 1 ! ≈ 4 • 10 !!" J, where ! ! = 1.38 ⋅ 10 !!" ! ⋅ 443! ! !!
449! !! ! !" = -! ! ! ! + ! ! (B1-1)

450!

The term of the left hand-side of Eq.B1-1 depicts the local time evolution of the bond 451! density (noted ! ! ). The first term of the right hand-side is proportional to the relative 452! expansion rate of the cell wall (! ! ) and accounts for the fact that the matrix can expand.

453!

The minus sign in front of it insures that spatial expansion of the matrix (! ! > 0) tends to

458!

Because the cell wall is constantly under tension, the spacing between molecules in the 459! solid phase and therefore the bond density can be assumed constant over time (i.e.

460!

!! ! !" = 0 in Eq.B1-1). This assumption implies that the two terms of the right hand-461! side of Eq.B1-1 compensate each other, leading to a direct relationship between the 462! expansion rate of the cell wall and the pectin chemical insertion rate:

463! ! ! ! ! = ! ! (B1-2)

464!

The proper derivation of chemical insertion rate (! ! ) from thermodynamical

103 ! 7 !

 1037 Pectins are produced in the Golgi and exocytosed within the cell wall. They are composed of four 104! different polysaccharides: Homogalacturonan, Rhamnogalacturonan I, Rhamnogalacturonan II 105! and Xylogalacturonan. Different pectin molecules correspond to various combinations of those 106! ! three residues covalently attached together. Hereby, homogalacturonans stand out, for they 107!

197 !

 197 Linking mechanics and biochemistry to geometry: spatial expansion198! of the pectin matrix under tensile force 199! One difference between the growth of living tissues and the plasticity of non-living materials is 200! the mass increase of the former. Therefore, it seemed logical to seek for a growth law based on 201! the equation of mass conservation of the solid phase of the matrix (see Box 1 for details). 202! Formalizing the fact that the local expansion of the matrix is fuelled by the insertion mechanism 203! previously exposed yields (combination of Eq. 2, 3 & B1-4): 204!

  of [45]: putting Chara 260! cells in a pectin-enriched solution increased their growing response to mechanical stress. 261! Another qualitative feature of Chara growth observed by Proseus et al and reproduced by our 262! model is its dependency on temperature: Fig 3c depicts how the pressure-growth rate relationship 263! is affected by a slight decrease of temperature (!), the shift from the plain curve to the dashed 264! one has to be compared to the shift between the two curves exposed on Fig.13 of [44]. These 265! qualitative agreements between our theoretical model and experimental data strengthen the idea 266! that force-driven assembly mechanisms of the cell wall matrix play an important role in cellular 267! growth. 268!

  et al 279! also describe a transduction mechanism between pectin chemistry and cell wall mechanics, 280! focusing their attention on the chemical control of the process. This shift of perspective between 281! their approach and ours (i.e. the question of causality between the chemical and the mechanical 282! states of the cell wall) is a molecular version of the debate between Linthillac and Schöpfer that 283! took place a few years back [46,47]: Is growth controlled by wall-loosening enzymes or by 284! turgor-induced stresses ? Our point is that what matters is the resulting equilibrium between both 285! ends which is formalized throught the critical concentration Eq.3. 286! Force-driven polymerization processes have been proposed to play a major part in several 287! biological functions relying on active gel dynamics and mechano-chemical transduction, from 288! actin filament formation and focal adhesion formation [9], [41], [48] (Kozlov et al) in animal 289! cells to cell wall instability and cell division in bacteria [49]. In our opinion, force-driven 290! polymerization processes could be as central in mechanobiology as Michaelis-Menten 291! mechanisms in biochemistry. 292! Indeed, force-driven polymerization can be seen as a transduction mechanism between a 293! mechanical force and a chemical reaction. From a theoretical perspective, in its full 3D version 294!

••••

  Cellulose microfibrils have historically drained most of the attention of biomechanical 422! studies of the cell wall. Recent results, however, suggest that the soft pectin matrix in 423! which they are embedded might also play a significant role in the physico-chemical 424! equilibrium of growing cells. 425! At the molecular scale, biophysicists have shown how mechanical forces applied on 426! molecular assemblies can modulate their chemical state and therefore initiate specific 427! biological responses. 428! At the tissular scale, numerical simulation tools borrowed from material sciences are 429! increasingly used in developmental biology. Morphogenesis in plant tissues is particularly 430!suited for this kind of approach. 431! Bridging the gap between the molecular and the tissular scales is a major challenge in 432! developmental biology, for it would help relating explicitly shape changes and specific 433!

Box 1 :

 1 and ! = 300!! respectively stand for the Boltzmann constant and the absolute temperature. 444! Work of a force: Considering an object in motion and summited to a force, the work developed 445! by this force corresponds to the scalar product between the force vector and the displacement 446From matter conservation to force-induced expansion 448!

  454! diminish the local bond density. Finally, the second term of the right hand-side (named 455! the source term and denoted ! ! ) grasps the chemical insertion mechanism that transfers 456! molecules from the liquid phase into the solid one when the two phases are out of 457! equilibrium.

  Figures legends

Figure 1 : 25 !Figure 2 :

 1252 Figure 1: Pectin strand dynamics & mechanical load distribution.

Figure 3 :

 3 Figure 3: Growth rate dependency on demethyl-esterified pectin

Figure 4 :

 4 Figure 4: Transducing chemical & mechanical inputs into

  Proceedings of the National Academy of Sciences102, 12383-12388 400! 42 Kha, H. et al. (2010) WallGen, Software to Construct Layered Cellulose-Hemicellulose 401! Networks and Predict Their Small Deformation Mechanics. Plant Physiology 152, 774-786 402! 43 Yi, H. and Puri, V.M. (2012) Architecture-Based Multiscale Computational Modeling of 403! Plant Cell Wall Mechanics to Examine the Hydrogen-Bonding Hypothesis of the Cell Wall

	404!
	405!	Network Structure Model. Plant Physiology 160, 1281-1292
	406!	44 Proseus, T.E. et al. (2000) Turgor, temperature and the growth of plant cells: using Chara
	407!	corallina as a model system. J. Exp. Bot. 51, 1481-1494
	408!	45 Proseus, T.E. and Boyer, J.S. (2007) Tension required for pectate chemistry to control
	409!	growth in Chara corallina. J. Exp. Bot. 58, 4283-4292
	410!	46 Wei, C. and Lintilhac, P.M. (2003) Loss of stability-a new model for stress relaxation in
	411!	plant cell walls. Journal of theoretical biology 224, 305-312
	412!	47 SCHOPFER, P. (2008) Is the loss of stability theory a realistic concept for stress relaxation-
	413!	mediated cell wall expansion during plant growth? Plant Physiology 147(3), 935-936 author
	414!	reply 936-938
	415!	48 Nicolas, A. et al. (2004) Cell mechanosensitivity controls the anisotropy of focal adhesions.
	416!	Proceedings of the National Academy of Sciences 101, 12520-12525
	417!	49 Jiang, H. et al. (2011) Mechanical Control of Bacterial Cell Shape. Biophysj 101, 327-335
	418!	!
	419!	!	!
	309!
		! !	19! 20!

Growth