J Martínez-Frutos

R Ortigosa

P Pedregal

F Periago

Robust optimal control of stochastic hyperelastic materials

Keywords: Robust optimal control, hyperelasticity, material uncertainty, turgor pressure, active fibers, soft robotics

Soft robots are highly nonlinear systems made of deformable materials such as elastomers, fluids and other soft matter, that often exhibit intrinsic uncertainty in their elastic responses under large strains due to microstructural inhomogeneity. These sources of uncertainty might cause a change in the dynamics of the system leading to a significant degree of complexity in its controllability. This issue poses theoretical and numerical challenges in the emerging field of optimal control of stochastic hyperelasticity. This paper states and solves the robust averaged control in stochastic hyperelasticity where the underlying state system corresponds to the minimization of a stochastic polyconvex strain energy function. Two bio-inspired optimal control problems under material uncertainty are addressed. The expected value of the L 2 -norm to a given target configuration is minimized to reduce the sensitivity of the spatial configuration to variations in the material parameters. The existence of optimal solutions for the robust averaged control problem is proved. Then the problem is solved numerically by using a gradientbased method. Two numerical experiments illustrate both the performance of the proposed method to ensure the robustness of the system and the significant differences that may occur when uncertainty is incorporated in this type of control problems.

Introduction

Since the early 1940s, the field of soft robotics continues exploring the creation of machines with biomimetic dexterous features superseding the capabilities of humans. This field represents a radically different paradigm with respect to the conventional field of hard robotics. The latter, which enabled an outstanding technological revolution in industrial automation, autonomous vehicles, etc. perform extraordinarily well for the specific tasks that they have been purposely designed for. The question that emerges then is therefore evident: if conventional robotics perform exceptionally, why to pursue a completely opposite robotic paradigm? Not to mention that soft robotics are made out of highly deformable materials such as elastomers, fluids and other soft matter, entailing a much higher degree of complexity in their controllability. The answer to this question lies on the programability of soft robotics for a wider range of tasks than their hard counterparts, along with their adaptability to rapidly changing uncertain conditions while performing these tasks [START_REF] Majidi | Soft-matter engineering for soft robotics[END_REF].

As already mentioned, the controllability of soft robots, potentially actuated by means of a wide spectrum of complex external stimuli (electric or magnetic field, mechanical pressure, osmotic pressure, etc.) is not a trivial task. This complexity reflects also in the numerical modelling of these soft materials, specially when aleatory sources of uncertainty such as random material properties, random geometry or random environmental conditions, are incorporated in the numerical model. This fact, together with the development of probabilistic uncertainty propagation methods, has fostered the interest in models represented by Partial Differential Equations (PDEs) with random inputs [START_REF] Babuska | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF][START_REF] Cohen | Approximation of high-dimensional parametric pdes[END_REF][START_REF] Gunzburger | Stochastic finite element methods for partial differential equations with random input data[END_REF]. In the context of stochastic hyperelasticity, large random variabilities observed in the bulk and shear moduli of Ogden-type hyperelastic materials, typically used to model elastomeric materials and soft biological tissues, are reported by recent works in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF][START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF]. Stochastic versions of the polyconvex energy functionals of these constitutive models have been proposed and validated against experimental data in those references.

The literature on optimal control of hyperelastic materials adopting a rigorous mathematical prism is relatively scarce even in the deterministic setting [START_REF] Günnel | Optimal control problems in finite strain elasticity by inner pressure and fiber tension[END_REF][START_REF] Lubkoll | An optimal control problem in polyconvex hyperelasticity[END_REF]. In these works it is wellaccepted to consider polyconvex strain energy functionals [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]11,[START_REF]Variational methods in nonlinear elasticity[END_REF] defining the constitutive model of the elastomer or soft material. Ball proves in his seminal paper [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF] that polyconvexity and coercivity (i.e. growth conditions) of the strain energy density entails the existence of minimisers of the total energy of the elastomer. Examples of isotropic polyconvex strain energy functions are the Odgen model, the Neo-Hookean model, the Mooney-Rivlin model, etc. The deterministic assumption made in most works concerning optimal control of highly deformable materials might be very convenient from the analysis and numerical standpoints. However, soft materials often exhibit intrinsic variability in their elastic responses under large strains leading to a significant degree of complexity in its controllability. Thus, accuracy and realism with respect to the physical response of a deformable soft material requires the incorporation of relevant sources of uncertainty in the numerical model. In this sense, a framework, based on information theory, for the modeling of stochastic Ogden class of stored energy functions is presented in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF][START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF]. The proposed models are almost surely consistent with the theory of nonlinear elasticity and are able to mimic the mean behavior and variability that are typically encountered in the experimental characterization of biological tissues.

The problem of optimal control of Stochastic PDEs has sparked significant interest last years. Existence results and, computationally efficient numerical resolution methods have been developed during the last two decades, but only for the case of random linear PDEs. In particular, in [START_REF] Zuazua | Averaged control[END_REF] the concept of averaged controllability was introduced and positive averaged controllability results have been recently obtained for parameter dependent evolution PDEs [START_REF] Lazar | Averaged control and observation of parameter-depending wave equations[END_REF][START_REF] Lü | Averaged controllability for random evolution partial differential equations[END_REF]. Although controlling the average of the state variable provides a first idea of robustness, it may be useless in scenarios where the dispersion of the system response is large. In the context of optimal control theory, the problem of controlling both the average and the variance of the state variable has been studied in [START_REF] Marín | Robust averaged control of vibrations for the bernoulli-euler beam equation[END_REF][START_REF] Martínez-Frutos | Robust optimal robin boundary control for the transient heat equation with random input data[END_REF][START_REF] Martínez-Frutos | Optimal control of PDEs under uncertainty. An introduction with application to optimal shape design of structures[END_REF] at the theoretical and numerical points of view. Significant differences between optimal controls minimizing only the average of the state variable and controls minimizing both the mean and the variance of the state variable were numerically observed in those references.

To the best of the author's knowledge, the optimal control of soft materials governed by stochastic hyperelastic constitutive models still remains an open issue. This paper addresses, both from the numerical and analysis standpoints, the robust optimal control of materials governed by stochastic hyperelastic constitutive models. A probabilistic framework is used to characterized the random parameters of a stochastic Ogden type stored energy function.

Inspired by [START_REF] Günnel | Optimal control problems in finite strain elasticity by inner pressure and fiber tension[END_REF], we consider two bio-inspired types of external actions (denoted as control in the sequel) on the elastomeric material, based on Turgor pressure and active fibers. As it will be shown in this paper, the first resembles the Turgor-based motion in the plant kingdom, whilst the second type of control simulates the electro-activation mechanism of muscles (e.g. the human heart). To take into account the variability of the deformed configuration with respect to desired target configuration, the deterministic cost functionals are averaged with respect to the underlying probability distributions.

The layout of the paper is as follows: Section 2 introduces the notions of hyperelasticity in the stochastic framework. In addition, it is shown how the two types of controls considered can be additively incorporated in the stochastic stored energy function of the material. Section 3 presents the optimal control problem in both deterministic and stochastic frameworks, showing the assumptions under which both control problems (associated with each type of control) are well-posed. Then, existence of optimal solutions for both problems is proved in Section 3.2. Section 4 describes the proposed, gradient-based, numerical resolution method and analyses two numerical experiments that show the importance of considering the uncertainty of the material parameters in the optimal control problem. Finally, Section 5 provides some concluding remarks and open problems.

Stochastic nonlinear continuum mechanics

Stochastic kinematics

Let D 0 ⊂ R N , N = 2, 3
, be an open, bounded and connected domain which represents the reference (or undeformed) configuration of an elastic body. Let (Ω, F, P) be a complete probability space, where Ω is the set of outcomes, F ⊂ 2 Ω is the σ-algebra of events, and P : F → [0, 1] is a probability measure. The random deformation of the body D 0 is defined through the mapping Φ : D 0 × Ω → R N , which, for every random even ω ∈ Ω, is assumed to be sufficiently smooth, injective and orientation preserving with respect to X ∈ D 0 . The mapping Φ links a material particle X ∈ D 0 to a particle x ∈ D (ω) ⊂ R N , the random deformed configuration, according to x = Φ (X, ω) (see Figure 1). As usual, u := x -X denotes the displacement field. Associated with the mapping Φ, the random deformation tensor F is defined as

F : D 0 × Ω → R N ×N , F = ∇ 0 Φ(X, ω),
where ∇ 0 (•) is the material gradient operator with respect to X ∈ D 0 . Associated with F , its co-factor H and its Jacobian J are defined as

H = (detF ) (F ) -T ; J = detF . (2.1)
The orientation-preserving condition is written as J (X, ω) > 0 a.e. X ∈ D 0 , and for every ω ∈ Ω. We assume that the boundary of D 0 is smooth and it is decomposed into two (also smooth) disjoint parts: Γ 0 D and Γ 0 N . On the Dirichlet boundary Γ 0 D , it is imposed Φ = Φ 0 for a given deformation Φ 0 : D 0 → R N , while on Γ 0 N , we impose a stress-free boundary condition. Notice that the above boundary conditions do not depend on ω ∈ Ω.

Stochastic polyconvex constitutive models

In this paper, we are interested in stochastic (also named random) versions of the classical deterministic polyconvex stored energy densities considerered in hyperelasticity theory. For the sake of clarity in the exposition, in the sequel, we focuss on the following compressible random Mooney-Rivlin model

W (ω, F ) = W(ω, F , H, J) = µ 1 (ω) ||F || 2 -3 + µ 2 (ω) ||H|| 2 -3 -2 (µ 1 (ω) + 2µ 2 (ω)) log J + µ 3 (ω) 2 (J -1) 2 , (2.2)
where the main novelty here is that the material parameters {µ 1 (ω), µ 2 (ω), µ 3 (ω)} are assumed to be random variables supported in R + * := {x ∈ R : x > 0}. As a consequence of the positivity of these material parameters, the density W (ω, F ) in (2.2) satisfies, for each fixed ω ∈ Ω, the following three crucial conditions (see [11,):

Polyconvexity: for the set V = {A, B, c} with A ∈ R N ×N , B ∈ R N ×N and c ∈ R + * , the following condition is satisfied W(ω, λV 1 + (1 -λ)V 2 ) ≤ λW(ω, V 1 ) + (1 -λ)W(ω, V 2 ), (2.3) ∀λ ∈ [0, 1] and ∀ V 1 , V 2 ∈ V. Coerciveness: W(ω, F , H, J) ≥ c 0 (ω) ||F || 2 + ||H|| 2 + J 2 + c 1 (ω); c 0 (ω) > 0, c 1 (ω) ∈ R. (2.4) Limit behavior 2 : lim J→0 + W(ω, F , H, J) = +∞. (2.5)
As in the deterministic situation, consistency with linearised elasticity in the stochastic framework entails the following relations between the above random material parameters and the shear µ(ω) and bulk κ(ω) moduli

µ 1 (ω) + µ 2 (ω) = 2µ(ω); 4µ 2 (ω) + µ 3 (ω) = κ(ω) - 2 3 µ(ω). (2.6)
Making use of the Maximum Entropy Principle, the authors in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF] show that the random material parameters {µ 1 (ω), µ 2 (ω), µ 3 (ω)} can then be given by

µ 1 (ω) = µ(ω)U 2 ; µ 2 (ω) = µ(ω)(1 -U ) 2 ; µ 3 (ω) = µ(ω) 2U - 8 3 + κ(ω), (2.7)
where {µ(ω), κ(ω)} are statistically independent random variables of Gamma-type, and U is a suitable Beta distributed random variable.

Bio-inspired controls

Fostered by [START_REF] Günnel | Optimal control problems in finite strain elasticity by inner pressure and fiber tension[END_REF], two types of bio-inspired controls are considered in the sequel. The effect of both controls can be modeled mathematically by additively decomposing the total stochastic stored energy density, denoted as W tot , as

W tot (ω, F , m) = W (ω, F ) + W m (X, F ); W m (m, F ) = m(X) Ŵm (F ), (2.8)
where the first term on the right-hand side of (2.8) corresponds to the purely mechanical contribution W in (2.2), and the second one, namely W m , represents the energy induced by the control m(X) : D 0 → R. Notice that there exists an implicit dependence of W m with respect to ω ∈ Ω through the deformation gradient tensor F = F (X, ω). The explicit expression of Ŵm (F ) for the two types of bio-inspired controls is presented next. (2.9)

This type of control resembles the turgor-driven movements of plants, characterised by efficiency and low power consumption. Based on this nature-inspired means of motion, research labs are in the process of developing new and biomimetic actuators [START_REF] Must | A variable-stiffness tendril-like soft robot based on reversible osmotic actuation[END_REF] (see Figure 2 a ).

It is shown next that the density W tot , as given by (2.8), enjoys the same polyconvex and coerciveness properties as the stored mechanical energy density W . Proposition 2.1. For each fixed m ∈ L 2 (D 0 ) and ω ∈ Ω, the total stored energy function (2.8) for the turgor pressure control is polyconvex and coercive.

Proof. The turgor-pressure contribution Ŵm (2.9) is linear with respect to the Jacobian J. Therefore, provided that W (ω, F ) is polyconvex, then W tot is polyconvex. With regards to coerciveness, from (2.4) it follows that

W tot (ω, F , m) ≥ c 0 (ω)(||F || 2 + ||H|| 2 + J 2 ) + c 1 (ω) -m(X)(J -1) ≥ c 0 (ω) 2 (||F || 2 + ||H|| 2 + J 2 ) + c 1 (ω) + c 0 (ω) 2 J 2 -m(X)(J -1).
(2.10)

The last two last terms in (2.10) are conveniently manipulated as

c 0 2 J 2 -m(X)(J -1) = c 0 2 J - m(X) c 0 2 + m(X) - m(X) 2 2c 0 ≥ m(X) - m(X) 2 2c 0 ,
which permits to conclude that

W tot (ω, F , m) ≥ c 0 (ω) 2 (||F || 2 + ||H|| 2 + J 2 ) + c 1 (X, ω); c 1 (X, ω) = c 1 (ω) + m(X) - m(X) 2 2c 0 (ω) ,
(2.11) and so W tot is coercive.

Tension fiber type control

In this case, the control exerts an inner fiber tension m(X) : D 0 → R along a unitary direction a(X) : D 0 → R N , with a(X) = 1, in the reference configuration. The energy Ŵm is then given by Ŵm (F ) = -

1 2 ||F a|| 2 .
(2.12)

There are relevant examples in nature where this type of actuation is harnessed as a means to achieve a desired motion, such as the complex electro-activation of the human heart, mathematically modelled as in (2.12). In that context, m represents the active cardiomyocyte contraction stress responsible for the deformation of the heart tissue (see Figure 2 b and [START_REF] García-Blanco | A new computational framework for electro-activation in cardiac mechanics[END_REF][START_REF] García-Blanco | Towards an efficient computational strategy for electro-activation in cardiac mechanics[END_REF]).

Proposition 2.2. Let µ 1 (ω) > 0 be the first parameter in the stochastic version of the Mooney-Rivlin model (2.2). Then, for each m ∈ L 2 (D 0 ) and ω ∈ Ω, with sup

X∈D 0 m(X) < 2µ 1 (ω), (2.13)
the total stored energy density function W tot in (2.8), with Ŵm given by (2.12), is polyconvex and coercive.

Proof. The energy density W tot can be re-arranged as

W tot (ω, F , m) = ψ 1 (ω, F ) + ψ 2 (ω, H, J); ψ 1 (ω, F ) = µ 1 (ω)||F || 2 - m 2 ||F a|| 2 ; ψ 2 (ω, H, J) = µ 2 (ω)||H|| 2 -2 (µ 1 (ω) + 2µ 2 (ω)) log J + µ 3 2 (J -1) 2 .
(2.14)

Clearly, the function ψ 2 (ω, H, J) is convex with respect to H and J. Therefore, polyconvexity of W tot requires convexity of ψ 1 (ω, F ) with respect to F . Let us define, for t ∈ R and F 0 , F ∈ R N ×N , the section

g(ω, t) = ψ 1 (ω, F 0 + tF ) = µ 1 (ω)||F 0 + tF || 2 - 1 2 m(X)|| (F 0 + tF ) a(X)|| 2 ,
The second derivative of g with respect to t yields

∂ 2 g ∂t 2 (ω, t) = 2µ 1 (ω)||F || 2 -m(X)||F a(X)|| 2 ≥ (2µ 1 (ω) -m(X)) ||F || 2 . (2.15)
It is therefore bounded from below provided that 2µ 1 (ω)-m(X) ≥ K > 0. As a result, ψ 1 (ω, F ) is convex with respect to F and the total contribution W tot , polyconvex. Coerciveness of the total contribution W tot is straightforward.

Stochastic state laws under the action of bio-inspired controls

Sections 2.3.1 and 2.3.2 have illustrated under what assumptions, stored energy functions of the type (2.8) with Ŵm , as those in (2.9) and (2.12), respectively, are polyconvex and coercive. As shown by Ball in his seminal paper [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF], this entails the existence of minimisers Φ(X, ω) for each m ∈ L 2 (D 0 ) and ω ∈ Ω for the total energy functional

Π (ω, Φ, m) := D 0 W tot (ω, F , m) dX, (2.16)
over the class of mapping fields

U := { Φ(ω) ∈ V D : H(ω) ∈ L 2 (D 0 ; R N ×N ), J(ω) ∈ L 2 (D 0
), and J(ω) > 0 a.e. in D 0 } ;

V D := Φ ∈ H 1 (D 0 ; R N ) : Φ = Φ 0 on Γ 0 D .
(2.17)

Remark 1. An additional term of the form

Π ext (Φ(X, ω)) = - D 0 f 0 (X, ω) • Φ(X, ω) dX - Γ 0 N t 0 (X, ω) • Φ(X, ω) dS, (2.18)
collecting the energy exerted by random, conservative volumetric f 0 and surface t 0 loads may be added to (2.16). For the sake of simplicity we do not consider such a term in what follows.

The stationary point Φ of the functional Π in (2.16) solves

∂Π ∂Φ (Φ(X, ω), m)(v) = D 0 P (X, ω) : ∇ 0 v dX = 0; v ∈ H 1 0 (D 0 ; R N ); (2.19) P (X, ω) : D 0 × Ω → R N ×N
being the random first Piola-Kirchhoff stress tensor, defined as

P (X, ω) = ∇ F W tot (ω, F , m) = ∇ F W (ω, F ) + ∇ F W m (F , m), (2.20)
where ∇ F W m (F , m) adopts the following expression for the Turgor-pressure and the tension fiber controls in (2.9) and (2.12), respectively,

∇ F W m (F , m) = -mH; ∇ F W m (F , m) = -m (F a ⊗ a) . (2.21) In equation (2.21), (a ⊗ b) ij = a i b j , ∀a, b ∈ R N .
Finally, integration by parts in (2.19) enables to obtain the strong or local form of the conservation of the linear momentum, i.e., DIV (P (X, ω)) = 0; in Ω 0 ;

(P (X, ω)) N = 0; on Γ 0 N ;

Φ = Φ 0 ; on Γ 0 D , (2.22)
where DIV(•) represents the material divergence operator, and N , the outward vector to the undeformed boundary Γ 0 .

Here is the issue of the potential non-uniqueness of solutions for (2.22). It is true that every minimizer in (2.16) is a solution of (2.22), but neither uniqueness of minimizers of (2.16) is guaranteed nor possible solutions of (2.22) that are not minimizers of (2.16). From an analytical point of view, we resolve this issue by considering, for each feasible m, all possible minimizers of (2.16). However, from the point of view of the numerical approximation, one has to work with (2.22), and non-uniqueness can be a real difficulty. It is nonetheless true that our numerical simulations of Section 4.2 have not shown any particular difficulty in this regard.

The control problem under uncertainty

The objective of the optimal control problem considered in this paper is two-fold: (a) to find a distribution of the control m(X) (Turgor pressure or fiber tension) that exerts a deformation on the undeformed domain D 0 , modifying its shape and deforming it into a desired target deformation Φ d ; (b) to minimise the effect of the uncertainty of the material parameters {µ i (ω), 1 ≤ i ≤ 3}, for the Mooney-Rivlin model (2.2) in the attainability of the target deformation.

Formulation of the optimal control problem

In this paper, we are mainly interested in the numerical approximation of robust optimal control problems for state laws of the form (2.22) . Thus, we introduce a discrete setting and so we assume that uncertainty in the material properties {µ 1 , µ 2 , µ 3 } is quantified by means of a multivariate and discrete random variable Z, with scenarios z i = (z i 1 , z i 2 , z i 3 ) ∈ R := R 3 + and probabilities π i , 1 ≤ i ≤ I. Continuous distributions may be recovered in the limit I → ∞ (see, e.g. [START_REF] Schrempf | Density approximation based on dirac mixtures with regard to nonlinear estimation and filtering[END_REF]).

The random variable Z induces a probability distribution P Z of the form

P Z = I i=1 π i δ z i , 0 < π i ≤ 1, I i=1 π i = 1,
where δ z i is a Dirac mass located at z i .

In the context of robust optimal control problems with continuous random variables and/or random fields, a similar discretisation process as the one proposed hereafter takes place (see [START_REF] Martínez-Frutos | Optimal control of PDEs under uncertainty. An introduction with application to optimal shape design of structures[END_REF]Chap. 4]). In this case, the robustness of the optimal control is achieved by incorporating a robustness criterion in the functional cost which typically involves the first two statistical moments of the performance. These statistical moments often require the numerical computation of multi-dimensional integrals that cannot be evaluated analytically. This issue has motivated the development of efficient methods to address the problem (such as dimension reduction methods and sparse grid collocation methods) which discretize the continuous random domain in a finite set of stochastic nodes z i and an associated weights π i .

For each scenario z i ∈ R, 1 ≤ i ≤ I, we denote by W i (F) the stored energy density

W i (F ) = z i 1 ||F || 2 -3 + z i 2 ||H|| 2 -3 -2 z i 1 + 2z i 2 (ω) log J + z i 3 2 (J -1) 2 , (3.1)
the total energy density W i tot (F , m) = W i (F ) + W m (X, F ), with W m (X, F ) given by (2.8), and finally the total energy functional

Π i (Φ, m) := D 0 W i tot (F , m) dX, Φ ∈ U. (3.2)
Let Γ 0 be a (measurable) part of the free boundary Γ 0 N in the undeformed configuration D 0 and let Φ d : Γ 0 → R N be a desired deformation. Then, for m ∈ H 1 (D 0 ), the following tracking-type in surface deformations functional is considered:

J (m) := I i=1 π i Γ 0 Φ i (X) -Φ d (X) 2 dS (X) + M 2 D 0 m(X) 2 dX + 2 D 0 |∇m(X)| 2 dX, (3.3
) where Φ i is a minimizer of (3.2). Eventually, the robust optimal control problem reads as:

Minimize in m ∈ H 1 (D 0 ) : J (m) subject to Φ i is a minimizer of (3.2), 1 ≤ i ≤ I. (3.4)

Existence of optimal controls

The goal of this section is to prove the following existence result:

Theorem 3.1. There is at least an optimal control m(X) for problem (3.4).

For the sake of clarity, this will be carried out in the deterministic case and later extended to the stochastic setting.

Deterministic case

Throughout this subsection it is assumed that I = 1.

Proof. We carefully describe all the steps necessary to prove Theorem 3.1.

(i) Let {m j } be a minimizing sequence for the control problem (3.4), and let {Φ j (X)} (Φ j ∈ U) be their corresponding mapping fields which are minimizers of the total energy functional (3.2), which throughout this subsection is denoted by Π (recall that I = 1). Let F j := ∇ 0 Φ j , and H j and J j be the co-factor and Jacobian of F j , respectively. The regularising part of the objective function J entails uniform boundedness of {m j } in H 1 (D 0 ). Hence, for a non-relabelled subsequence, m j m in H 1 (D 0 ) and m j → m in L 2 (D 0 ), for some feasible m. (ii) Take any feasible Φ ∈ U and let F := ∇ 0 Φ. Since Φ j is a minimiser of Π(m j ), then Π(Φ j , m j ) ≤ Π( Φ, m j ).

(3.5)

The uniform boundedness of {m j } in L 2 (D 0 ) implies the existence of a constant D independent of j which bounds the right-hand side of (3.5). Coerciveness of W tot then implies that Φ j is uniformly bounded in H 1 (D 0 ). Hence, by [11,, there exists Φ ∈ H 1 D 0 ; R N , such that for a not relabelled subsequence,

F j F , H j H, J j J in L 2 (D 0 ), (3.6)
where F := F (∇ 0 Φ), H := H(∇ 0 Φ), J := J(∇ 0 Φ). (iii) Inequality (3.5) may be rewritten as

Π(Φ j , m) + D 0 (m j -m) Ŵ (F j ) dX ≤ Π( Φ, m) + D 0 (m j -m) Ŵ ( F ) dX. (3.7)
Using the strong convergence m j → m in L 2 (D 0 ), the weak convergences (3.6), the polyconvexity of W tot and passing to the limit in (3.7) we conclude that

Π(Φ, m) ≤ Π( Φ, m) (3.8)
The arbitrariness of Φ implies that Φ is a minimiser of Π(m). By using Mazur's theorem, the compactness of the trace operator and following the same lines as in the proof of [11,, one proves that J > 0 a.e. in D 0 and that Φ = Φ 0 on Γ 0 D . Therefore, Φ ∈ U is a minimiser of Π for the limit control m. (iv) Finally, the weakly lower semicontinuity of the three contributions in the cost functional (3.3) implies that m is indeed a minimiser of problem (3.4).

Remark 2. Notice that the term involving the gradient of the control m(X) in J ensures the strong convergence m j → m, which is of crucial importance in the previous proof. The same conclusion may be obtained if lower and upper point-wise constraints on m, and boundedness of ∇m in L 2 (D 0 ) N are imposed.

Stochastic case

Proof. The proof follows along the same lines as in the deterministic case, given that we only have a finite number of possible scenarios. We highlight the specific ingredients brought by the probabilistic framework. Let {m j (X)} be minimizing for problem (3.4) and let {Φ i j (X}, 1 ≤ i ≤ I, be their corresponding mapping fields which are minimizers of the functional Π i in (3.2) for m = m j . Similarly to the deterministic case, m j m in H 1 (D 0 ) and m j → m in L 2 (D 0 ), for some m. From the proof of the deterministic case, it follows that for each 1 ≤ i ≤ I, there exists Φ i , weak limit in V D of (a subsequence of) Φ i j , that minimizes Π i in (3.2) for the limit control m(X). Since there is only a finite number of possibilities 1 ≤ i ≤ I, we can select a subsequence j, not relabelled, which is valid for all such i simultaneously, and

m j m, Φ i j Φ i , in H 1 (D 0 )
for all i. The result then follows taking into account that all the terms involved in the cost functional J (m), as given by (3.3), are weakly lower semicontinuous.

Numerical resolution method and numerical experiments

Gradient-based optimisation algorithm

In this paper, we advocate for a gradient-based optimisation method for the numerical approximation of the robust optimal control problem (3.4). As customary in this type of methods, in order to compute a descent direction, we use the standard Lagrangian method. To this end, let us consider the Lagrangian L defined as andP i is the first Piola-Kirchhoff stress tensors associated with scenario i, 1 ≤ i ≤ I, as given by (2.20) with ω replaced by z i and F by

L(Φ, p, m) = I i=1 π i Γ 0 Φ i (X) -Φ d (X) 2 dS (X) + M 2 D 0 m(X) 2 dX + 2 D 0 |∇m(X)| 2 dX - I i=1 D 0 P i : ∇ 0 p i dX, (4.1) where Φ = Φ 1 , • • • , Φ I , p = p 1 , • • • , p I ∈ H 1 D 0 ; R ,
F i = ∇ 0 Φ i .
Notice that in (4.1), (Φ, p, m) are considered as independent variables. Let (Φ, p, m) be a stationary point of L. The stationary condition of L with respect to p i coincides with the stationary condition of the functional Π (2.16) with respect to Φ (see equation (2.19)) for each scenario z i , 1 ≤ i ≤ I, i.e.

∂L ∂p i (Φ, p, m)(v) = - D 0 P i : ∇ 0 v dX = 0, ∀v ∈ H 1 0 (D 0 ). (4.2)
Equation (4.2) is nonlinear. A consistent linearisation of (4.2) is carried out by means of the standard Newton-Raphson method. The stationary condition of the Lagrangian L with respect to Φ i yields

D 0 ∇ 0 p i : C i : ∇ 0 v dX = I i=1 π i Γ 0 Φ i -Φ d • v dS (X) , ∀v ∈ H 1 0 (D 0 ), (4.3)
C i = C i (X, z i , m
) being the fourth order elasticity tensor, defined as

C i (X, z i , m) = ∇ 2 F F W (z i , F i ) + ∂ 2 F F W m (X, F i ), (4.4)
where the expression for ∂ 2 F F W m (X, F ) for both the Turgor-pressure and fiber tension control contributions are

∂ 2 F F W m (X, F ) = - m J (H iI H jJ -H iJ H jI ) ; ∂ 2 F F W m (X, F ) = -δ ij ma I a J , (4.5)
respectively, and where δ ij denotes the ijth component of the Kronecker Delta tensor.

Finally, the directional derivative of the Lagrangian L with respect to the control m in the direction m permits to obtain the descend direction. Precisely,

∂L ∂m (Φ, p, m)( m) = I i=1 D 0 ∂(∂ F Wm)(X,F i ) ∂m ( m) : ∇ 0 p i dX +M D 0 m m dX + ε D 0 ∇ 0 m • ∇ 0 m dX, (4.6)
where the expression of ∂(∂ F Wm(X,F )) ∂m ( m) for both the Turgor-pressure and fiber tension control contributions in (2.9) and (2.12) are, respectively

∂ (∂ F W m ) ∂m ( m) = -mH; ∂ (∂ F W m ) ∂m ( m) = -m (F a ⊗ a) . (4.7)

Numerical experiments

To keep track of the evolution of simulations during the optimisation procedure, we have introduced a sequence of discrete optimisation iterations akin to a pseudo-time parameter τ = {τ 0 , . . . , τ m }. At each discrete pseudo-time τ , the pseudo-time evolving control m(X, τ ) induces a random deformation on the undeformed configuration D 0 . For a random realization z i of the material parameters, the undeformed configuration is transformed into D

i (m(X, τ )) according to D i (m(X, τ )) = Φ i (m(X, τ ))(D 0 ),
where Φ i = Φ i (m) is a minimiser of (2.16) associated with m, as described earlier (see Figure 3). In the two numerical experiments that follow, for simplicity, the material parameter µ 2 in (2.2) has been taken igual to zero. µ 1 and µ 3 are modeled as two independent random variables following a Gamma distribution with means 1 and 3, and standard deviations 0.15 and 0.45, respectively. The continuous Gamma distributions are discretized by means of an isotropic sparse grid stochastic collocation method [START_REF] Martínez-Frutos | Robust optimal shape design for an elliptic pde under uncertainty in its input data[END_REF] leading to a finite set of scenarios z i depicted in Figure 4. Very similar cualitative results are obtained when the number stochastic nodes is increased.

Numerical experiment 1: Turgor pressure control

In this example, a Turgor-pressure type control (see equation (2.9)) is considered. Lower and upper bounds are imposed on the control variable, namely, -1 ≤ m(X) ≤ 1 for all X ∈ D 0 . The undeformed configuration D 0 is shown in Figure 5. The prismatic-shaped domain D 0 has lengths {L X , L y , L z } = {10, 3, 1} parallel to the {OX, OY, OZ} axis. The design domain is tessellated using a 24 × 4 × 2 regular mesh of linear hexahedral elements. D 0 is completely fixed in the surface corresponding to the minimum X coordinate. The image of the target configuration Φ d is also shown in Figure 5 (see black surface labeled as Target Surface (Spatial configuration)). Γ 0 is denoted by Material configuration in Figure 5.

Figure 6 represents the deformed configuration obtained by means of the deterministic control, i.e. taking the material parameters µ 1 = 1 and µ 3 = 3. In this case, the controllable part Γ 0 of the boundary perfectly matches the target configuration (both are overlapped) and Figures 7 a and 7 b show the mean deformed configuration for the deterministic control after computing all the deformed configurations associated with all the realizations z i . The value of the averaged objective function J in this case is J = 0.0193. Clearly, the mean deformed configuration of the Γ 0 is now considerably far from the target configuration. On the contrary, Figures 7 c and7 d show the mean deformed configuration for the robust control after computing all the deformed configurations associated with all the scenarios z i . In this case, the mean deformed configuration of Γ 0 is extremely close to the target (they seem overlapped in the Figure ), resulting in a value of the averaged cost functional of J = 0.0107. Crucially, this indicates that, as expected, the robust control (which takes into account the stochastic nature of the material parameters is less sensitive to dispersion in their random values. This is also corroborated by the lower value of the variance of the L 2 norm of the displacement field u := x -X in Figure 7 d (for robust control) and opposed to Figure 7 b (deterministic control).

Furthermore, Figure 8 shows the contour plot distribution of the control for the deterministic and robust cases from two different geometrical perspectives.

Numerical experiment 2: Active Fibers control

In this experiment, a fiber tension type control (see equation (2.12)) is considered. Fibers are oriented along the longitudinal direction (see Figure 11). The control variable is constrained according to 0 ≤ m(X) ≤ 1. The undeformed configuration D 0 , shown in Figure 9, is composed of four prismatic-shaped branches with lengths {L X , L y , L z } = {10, 3, 1} parallel to the {OX, OY, OZ} axis. Similarly to experiment 1, each branch of the domain is tessellated using a 24 × 4 × 2 regular mesh of linear hexahedral elements. The target surface (black surface labeled as Target Surface (Spatial configuration)) and zero Dirichlet boundary conditions (Fixed surface) are shown in Figure 9.

Figures 10 a and 10 b show the mean deformed configuration for the deterministic control after computing all the deformed configurations associated with all the realizations z i . The value of the averaged objective function J in this case is J = 0.0140. Clearly, the mean deformed configuration of the controllable portion Γ 0 is now considerably far from the target configuration. On the contrary, Figures 10 c and 10 d show the mean deformed configuration for the robust control after computing all the deformed configurations associated with all the realizations z i . In this case, the mean deformed configuration of Γ 0 is extremely close to the target, resulting in a value of the averaged cost functional of J = 0.0074. As in the previous example, but now with a different type of control, this indicates that the robust control is less sensitive to dispersion in the values of the material parameters. This is also corroborated by the lower value of the variance of the L 2 norm of the displacement field in Figure 10 d (for robust control) and opposed to Figure 10 b (deterministic control).

Finally, Figures 11 a and 11 b show the mean deformed configuration for the deterministic and robust controls, respectively after computing all the deformed configurations associated with all the realizations z i , seeing from a different geometrical perspective as in Figure 10. Notice also that the results have been rotated after applying rotations of {π/2, π, -π/2} with respect to the OZ axis, resulting in the four-legged actuator shown in Figure 11. In this Figure, the contour plot distribution of the control (in the deformed configuration, i.e., m (Φ -1 (x))) for both deterministic and robust settings is shown.

Conclusions and perspectives

This paper is a first step towards the very challenging problem of mathematical analysis and numerical approximation of robust optimal control problems constrained by the system of nonlinear elasticity with uncertainty in the material parameters. Precisely, we have analysed two types of bio-inspired controls: turgor pressure and fiber tension. Although the deterministic situation was numerically explored for the first time in [START_REF] Günnel | Optimal control problems in finite strain elasticity by inner pressure and fiber tension[END_REF], a rigorous mathematical proof of existence of optimal solutions for both problems in the deterministic and stochastic frameworks is presented in this manuscript. Although we have only considered a tracking-type in surface deformations cost functional, the theoretical results and numerical methods here exposed may be extended to other cost functional like the desired direction functional introduced in [START_REF] Günnel | Optimal control problems in finite strain elasticity by inner pressure and fiber tension[END_REF].

The analysis carried out in this paper focuses on the case of discrete distributions for the material parameters. A more sophisticated study is required in order to examine the more challenging cases of continuous distributions and random fields. A number of difficulties appear in the case of continuous randon variables and fields. Among them, the smooth selection of minimisers (in particular, measurability of the equilibrium configurations Φ (X, ω) with respect to the random parameter ω, which is required to be able to compute statistical moments of such configurations) and the existence itself of optimal solutions for the underlying optimal control problem.

Also, as is well-known in the case case of linear PDEs [START_REF] Babuska | A stochastic collocation method for elliptic partial differential equations with random input data[END_REF] convergence of stochastic collocation methods is a consequence of the regularity of the solution of the underlying PDE with respect to the random paramter. This regularity issue is a very challenging and open problem in the nonlinear situation considered in this paper. Nonetheless, the examples analysed in Subsection 4 show numerical evidence that convergence of stochastic collocation methods holds, as the consideration of finer discrete probability distributions do not introduce numerical instabilities.

Figure 1 :

 1 Figure 1: Graphical illustration of the random deformation mapping Φ (•, ω) : D 0 → D (ω) between reference D 0 and random deformed D (ω) configurations.

Figure 2 :

 2 Figure 2: (a) Tendril-like soft robot based on an osmotic actuation [19]. (b) Fiber structure of heart (courtesy of "The Correlation of 3D DT-MRI Fiber Disruption with Structural and Mechanical Degeneration in Porcine Myocardium")

Figure 3 :

 3 Figure 3: Undeformed configuration D 0 ; Deformed configurations at pseudo-time τ j for random realizations of z (blue dotted lines), reached after the application of the control m(τ )| τ =τj ; Target or desired configuration.

Figure 4 :

 4 Figure 4: Nodes of the considered isotropic sparse grid for uncertainty quantification in the material parameters µ 1 and µ 3 in (2.2).

Figure 5 :

 5 Figure 5: Numerical experiment 1: Problem configuration.

Figure 6 :

 6 Figure 6: Numerical example 1. Deformed configuration after application of the optuimal control in the deterministic setting. Contour plot distribution of the L 2 norm of the displacement field u := x -X.

Figure 7 :

 7 Figure 7: Numerical example 1. Mean deformed configuration after evaluating all the realizations z i of the material parameters for deterministic control (first row) and robust control (second row). The transparent configurations represent the deformed configurations associated with all the realizations z i . Control plot distribution of: mean of the L 2 norm of u for (a) deterministic and (c) robust control; variance of the L 2 norm of u for (b) deterministic and (d) robust control.

Figure 8 :

 8 Figure 8: Numerical example 2. First row: mean deformed configuration for deterministic control; Second row: mean deformed configuration for robust control. Top (a) and bottom (b) view for the deterministic control; Top (c) and bottom (d) view for the robust control.

Figure 9 :

 9 Figure 9: Numerical experiment 2: Problem configuration.

Figure 10 :

 10 Figure 10: Numerical example 2. Mean deformed configuration after evaluating all the realizations z i of the material parameters for deterministic control (first row) and robust control (second row). The transparent configurations represent the deformed configurations associated with all the realizations z i . Control plot distribution of: mean of the L 2 norm of u for (a) deterministic and (c) robust control; variance of the L 2 norm of u for (b) deterministic and (d) robust control.

Figure 11 :

 11 Figure 11: Numerical example 2. Mean deformed configuration for deterministic control; Second row: mean deformed configuration or robust control. Control plot distribution of the (a) deterministic and (b) robust control along the fibers in the deformed configuration.

The underlying justification for this condition is that it enforces the orientation-preserving condition, namely J(X, ω) > 0 a.e X ∈ D 0 , and hence the local invertibility and local injectivity of the mapping Φ.

Acknowledgements. This research was partially supported by the AEI/FEDER and UE under the contracts DPI2016-77538-R and MTM2017-83740 and by the Fundacion Séneca (Agencia de Ciencia y Tecnología de la Región de Murcia (Spain)) under the contract 20911/PI/18. The second author acknowledges the support provided by Ministerio de Ciencia, Innovación y Universidades, for the award of a Juan de la Cierva Formación Fellowship.


Szymon Greś

Michael Döhler

Palle Andersen

Laurent Mevel Context

Kalman filter-based stochastic subspace identification under mixed stochastic and periodic excitation

• Subspace-based system identification from output-only vibration measurements collected from structures in-operation -Modal analysis of civil, mechanical or aeronautical structures -Vibration modes are identified from eigenstructure of LTI system

• Intrinsic nature of the excitation may pose difficulties, e.g. presence of periodic inputs originating from rotating components of the structure, in addition to stochastic inputs -Identified eigenstructure contains both system and periodic modes -Consistency of the covariance-based subspace identification for measurements with oscillatory components showed in [START_REF] Favaro | Consistency of subspace methods for signals with almost-periodic components[END_REF] -In practice periodic modes often disturb the estimation of close structural modes AIMS • Reconstruction of output signal where the periodic part is removed as preprocessing for engineering applications

• Identification of the eigenstructure of the stochastic system part only MODELING

• System states are physical quantities of mechanical system i.e. displacements and velocities

• Periodic states represent the periodic excitation u(t)

Mechanical model

Mz(t) + C ż(t) + Kz(t) = f (t) + bu(t) y(t) = Lz(t) + v(t) u(t) = h i=1 e i sin(ω i t + g i ) A sys :λ sys i , ϕ sys i → f sys i , ζ sys i , ϕ sys i A per :λ per i , ϕ per i → f per i , ζ per i , ϕ per i ⇒ ⇐ Continuous-time combined state-space model ẋsys (t) = A sys c x sys (t) + bu(t) + w(t) y(t) = C sys x sys (t) + du(t) + v(t) ⇓ Discrete-time combined state-space model x sys k+1 x per k+1 = A sys A b 0 A per =A x sys k x per k + w k 0 y k = [C sys C per ] =C x sys k x per k + v k
• Proposed modeling is equivalent to [START_REF] Favaro | Consistency of subspace methods for signals with almost-periodic components[END_REF] up to similarity transform

• Identification of both system and periodic modes with the UPC algorithm [START_REF] Van Overschee | Subspace Identification for Linear Systems[END_REF]:

Left: Discrete-time poles of a 6 DOF chain system, Right: Estimated frequencies f sys i ,f per i and damping ratios ζ sys i , ζ per i of 6 DOF chain system under combined random and periodic excitation.

REMOVAL OF THE PERIODIC SUBSIGNAL BY ORTHOGONAL PROJECTION AND SYSTEM IDENTIFICATION

• The goal is to reconstruct responses where the periodic signal is discarded and then to identify the observability matrix of the structural system: Γ sys Γ sys = (C sys ) T (C sys A sys ) T . . . (C sys A sys p ) T .

Algorithm:

• Prediction of the periodic subsignal with the non steady-state Kalman filter Similarity transform of the innovation state-space model into modal basis

xV k+1 = A V xV k + K V k e k , y k = C V xV k + e k , with V = [ (Ψ) (Ψ)]
, where Ψ = ϕ sys 1 . . . ϕ sys m ϕ per 1 . . . ϕ per h , and

xV k = V -1 xk , A V = V -1 AV , C V = CV , K V k = V -1 K k .
This allows to select states corresponding to the periodic modes with a selection matrix S and subsequently approximate the periodic subsignal

ŷper k = C V S xV k .
• Reconstruction of the output system response by projection of the raw output measurements onto the orthogonal complement of the periodic subsignal estimate

The orthogonal projection of the raw data Hankel matrix Y raw onto the data Hankel matrix of the predicted periodic subsignal Y per yields the decomposition • Frequencies and damping ratios estimated for model orders 10-40 and compared to estimates from random vibration experiment:

Y - pro = Y - raw /Y - per ⊥ = Γ sys Z -+ KE -+ E - K , Y + pro = Y + raw /Y + per ⊥ = Γ sys Z + + KE + + E + K . Γ sys can
Left: Raw measurements, Right: Reconstructed measurements.

CONCLUSIONS

• Extension of the output-only stochastic and periodic state-space modeling to the UPC algorithm • Consistent algorithm to reconstruct the stochastic part of the system response to mixed random and periodic excitation • In practical application parameters estimated from the reconstructed measurements are closer to estimates from a comparable random vibration experiment

  be factorized from the projection similar to UPC

	Y + pro /Y -pro = Γ sys Z .
	NUMERICAL VALIDATION
	• Comparison of the exact and estimated states and system response of 6 DOF chain system
	under mixed periodic and random excitation:
	Top: Periodic states, Bottom: System response.
	APPLICATION: MODAL ANALYSIS OF A PLATE
	• Rectangular aluminum plate excited with a shaker
	• Vibration measurements with 16 sensors at 4096 Hz
	with 491,520 samples
	• Two experiments:
	-Random vibration
	-Random vibration with added sinusoidal signal of
	370 Hz induced by the shaker
	• Stochastic subspace identification and
	orthogonal projection of predicted periodic subsignal

• Stabilization diagrams of modal frequencies with output measurements spectra:

Left: Raw measurements, Right: Reconstructed measurements.


E Farge

M Bitbol

P E Devaux

Biomembrane elastic response to intercalation of amphiphiles

Keywords:

A model of the elastic behavior of a biomembrane in response to intercalation of amphiphiles into the bilayer is developed. This model takes into account the bilayer couple hypothesis (Sheetz and Singer 1974), and assumes that incorporation of amphiphiles into one layer of the membane exerts mechanical work on the elastic biomembrane. The model accounts for an apparent experimental discordance noted by several authors: the variation in area observed upon incorporating amphiphiles is smaller by a factor of about 2 than the variation expected using previous models.

Introduction

In good qualitative agreement with the bilayer couple hypothesis [START_REF] Sheetz | Biological membranes as bilayer couples. A molecular mechanism of drug-erythrocyte interactions[END_REF] and geometrical analysis of red blood cell shape transformations [START_REF] Beck | Relations between membrane monolayers in some red blood cell shape transformations[END_REF], intercalation of drugs or exogenous phospholipids in the outer leaflet of the erythrocyte bilayer results in echinocytosis [START_REF] Lange | Interaction of cholesterol and lysophosphatidylcholine in determining red cell shape[END_REF][START_REF] Kuypers | Shape changes in human erythrocytes induced by replacement of the native phosphatidylcholine with species containing various fatty acids[END_REF][START_REF] Ferrell | Membrane bilayer balance and erythrocyte shape: a quantitative assessment[END_REF][START_REF] Allan | Echinocytosis and microvesiculation of human erythrocytes induced by insertion of merocyanine 540 into the outer membrane leaflet[END_REF]. Conversely, depletion of the outer layer, or intercalation of amphiphiles in the inner layer, produces stomatocytosis [START_REF] Isomaa | Shape transformations induced by amphiphiles in erythrocytes[END_REF][START_REF] Rosso | Influence of chlorpromazine on the transverse mobility of phospholipids in the human erythrocyte membrane: relation to shape changes[END_REF][START_REF] Haest | Selective removal of lipids from the outer membrane layer of human erythrocytes without hemolysis. Consequences for bilayer stability and cell shape[END_REF]). However, quantitative assessments have shown that the expansion (or the constriction) of the treated layer is not as large as one would expect taking into account the number of molecules incorporated (or removed), and the average area per molecule [START_REF] Ferrell | Membrane bilayer balance and erythrocyte shape: a quantitative assessment[END_REF]Haest etal. 1981;[START_REF] Allan | Echinocytosis and microvesiculation of human erythrocytes induced by insertion of merocyanine 540 into the outer membrane leaflet[END_REF]). The authors observed variations in surface area which were half the expected values. Several explanations have been put forward in the above quoted papers to explain this differ-Offprint requests to. E. Farge ence. Some of these explanations speculated that cholesterol redistributed across the bilayer. Other explanations concern the mechanical properties of the membrane. In the present paper, we propose a simple mechanical analysis of the membrane elastic response to intercalation (or depletion) of amphiphiles into the bilayer, which accounts for the numerical discrepancy between observed and predicted variations of area.

Theory

We consider that the membrane behaves as if it were made of two coupled layers. For two coupled layers, bending (or shape change) is associated with the compression or dilation of each leaflet. So each layer is described as a two-dimensional elastic surface, with its own set of elastic moduli: we shall only make use of the elastic area compressibility moduli Ki (K1 for the inner layer, and K 2 for the outer layer). For the sake of simplicity, we shall focus our attention on the situation wherein amphiphiles (drugs, phospholipids or lysophospholipids) are intercalated in the outer layer. The results obtained in this case will easily be transposed to the other situation where amphiphiles are depleted from the outer layer.

The mechanical work needed to vary the area of each elastic surface can be written as follows (Evans 1974):

Wl = 1/2K~e2; W 2 = I/2K 2e22 (1)
where

S 1 -Sol $2 -S02 ~1 - So 1 and a2- So 2 (2)
S01 and So2 are respectively the areas of the inner and the outer leaflet in their unconstrained state, whereas S 1 and S 2 are the actual areas of same leaflets in any given state.

Intercalation or depletion of amphiphiles in the outer layer of the bilayer provides an external mechanical work dW ex~, producing the bending of the membrane and its shape change (echinocytosis in the case we are interested in). This work is due to the change of the unconstrained outer area 8o2 and can be written:

dWeXt ~ W2 -~S~o2 ~; dS°2' (3)
received by the two The mechanical work dW rec leaflets expresses the resistance of the membrane to bending. It is due to compression (or extension) of their surface relative to the unconstrained area, and is the following:

dW rec 0W1 dcq ~W2 = e~, So, + ~ So~ da2. ( 4
)
The bending of the bilayer as a consequence of the unconstrained area So2 variation is a purely mechanical and adiabatic transformation. Therefore, the energy conservation leads to a conservation equation for the mechanical work wherein external given work must be equal to received work:

0W2 8W1 de 1 + ~ W2 d~ 2 :- dS02 ~0~1 Sot ()(X2 SO2 ~802 ~t2 or, taking (1) into account: dS02 K1 ~1 dczl + K2 °~2 d~2 ---K2 ~2(~z2 + 1) So~ (5)
One may go one step further, by noticing that, since the two layers are coupled, extension of So2 leads the membrane to a bending equilibrium state. At this bending equilibrium, the stresses of each layer must be equal:

Kt ~1 = K2 °~2 or K1 d~l = K2 d~2 • Then, (5) becomes: dc~ 2 - K1 dS°2 (6) e2 + 1 K l + K 2 802
Finally, taking (2) into account, (6) can be written:

dS2 _ K: dSoa (7) $2 K2+K1 So2
Relative variations of the area are on the order of a few percents: this justifies the infinitesimal notation employed in (7). Formula (7) shows that the relative variation of the outer area dS2/$2 is not equal to the relative variation of the corresponding unconstrained area dSoz/Soz. In the situation we are interested in, the variation of the unconstrained area actually equals the average area per molecule times the number of incorporated, or extracted, molecules. Thus, we find that, taking into account the elastic properties of both layers of the membrane, the variation of the outer layer area is necessarily smaller than the average area per (added or extracted) molecule times the number of molecules.

Moreover, if K 1 = K 2, we have dS2/S 2 = 1/2 dSo2/So2, which is in good agreement with the measurements performed by [START_REF] Ferrell | Membrane bilayer balance and erythrocyte shape: a quantitative assessment[END_REF] and by [START_REF] Allan | Echinocytosis and microvesiculation of human erythrocytes induced by insertion of merocyanine 540 into the outer membrane leaflet[END_REF].

Discussion

Several explanations have been put forward in order to account for the difference between the experimental evaluations of the outer leaflet expansion, and the theoretical estimates based on the cross-sectional area of the intercalated molecules. We discuss hereafter two major proposals:

The cholesterol hypothesis

It has been suggested [START_REF] Haest | Selective removal of lipids from the outer membrane layer of human erythrocytes without hemolysis. Consequences for bilayer stability and cell shape[END_REF]) that cholesterol translocation could compensate for the imbalance induced by intercalation (or depletion) of amphiphiles in the outer leaflet. This model is based on the assumption that the transbilayer movement of cholesterol is very fast (see [START_REF] Lange | The rate of transmembrane movement of cholesterol in the human erythrocyte[END_REF]. But these data have recently been challenged [START_REF] Van Meer | Plasma membrane cholesterol pools[END_REF]. Moreover, a direct study of the effect of cholesterol on lysolecithin-induced red cell crenation [START_REF] Lange | Interaction of cholesterol and lysophosphatidylcholine in determining red cell shape[END_REF] showed that the amount of lysolecithin required for echinocytosis decreases as the membrane cholesterol content increases.

One would expect the opposite if cholesterol redistribution accross the bilayer would occur, and act as a partial compensatory mechanism for the lysolecithin-induced excess of the outer leaflet area: the less cholesterol present in the membrane, the less efficient should be the compensation.

The dilation hypothesis

The other main explanation is purely mechanical, involving an expansion of the inner lipid leaflet in response to the outer leaflet expansion [START_REF] Allan | Echinocytosis and microvesiculation of human erythrocytes induced by insertion of merocyanine 540 into the outer membrane leaflet[END_REF]. But here, one must realize that the outer leaflet expansion following intercalation of amphiphiles is related to an increase of the unconstrained area So2, whereas the actual area $2 is likely to increase less than So2. As for the hypothetical inner leaflet expansion, it can be related to an increase of the actual area $1 with respect to a fixed Sol. In this case, one would have (~2 <0 and ~1 :>0, and thus, the tensions in the two leaflets would have opposite directions. Hence, a torque would persist, and the bilayer would not be in mechanical equilibrium. Bending would result, and lead the membrane to an equilibrium state with both layers compressed.

Intuitively, if the surface does not extend -or diminish -as much as it should, elastic compression -or dilation -may effectively occur, and one must take into accout the elastic properties of the biomembrane. We have developed a model of elastic response of a biomembrane to amphiphiles intercalation, which resolves quantitatively experimental discordances in the above cited publications. The present model allows us to make some considerations on an hypothetic mechanical asymmetry in biological membrane, on its consequences and its origin.

The bilayer elastic response

In our model, the bilayer couple hypothesis is expressed by the equilibrium condition K 1% =K 2 e2. This condition expresses equality between the tensions T~ and T 2 of the two layers, which is the definition of mechanical equilibrium between two coupled layers related by shape changes: this is the bending equilibrium state. Then we write the external mechanical work received by the bilayer because of amphiphiles incorporation into a leaflet. Writing the energy conservation equation as an equality between received and external work, and introducing the bilayer couple hypothesis leads us to show that when 802 is increased by intercalation of amphiphiles, then both e, and ~2 are negative, namely both leaflets are compressed with respect to their unconstrained area. Depletion of the outer leaflet and decrease of 802 would rather lead to e~ and % positive and to an expansion with respect to the unconstrained area.

This model allows us to give a quantitative answer to the publications mentioning a discrepancy of about 2 between d8o2/8o2 (the expected value of area variation) and dS2/$2 (the observed value). Assuming K 1 =/£2, and introducing this condition in (7), one finds: dS2/$2 = (1/2) dSo2/So2, which is the result obtained in the experiments, account being taken of the inaccuracies.

A putative mechanical asymmetry

Moreover, if one considers the fact that in the quoted experiments the factor is not exactly equal to 2, our model provides the possibility of estimating a putative mechanical asymmetry of compressibility moduli across the membrane. The asymmetric transverse repartition of the phospholipids is already known to involve physical asymmetric properties across the bilayer: this is the case for the lateral diffusion of phospholipids in the membrane which is faster by about a factor 5 in the inner layer than in the outer layer of the red blood cell membrane (For a review see [START_REF] Zachowski | Bilayer asymmetry and lipid transport across biomembranes[END_REF]. Such an asymmetry could be the cause of differences in the elastic physical properties between the two layers. It has been shown in model membranes that there is a close correlation between phospholipid composition and elastic compression moduli [START_REF] Evans | Physical properties of surfactant bilayer membranes: thermal transitions, elasticity, rigidity, cohesion and colloidal interactions[END_REF]. Each layer, characterized by its own phospholipid composition, could consequently have its own elastic behavior, and thus different elastic moduli K1 and K 2 .

Deuling and Helfrich have pointed out that such an asymmetry may play a role in the spontaneous curvature of the membrane [START_REF] Deuling | The curvature energy of fluid membranes: a catalogue of vesicule shapes[END_REF]. Also it has been shown theoretically, with a bending energy minimisation method [START_REF] Svetina | Membrane bending energy and shape determination of pbospholipid vesicles and red blood cells[END_REF], that the value of the difference between external area of the two layers can be correlated to some red blood cell shape changes, like discocyte or cup formation. A similar hypothesis has been put forward to explain vesicle shape changes [START_REF] Sackmann | Membrane bending elasticity and its role for the shape fluctuations and the shape transformations of cells and vesicles[END_REF]).

Thus, one could imagine that mechanical asymmetry could control changes in membrane curvatures, and pos-sibly play a role in phenomena involving local shape changes which are involved in the early stage of endocytosis. In conclusion, modulation of the difference in area between the inner and outer layers of eukaryotic cell membranes via the putative aminophospholipid translocase may play an important physiological role [START_REF] Devaux | The aminophospholipid translocase: a transmembrane lipid pump; physiological significance[END_REF].

The mechanical asymmetry could either be purely related to the properties of the lipid bilayer or involve the cytoskeleton. As an instance of the second case, one can consider an experiment [START_REF] Haest | Stabilisation of erythrocyte shape by a chemical increase in membrane shear stiffness[END_REF]) wherein the addition of diamide cross-links the spectrin network and decreases its compliance, increasing the K1 value. In this situation, the effective area variation under amphiphile depletion of the outer layer (dS2/$2) tends to 0. Using formula (7) with K1 ~K2 leads to the same result. Our model, taking into account the elastic asymmetry is, thus, in good agreement with this diamide experiment too, where the cytoskeleton is involved as the internal layer. Noticing that echinocytic shape change can be conceptualised as an expansion of the bilayer in respect to the cytoskeleton, properties of the cytoskeleton may also be involved in mechanical asymmetry for such a transformation.

A consequence of our model is that shape changes induced by incorporation or depletion of amphiphiles respectively lead the membrane to a compressed or dilated equilibrium state. Consequently, the surface tension of the membrane should be non-zero. A non-zero tension is known to inhibit flickering phenomenon, the thermal fluctuations of the membrane [START_REF] Brochard | Frequency spectrum of the flicker phenomenon in erythrocytes[END_REF]. According to our model, shape change should then imply inhibition of flickering. However, because of the very low value of the bending elastic modulus -about ]0 -2° N -m -tensions due to bending are probably to small to involve any detectable inhibition of flickering. This could be the reason why the discocyte-echinocyte shape transformation does not experimentally exhibit any inhibition of the flickering [START_REF] Zeman | Bending ondulations and elasticity of the erythrocyte membrane: effects of cell shape and membrane organisation[END_REF]).

Acknowledgements. This work was supported by grants from the Centre National de la Recherche Scientifique (UA 526) and the Universit6 Paris VII. E. E is a recipient of a followship from the Ministare de la Recherche et de la Technologic.


Farge

P F Devaux

L Mathivet

S Cribier

E Farge

How lipid asymmetry can make vesicles fusion-competent by inhibition of the thermal undulations

How lipid asymmetry can make vesicles fusion-competent by inhibition of the thermal undulations

Introduction

The first step in any fusion event involving biological membranes or liposomes is the close contact between the two lipid bilayers. Hefore the local collapse which gives rise to a single membrane from two membranes, the two phospholipid interfaces have to approach each other. Different interacting forces are involved. Van der Waals forces are attractive and can extend to distance of the order of 10 nm. Electrostatic forces between charged lipid head groups belonging to different membranes are generally repulsive, but the high dielectric constant of water and the screening by ions can minimize this repulsion. Moreover the lateral diffusion of phospholipids enables the charged lipids to slide away from the contact points. The so-called hydration forces are responsible for a potential barrier with a thickness of about 2-3 nm, depending on the nature of the head groups; aminophospholipids, being less hydrated than choline-containing phospholipids, are more favourable for fusion [ 1,2].

In addition to these terms, Helfrich's group have shown that the thermal undulations of lipid bilayers correspond to a steric interaction, which can be accounted for by a repulsive force between membranes, with a typical length of interaction of the order of 10 nm [3]. Thus, van der Waals attractive forces and steric repulsions due to thermal undulations are the most important and opposing terms in fusion events and to adhere, membranes have to overcome this steric repulsion.

The thermal undulations correspond to the superposition of surface fluctuations with different wavelengths. The amplitude of these waves gives an indication of the distance of steric repulsion. With a light microscope very long waves with a large amplitude are visible, but typical amplitudes are usually below the resolution of a light microscope.

The thermal undulations can be inhibited in a variety of ways, for example by reducing the temperature or by stretching the membranes. Indeed, phospholipid membranes undulate because the Abbreviations used: a, areal strain; dS, surface elements; AN, difference in number of lipids in both layers; BSA, bovine serum albumin; EPC, egg phosphatidylcholine; h, thickness of bilayer; No, average number of lipids per leaflet; Ro, average radius. lateral tension is normally approximately zero [4].

Inducing a lateral tension can inhibit the fluctuations and hence permits membrane cohesion. Servuss and Helfrich have shown that cohesion can be 'turned on' by an osmotically-induced lateral tension [5]. Rending the membranes is also a way of reducing the undulations: firstly because very long wavelengths can no longer exist, thus the frequency spectrum is reduced and secondly because bending itself can induce lateral tension.

In this paper we demonstrate that an asymmetrical distribution of phospholipids, which is generally associated with membrane curvature, can favour membrane adhesion, probably by inhibition of the thermal fluctuations.

Theory

In the following we will consider only phospholipid bilayers made of a single phospholipid species, so that the areal asymmetry is associated with a proportional difference in the number of lipids between leaflets: AN/No, where A N is the difference in number of lipids in both layers and No the average number of lipids per leaflet. It has been shown previously that the transition from a unilamellar liposome with a discoid shape, characterized by an average radius &, to a liposome made of two connected spheres of approximately the same radius, necessitates a change in lipid asymmetry, where AN, and ANin,, are the difference in the number of lipids in both layers in the final state and in the initial state, respectively, h is the thickness of the bilayer and Iz a constant equal to 0.83 [6]. This transition from a discoid to two spheres is sometimes called a 'vesiculation', although it does not imply the separation into two independent vesicles. From (l), it can be deduced that the quantity h/Ro is a scaling factor for the shape changes reduced by a compositional asymmetry. For giant liposomes R, is of the order of 10 p m and h is 8 nm, thus the percentage of lipid transfer necessary for a shape change is very small (of the order of 0.1%). For the so-called large unilamellar vesicles (LUV) with a radius of about 45 nm, the percentage of lipid trans-fer required for similar shape change is 15% and it would be more than 30% for sonicated vesicles with a radius of 20 nm. Thus in the two latter cases a redistribution of a few per cent of the lipids does not trigger any significant shape change. On the other hand, the lipid asymmetry builds up a lateral tension faster in small vesicles than in giant vesicles.

The lateral tension of a liposome which is not under osmotic stress and which has a large radius is very small. An expansion or contraction of the surface area, SS, associated with the areal strain,

6S a = -S"'

can generate a tension z. When considering the inhibition of the undulations, the relationship between z and a has been established by Helfrich and Servuss; it is not a simple linear relation

[3].
The value of a must be approx.

or larger to inhibit the undulations partially.

The asymmetrical insertion (or depletion) of molecules is equivalent to the addition (or depletion) of a surface element 6s' from one monolayer.

We have shown that, because of the coupling between the two elastic monolayers, the actual surface element incorporated (or subtracted) is

6S= 112 6S* [7], so that 6S* 2 s a =--.
The transfer of a surface element from one monolayer to the other can be decomposed into two steps: (1) its depletion from one monolayer, which builds up a tension z; (2) its insertion into the opposite monolayer, creating an opposite tension which cancels out the tension created by (1). However if the membrane has an average local radius of curvature R,, the relative surface area taken from one side is different from the relative surface area added to the other side. It has been shown [6] that the areal strain a is related to the lipid asymmetry by the following relation In this equation h/R, appears also as a scaling factor. For giant liposomes, h/R, z 4 x and a z 5 x lo-', a value too small to inhibit the thermal undulation [4]. On the other hand for large unilamellar vesicles, a lipid translocation of a few per cent generates a tension sufficient to block the undulations.

The Experimental section gives a few examples of inhibition of thermal undulation accompanied by vesicle aggregation.

277

Experimental section

Giant liposomes of egg phosphatidylcholine (EI'C) with a diameter in the range 10-50 p m were obtained by two different procedures. The first procedure was described in [8] and involves incubating a lipid film overnight in water, at 45"C, under an argon atmosphere. This technique gives an inhomogeneous population of liposomes, some of which are quasi-spherical, others elongated with a (a)

discoid shape. The unilamellarity can be verified by the micropipette technique of Kwok and Evans [9].

In the second procedure, the lipid film is submitted to very low frequency electric field for approx. 2 h, using the technique described by Angelova et al.

[ lo]. The latter technique generates unilamellar spherical liposomes with a rather homogeneous size distribution. The liposomes were used either immediately or after several days of incubation at 4°C.

This storage period was accompanied by partial lipid degradation with formation of lysophosphatidylcholine and with some lipid peroxidation. Addition of fatty acid-free bovine serum albumin (HSA) had no effect on freshly-prepared liposomes, but generated important modifications in the 'aged liposomes, as indicated below.

Discoid liposomes

We will first review briefly previous experiments from this laboratory [ 81. Lysophosphatidylcholine was injected through a micropipette in the vicinity of a non-spherical giant liposome. This liposome underwent a sequence of shape change from discoid to two connected spheres (see Figure 5 (a)-(c) in [8]), followed rapidly by the adhesion of the two spheres over a large fraction of their surface. The same result was obtained recently with the charged molecule lysophosphatidylglycerol (data not shown). We interpret these observations by assuming that the addition of the lyso-compound with a micropipette, which injects a large quantity of molecules onto the outer monolayer of the liposome, creates a surface tension large enough to block the undulations. The two vesicles composing the giant liposome then attract each other by van 2 of reference [8]. In the latter case, the fraction of lipid which is reoriented is too small to induce a significant surface tension. Figure l a shows several giant liposomes with a discoid shape. These liposomes were stored for several days and then BSA was injected through a micropipette in the vicinity of these liposomes. The same liposomes, are shown in Fig. 16 and are now aggregated. BSA, by depleting the outer monolayer of lyso-compounds and oxidized lipids, not only triggers shape changes [8], but also, by inducing lipid asymmetry, generates a surface tension that blocks the thermal undulations. Note that the vesicles become 'fusion competent', but do not fuse.

Spherical liposomes

Figure 2 shows a populaton of spherical EPC liposomes obtained by the electric-field technique (the second procedure described previously). It is not known why this technique generates unilamellar liposomes that are so perfectly spherical. BSA added to these liposomes when freshly-prepared had no effect. After several days, the undulations became clearly visible as indicated by the less spherical pattern of the vesicles in Figure 3a. At this stage, HSA provokes aggregation of the vesicles to form typical honeycomb regions (Figure 3b). Again the most probable interpretation is that a fraction of the lipids is degraded during the storage, so that BSA extracts degraded lipids from the outer mono-

The same population of liposomes as in Figure 2, after 5 days storage ( 0 ) before and (b) after addition of BSA

Observation was with a Normarsky contrast interdifferential (Zeiss inverted microscope) (4 layer selectively and creates an imbalance between the inner and outer monolayers.

Discussion 279 The results presented above show that the selective addition or removal of lipids from one monolayer can be accompanied, not only by shape changes, but also by vesicle adhesion. In fact in the case of spherical vesicles the shape changes are rarely observable. In the experiments reported here, only lipid insertion or lipid depletion efficiently created a tension. The redistribution of lipids from one monolayer to the other did not cause adhesion. W e have explained in the Theory section why lipid reorientation in giant liposomes with a typical diameter of 10-30 p m should not be expected to lead to adhesion. Indeed equation ( 3 ) reveals that a scaling factor is involved, so that the redistribution of phospholipids necessary to allow membrane adhesion varies with the size of the vesicle.

With vesicles the size of a cell (several pm), phospholipid redistribution driven, for example, by a lipid pump such as the aminophospholipid translocase [ 121 triggers essentially shape changes. However with organelles the size of endosomes or granules (about 100 nm) the aminophospholipid translocase is more likely to induce or maintain a surface tension than to trigger shape changes. In turn, surface tension will inhibit surface undulations making these small vesicles fusion competent.

Hydrated bilayers of phospholipid normally repel each other strongly, as a consequence of hydration repulsion arising from solvation of the polar head groups. When CaL+ ions are added to phosphatidylserine, however, the surface water that gives rise to the repulsive forces is completely displaced by the bound CaL+ and the phospholipid bilayers In common with most cell types, erythrocytes do not exhibit spontaneous cell-fusion reactions. This would seem to be consistent with the topological distribution of phosphatidylserine in their plasma membranes, since the phospholipids in the outer leaflet of the human-erythrocyte plasma membrane are primarily phosphatidylcholine and sphingomyelin and the inner leaflet contains the acidic phospholipids [S]. It is therefore important to know if there is a requirement for phosphatidylserine to be on the outside surfaces of cells for cell fusion to occur and, in my laboratory, the human erythrocyte has been employed in work on this question. W e have investigated whether changes in phospholipid asymmetry occur in human erythrocytes that are induced to fuse artifically by two different protocols: (1) osmotic swelling in the presence of Ca'+ and (2) exposure of the cells to electrical-breakdown pulses in the presence of several different divalent cations.

Abbreviation used: PEG, poly(ethy1ene glycol)

Exposure of phosphatidylserine is associated with osmotically-induced cell fusion

In our initial work we investigated the fusion of human erythrocytes that occurs when monolayers of the cells are induced to swell in the presence of 10 mM Ca2+ by the entry of small molecules such as malonamide or small poly(ethy1ene glycol)s, for example, PEG 400. Under these conditions a proportion of the cells fuse before they lyse, whereas in the absence of Caz+, lysis occurs without any cell fusion 161. Exposure of phosphatidylserine on the outside of platelets [7, 81 and sickled erythrocytes [ 91 yields procoagulant surfaces that facilitate the conversion of prothrombin into thrombin. The catalytic potential of such surfaces has been considered to provide a semi-quantitative measure of the exposed phosphatidylserine [7, 81. Using this assay we found that the cell fusion which accompanied cell swelling induced by the entry of small permeant molecules in the presence of CaL+ was associated with the exposure of phosphatidylserine in the outer leaflet of the plasma membrane of erythrocytes [ 101.

It has been observed now that the fusion of monolayers of human erythrocytes which accompanies cell swelling induced by the entry of malonamide or of PEG 300 is preceded by hemi-fusion (that is, membrane fusion without cytoplasmic fusion). This was demonstrated by the rapid movement of a fluorescent membrane probe (carbocyanine) from unlabelled to labelled cells [ 111.

Unlike complete cell fusion, hemi-fusion occurred both in the presence and in the absence of added Caz+, but was inhibited by the addition of 10 mM EGTA. W e also found that hemi-fusion (and hence complete fusion) in these experiments was prevented strikingly by pre-incubation of the cell

Figure I

 I Figure I Suspension of EPC liposomes seen by phase contrast microscopy with an inverted optical microscope (Leitz)The bars correspond to 20 pm. ( a ) Before addition of BSA: ( b )after addition of BSA.

Figure 2 Spherical

 2 Figure 2 Spherical unilamellar liposomes obtained by submitting a lipid film of EPC to a low-frequency electric field for 2 h The bar corresponds t o 20 p m . Observation was with a Normarsky contrast interdifferential (Zeiss inverted microscope)

  achieve molecular contact [ 1, 21. Caz+ ions are much more effective than other cations (for example, Mg'+) in inducing membrane fusion in model membrane systems based on phosphatidylserine [ 31 because the corresponding complex of phosphatidylserine with MgL+ does not lose its layer of hydration [4].

The work described here was supported by grants from the Centre National de la Kecherche Scientifique (UA526), the Institut National de la SantC et de la Kecherche Medicale (CRE 90 0104), the liniversitt-Paris VII and the Direction des Kecherches, Etudes et Techniques of the LXItgation GtnCrale pour I'Armement (9 1.079); L.M. holds a fellowship from the 1)ClCgation GCnCrale pour I'Armement.