Anne-Florence Bitbol

Luca Peliti

Jean-Baptiste Fournier

Membrane stress tensor in the presence of lipid density and composition inhomogeneities

We derive the expression of the stress tensor for one and two-component lipid membranes with density and composition inhomogeneities. We first express the membrane stress tensor as a function of the free-energy density by means of the principle of virtual work. We then apply this general result to a monolayer model which is shown to be a local version of the area-difference elasticity (ADE) model. The resulting stress tensor expression generalizes the one associated with the Helfrich model, and can be specialized to obtain the one associated with the ADE model. Our stress tensor directly gives the force exchanged through a boundary in a monolayer with density and composition inhomogeneities. Besides, it yields the force density, which is also directly obtained in covariant formalism. We apply our results to study the forces induced in a membrane by a local perturbation.

Introduction

Biological membranes consist principally of a lipid bilayer. In each monolayer, the molecules are arranged with their hydrophobic ends directed inward toward the second monolayer, and their hydrophilic ends directed outward toward the surrounding water [START_REF] Alberts | Molecular Biology of the Cell[END_REF][START_REF] Mouritsen | Life-as a matter of fat (The frontiers collection[END_REF]. Membranes are fluid and possess a weak resistance to bending, so they fluctuate strongly under the effect of thermal excitations [START_REF] Helfrich | [END_REF]. In many relevant situations, their total area may be treated as a constant: this is due to the existence of a practically constant area per lipid and to the fact that lipids are virtually insoluble in water [START_REF] Mouritsen | Life-as a matter of fat (The frontiers collection[END_REF]. As a consequence, membranes develop a lateral tension in response to external forces [4,5].

Although the constant area approximation works finely for many fundamental aspects (see, e.g., Refs. [4,[6][7][8][9][10][11][12]), there are important situations in which it is essential to take into account the variability of the lipid density in each monolayer. Indeed, an essential parameter in the phase diagram of vesicle shapes is the difference between the areas of the two monolayers [14][15][16]18], and the most accurate equilibrium description of the elasticity of membranes thus involves a monolayer "area difference elasticity" (ADE) term in addition to the usual bending terms [17][18][START_REF] Seifert | The Structure and Conformation of Amphiphilic Membranes[END_REF]. There are also spectacular instabilities that may be triggered by a local expansion of one monolayer, yielding cristaelike invagination [21], ejection of tubules [22] or bursting and curling of polymersomes [23]. Finally, the dynamics of membranes at the sub-micron scale is controlled by the a Author for correspondance.

intermonolayer friction resulting from local velocity differences between the monolayers fluids that must be treated as compressible [24,25].

In this paper, we derive a fundamental tool that will help to rationalize these phenomena: the stress tensor acting in a membrane monolayer of variable shape, variable lipid density and variable composition. Understanding forces in complex membranes with various degrees of freedom is crucial to understand their equilibrium shape and their dynamics [26,27]. The divergence of the stress tensor gives the density of elastic forces in the membrane, which is the basis of a dynamical description. It is very useful in the study of shape instabilities [28]. Moreover, the stress tensor provides the forces directly exchanged through a boundary. This is a valuable information that can be used, e.g., in the calculation of membrane-mediated interactions [29,30], or in the study of the spatial distribution of forces on the edge of an object embedded in the membrane.

The outline of our paper is the following: In Sec. 2, we present a formal derivation of our stress tensor from the principle of virtual work. In Sec. 3, following Ref. [25], we construct a free-energy density for one and two-component monolayers that extends that of the ADE model. In Sec. 4, we derive the stress tensor associated with this free-energy density, and in Sec. 5 we discuss the associated force density. Finally, in Sec. 6, we show how our results enable to understand the forces and the dynamics in the case of a local perturbation of a membrane by the microinjection of a reagent.

Stress tensor formal derivation

Let us consider one monolayer of the bilayer. We describe its shape in cartesian coordinates by the equation z = h(x, y) of a surface S parallel to its hydrophobic interface with the other monolayer. Such a description, based on the height of the membrane with respect to a reference plane, is often referred to as the Monge gauge. We do not assume yet that the membrane is weakly deformed. Let ρ(x, y) be the projected mass density, i.e., the lipid mass per unit area of the reference plane (x, y). In order to study the case of a two-component monolayer, let us denote by φ(x, y) the local mass fraction of one of the two lipid species, say species number 1. The case of a onecomponent monolayer can be obtained by setting φ = 0.

Let f (ρ, φ, h i , h ij ) be the projected free-energy density of the monolayer (i.e., the free energy per unit area of the reference plane). Here and in the following, Latin indices represent either x or y (not z) and h i ≡ ∂ i h, h ij ≡ ∂ i ∂ j h, etc. Note that we are assuming that the free energy depends only on the mass density, on the local lipid composition, and on the slope and curvature of the monolayer. We have thus neglected the gradients of the curvature and of the lipid density and composition. The former approximation has already proved successful [START_REF] Helfrich | [END_REF]18]; the latter is justified by the fact that the correlation length of the density fluctuations should not be larger than the monolayer thickness (far from a critical point).

Let us consider an infinitesimal cut with length ds separating a region A from a region B in the monolayer, and let us denote by m the normal to the projected cut directed toward region A. The projected stress tensor Σ relates linearly the force df that region A exerts onto region B to the vectorial length m ds of the projected cut through df = Σ • m ds .

(

This defines the six components of the projected stress tensor: Σ ij and Σ zj , where i ∈ {x, y} and j ∈ {x, y} [32].

To determine the projected stress tensor, we shall follow the method presented in Ref. [32], which is based on the principle of virtual work. Let us consider a monolayer patch standing above a domain Ω of the reference plane. This patch is supposed to be a closed system with fixed total mass of each lipid species. Its free energy reads

F = ∫ Ω d 2 r f (ρ, φ, h i , h ij ) .
(

) 2
In order to deal with arbitrarily deformed states of the monolayer patch, we assume that in addition to the boundary forces (and torques) exerted by the rest of the monolayer, the patch is submitted to a surface density w(x, y) of external forces, and to individual external forces acting on the molecules and deriving from a one-body potential energy v α (x, y) for the lipid species α ∈ {1, 2}. The former forces control the shape of the patch and the latter control the mass density distribution of the lipids within the patch. The total potential energy corresponding to these latter forces can be written as

E p = ∫ Ω dn 1 v 1 + dn 2 v 2 = ∫ Ω d 2 r [ ρφ v 1 µ 1 + ρ(1 -φ) v 2 µ 2
] ,

(3) where d 2 r = dx dy, and µ α denotes the mass of one lipid of the species α.

Introducing v = v 2 /µ 2 and u = v 1 /µ 1 - v 2 /µ 2 ,
this potential energy can be rewritten as:

E p = ∫ Ω d 2 r [ρv + ρφu] . ( 4
)
At equilibrium, the lipid density ρ, the composition φ and the shape h of the monolayer are controlled by the external actions represented by u(x, y), v(x, y) and w(x, y).

Let us study a small deformation of the monolayer patch at equilibrium: Ω → Ω +δΩ, h → h+δh, ρ → ρ+δ ρ and φ → φ + δφ. Each element of the patch, initially at position (x, y), undergoes a displacement δa(x, y), with δa z = δh + h k δa k [32]. The variation of the free energy of the monolayer patch during the deformation reads

δF = ∫ Ω d 2 r [ ∂ f ∂ ρ δ ρ + ∂ f ∂φ δφ + ∂ f ∂h i δh i + ∂ f ∂h ij δh ij ] + ∫ δΩ d 2 r f , (5)
We now perform two integrations by parts, and we use the relation

∫ δΩ d 2 r = ∫ B ds m i δa i , ( 6
)
where B denotes the boundary of Ω. Assuming that the translation of the monolayer edges is performed at a fixed orientation of its normal, so that δh j = -h jk δa k along the boundary [32], we obtain

δF = ∫ Ω d 2 r [ ∂ f ∂ ρ δ ρ + ∂ f ∂φ δφ + δF δh δh ] + ∫ B ds m i { f δa i + [ ∂ f ∂h i -∂ j ∂ f ∂h ij ] δa z (7) + [( ∂ j ∂ f ∂h ij - ∂ f ∂h i ) h k - ∂ f ∂h ij h jk ] δa k } , where δF δh = ∂ k ∂ j ∂ f ∂h jk -∂ j ∂ f ∂h j . ( 8
)
The potential energy variation during the deformation is

δE p = ∫ Ω d 2 r [(v + uφ) δ ρ + uρ δφ] + ∫ B ds m i δa i [v ρ + uρφ] . ( 9
)
The total variation δF + δE p of the energy of the system must be balanced by the work δW done by the surface force density w and by the boundary forces exerted by the rest of the membrane on our patch. Since the translation of the monolayer edges is performed at a fixed orientation of its normal, the torques produce no work. We may write

δW = ∫ Ω d 2 r [w k δa k + w z δa z ] + ∫ B ds m i [Σ ki δa k + Σ zi δa z ] . ( 10
)
As the monolayer patch is considered as a closed system, the total mass of each lipid species in the patch is constant during our deformation:

∫ Ω d 2 r ρ = M and ∫ Ω d 2 r ρφ = M 1 , ( 11
)
where M and M 1 are constants. Let us introduce two constant Lagrange multipliers λ and µ to implement these two global constraints. The relation δF + δE p -δW + λδM + µδM 1 = 0 must hold for any infinitesimal deformation of the monolayer patch. The identification of bulk terms in this relation yields

w z = δF δh and w k = -h k w z , ( 12
)
∂ f ∂ ρ = -(v + λ) -(u + µ)φ , ( 13
)
∂ f ∂φ = -(u + µ)ρ . ( 14
)
By identifying the boundary terms and using (13), we obtain the components of the membrane stress tensor:

Σ ij = ( f - ρ ∂ f ∂ ρ ) δ ij - ( ∂ f ∂h j -∂ k ∂ f ∂h kj ) h i - ∂ f ∂h kj h ki , ( 15
)
Σ zj = ∂ f ∂h j -∂ k ∂ f ∂h kj , ( 16
)
which generalizes the result of Ref. [32] to the case where there are inhomogeneities in ρ and φ. Note that the fraction of each lipid species does not appear explicitly in this result. Therefore, Eqs. ( 15)-( 16) hold both for onecomponent monolayers and for two-component monolayers. Note also that Σ does not depend directly on the external actions, which confirms its intrinsic nature.

Comparing our result with Ref. [32] shows that taking into account lipid density variations only changes the isotropic term of the stress tensor, which now reads fρ ∂ f /∂ ρ. This term is reminiscent of minus the pressure of a two-dimensional homogeneous fluid binary mixture with area A described by a free energy F (T, A, N 1 , N 2 ) = Af (T, ρ, φ):

-P = ∂F ∂A T,N1,N2 = f -ρ ∂f ∂ρ T,φ . ( 17
)
While the last expression can be used locally in a nonhomogeneous fluid mixture, the case of the membrane is more complex since its free energy depends on the curvature. The "surface pressure" in a membrane is sometimes defined as ρ ∂g/∂ρ -g where g is the part of the membrane free-energy density f that depends only on ρ [34] or, equivalently, what remains of f for a planar membrane [35]. Interestingly, we find that the isotropic part of the membrane stress tensor does not identify to minus this surface pressure, since it is the complete, curvaturedependent projected free-energy density f that appears in f -ρ ∂ f /∂ ρ.

The divergence of the stress tensor gives the force per unit area, p, exerted by the rest of the monolayer on the patch. By direct differentiation, we obtain:

p z = ∂ j Σ zj = - δF δh , ( 18
)
p i = ∂ j Σ ij = δF δh h i -ρ ∂ i ∂ f ∂ ρ + ∂ f ∂φ ∂ i φ, ( 19
)
where we have used

∂ i f = (∂ f /∂ ρ)∂ i ρ + (∂ f /∂φ)∂ i φ + (∂ f /∂h j )h ij + (∂ f /∂h jk )h ijk .
At equilibrium, we can use ( 12)-( 14) to express p, which yields

p z = -w z , ( 20
) p i = -w i + ρ∂ i v + ρφ∂ i u . ( 21
)
These relations constitute the balance of surface force densities for the monolayer at equilibrium. In particular, p vanishes at equilibrium when the membrane is submitted to no external actions (i.e. w = 0 and u = v = 0).

Monolayer model

Free-energy density in terms of local variables

One-component monolayer

Let us derive the elastic free energy of a monolayer in a bilayer from basic principles, first for a one-component monolayer. We shall recover and extend the model of Ref. [25], which is actually a local version of the ADE model. We assume that the free-energy density f per unit area of the monolayer depends only on the mass density ρ and on the local principal curvatures c 1 and c 2 of this monolayer. As in the previous section, the gradients of the curvature and of the density are neglected in our description. Note that, unlike f and ρ, f and ρ are the free energy and the mass per actual unit area of the monolayer, and not per projected unit area. We use for both monolayers the density and the principal curvatures defined on the same surface S of the bilayer, so that the curvatures are common to the two monolayers.

We will consider the physically relevant regime of curvature radii much larger than the membrane thickness. We will also restrict ourselves to small variations of the density around a reference density ρ 0 . Note that it can be convenient to take ρ 0 different from the equilibrium density ρ eq of a plane monolayer with fixed total mass, for instance to study a monolayer under tension. Let us define

r = ρ -ρ 0 ρ 0 = O ( ) , (22)
H = (c 1 + c 2 ) e = O ( ) , ( 23
) K = c 1 c 2 e 2 = O ( 2 ) , ( 24
)
where e is a small length in the nanometer range that allows to define the scaled total curvature H and the scaled Gaussian curvature K. Since we typically expect 10 -4 ≤ |r| ≤ 10 -2 and 10 -4 ≤ |c i e| ≤ 10 -2 , it is sensible to assume that r and H are O(

) while K is O ( 2 )
. The free-energy density f is a function of these three nondimensional small variables.

To study small deformations, we write a second-order expansion of f :

f (r, H, K) = σ 0 + A 1 H + A 2 (r -H) 2 + A 3 H 2 + A 4 K + O ( 3 ) . ( 25
)
Three comments are due here. i)We have not included any term linear in r in this expansion. Indeed, the total mass of the monolayer, i.e. the integral of ρ = ρ 0 (1 + r), is assumed to be constant, so including a term linear in r is equivalent to redefining the constant term σ 0 . ii) The freedom associated with the choice of e allows to set the coefficient of -rH equal to twice that of r 2 . iii) All the coefficients in (25) depend on the reference density ρ 0 . We shall come back in the following on the constitutive relation σ 0 (ρ 0 ). Defining the constants κ, κ, c 0 and k through

A 1 = - κc 0 2e , A 2 = 1 2 k, A 3 = κ 4e 2 , A 4 = κ 2e 2 , ( 26
)
and setting c = c 1 +c 2 , we obtain the following expression, which generalizes those of Refs. [25,35]:

f = σ 0 + k 2 (r -ec) 2 - κc 0 2 c + κ 4 c 2 + κ 2 c 1 c 2 + O ( 3 ) . ( 27
)
Note that all the above terms have the same order of magnitude. Indeed, since typically

k ≈ 10 -1 J/m 2 , κ ≈ κ ≈ 10 -19 J, e ≈ 1 nm and c -1 0 ≈ 50 nm (see, e.g., Refs. [2, 13, 37]), we have A 2 ≈ A 3 ≈ A 4 ≈ 100 A 1 .
The advantage of the procedure we have employed is that we control precisely the order of the expansion.

In ( 27), e can be interpreted as the distance between the surface S where c 1 , c 2 and r are defined and the neutral surface of the monolayer [START_REF] Safran | Statistical Thermodynamics of Surfaces, Interfaces and Membranes[END_REF][START_REF] Petrov | [END_REF]. As a matter of fact, the density on a surface parallel to S can be expressed as a function of the distance between S and this surface as 2 , where the sign depends on the orientation. We choose the minus sign here, keeping in mind that the second monolayer then has the plus sign (see figure 1). Let us now consider the surface such that = e. If f is written as a function of r(e) ≡ r n and of the curvatures, it features no coupling between these variables. This corresponds to the definition of the neutral surface [START_REF] Safran | Statistical Thermodynamics of Surfaces, Interfaces and Membranes[END_REF], which means that e is the distance between S Fig. 1. Schematic drawing of a lipid bilayer. The principal curvatures c1, c2 and the scaled densities r ± of both monolayers are defined on S. The distances between S and the neutral surfaces N ± of monolayers ± are denoted e ± . If the orientation convention is chosen in such a way that c < 0 on the drawing, the densities on

r( ) = r ± c + O( c)
N ± are r ± n = r ± ±e ± c+O ( 2
)
. In this example, monolayer + is constituted of two different lipid species.

and the neutral surface of the monolayer. On this surface, the density which minimizes f for any given membrane shape is ρ 0 (at first order in ).

Let us examine the case of a plane monolayer with fixed total mass. Its equilibrium density ρ eq = ρ 0 (1 + r eq ) can be obtained by minimizing its free energy per unit mass

f ρ = σ 0 ρ 0 (1 -r) + k + 2σ 0 2ρ 0 r 2 + O ( 3 ) (28)
with respect to r. We obtain

r eq = σ 0 k + 2σ 0 , (29)
and thus

σ 0 (ρ 0 ) = k ρ eq -ρ 0 ρ eq + O ( ρ eq -ρ 0 ρ eq ) 2 . ( 30
)
In particular, σ 0 vanishes when ρ 0 = ρ eq . Note that we have assumed that our reference density ρ 0 was sufficiently close to ρ eq for the second-order expansion of f to be valid at ρ eq .

Two-component monolayer

Let us now consider the case of a two-component monolayer. The mass fraction φ of one lipid species must be taken into accout in our monolayer model. In order to study small variations of φ around a reference value φ 0 , we introduce a fourth small variable

ψ = φ -φ 0 φ 0 = O( ) . ( 31
)
The expansion of f can now be written as:

f = σ 0 + σ 1 ψ + σ 2 2 ψ 2 + k 2 (r -ec) 2 - κ 2 (c 0 + c0 ψ) c + κ 4 c 2 + κ 2 c 1 c 2 + O ( 3 ) , ( 32
)
where all the coefficients depend on φ 0 as well as on ρ 0 . It is not necessary to include a term in rψ in this expansion. As a matter of fact, the conservation of the total mass of each lipid species in the monolayer entails that the integral of ρφ is a constant as well as the one of ρ, so a linear term in rψ would be redundant with the one in ψ.

The equilibrium density ρ eq for a flat monolayer with a fixed total mass and a uniform lipid composition φ = φ 0 has the same expression (29) as in the case of a singlecomponent monolayer, but its value depends on φ 0 since σ 0 and k do.

Expanding f in terms of small variables relies on the assumption that f is analytical. However, in the case where φ is very small, the free energy contains a non-analytic part, which reads per unit surface σ(1 + r)φ ln φ, with σ = k B T ρ 0 /µ 1 , where µ 1 is the mass of one lipid of the species with mass fraction φ (see, e.g., the entropy of mixing in Ref. [START_REF] Doi | Introduction to polymer physics[END_REF]). In order to include this term, let us write f as:

f = σ 0 + σ 1 ψ + σ 2 2 ψ 2 + σ (1 + r) φ ln φ + k 2 (r -ec) 2 - κ 2 (c 0 + c0 ψ) c + κ 4 c 2 + κ 2 c 1 c 2 + O ( 3 ) , ( 33
)
where the term between double square brackets must be taken into account only if φ is very small, in which case ψ stands for φ and not for (φ -φ 0 )/φ 0 . In the following, the double square brackets will always be used with this meaning.

Consistency with the ADE model

The area-difference elasticity (ADE) model [17][18][START_REF] Seifert | The Structure and Conformation of Amphiphilic Membranes[END_REF] can be deduced from this model by considering a membrane made of two one-component monolayers (denoted + and -as on figure 1) with fixed total masses M ± , and by eliminating the densities by minimization [18]. Choosing ρ ± 0 = ρ ± eq , which implies σ ± 0 = 0, we have

f ± = k ± 2 ( r ± ± e ± c ) 2 ± κ ± c ± 0 2 c + κ ± 4 c 2 + κ± 2 c 1 c 2 . ( 34
)
Since c 1 and c 2 are defined on the same surface S of the bilayer, and since a single orientation convention is adopted for both monolayers, the curvatures are common to the two monolayers. Here e + (resp. e -) is the (positive) distance between S and the neutral surface of bilayer + (resp.

-). The ± sign in front of c ± 0 ensures that two identical monolayers forming a bilayer would share the same value of c 0 . Note that the free energy density f defined in the previous section corresponds to f -.

Let us minimize with respect to r ± the free energy of each monolayer with fixed total mass M ± and total area A (defined on S, like r ± and c 1 and c 2 ). To take into account the constraints associated with the total masses, we introduce two Lagrange multipliers λ ± . We thus minimize f ± -λ ± ρ ± 0 (1 + r ± ) with respect to r ± , which gives

r ± = λ ± ρ ± 0 k ± ∓ e ± c . ( 35
)
Using the constraint

∫ A dA ρ ± 0 (1 + r ± ) = M ± yields λ ± ρ ± 0 k ± = M ± ρ ± 0 A -1 ± e ± A ∫ A dA c . ( 36
)
Let us define the relaxed area A ± 0 of a monolayer as the area (defined on S) that it would spontaneously adopt in the absence of any constraint:

A ± 0 = M ± /ρ ± 0 , ( 37
)
and its actual area A ± measured on the neutral surface of the monolayer, which verifies:

A ± A = 1 ∓ e ± A ∫ A dA c + O ( 2 ) . ( 38
)
We may then rewrite [START_REF] Doi | Introduction to polymer physics[END_REF] as

λ ± ρ ± 0 /k ± = (A ± 0 -A ± )/A, so that (35) becomes r ± ± e ± c = A ± 0 -A ± A . ( 39
)
Thus, after this partial minimization with respect to r ± , the monolayer free-energy density can be written as

f ± = k ± 2 ( A ± -A ± 0 A ) 2 ± κ ± c ± 0 2 c + κ ± 4 c 2 + κ± 2 c 1 c 2 .(40)
The total free energy F b of the bilayer is obtained by integrating f + + f -over S. Defining the bilayer elastic constants and spontaneous curvature as

κ b = κ + + κ - 2 , κb = κ+ + κ- 2 , c b 0 = κ -c - 0 -κ + c + 0 κ + + κ -, (41)
we obtain

F b = F b A + κ b 2 ∫ A dA c 2 -κ b c b 0 ∫ A dA c + κb ∫ A dA c 1 c 2 , ( 42
)
where

F b A = [k + (A + -A + 0 ) 2 + k -(A --A - 0 ) 2 ]/(2A).
If the position of S in the bilayer, which was arbitrary until now, is chosen in such a way that

e -k -= e + k + , ( 43
)
F b A can be written as

F b A = k b A (A -A 0 ) 2 + K b 4A (∆A -∆A 0 ) 2 , ( 44
)
with

∆A = A --A + , ∆A 0 = A - 0 -A + 0 , and
k b = k + + k - 2 , K b = 2k + k - k + + k -, A 0 = k + A + 0 + k -A - 0 k + + k - .
(45) Condition (43) expresses that S is the neutral surface of the bilayer [20]. Thus, e + (resp. e -) is the distance between the neutral surface of the monolayer + (resp. -) and the neutral surface of the bilayer.

The total free energy F b of the bilayer ( 42)-( 44), obtained from our monolayer free-energy density, corresponds to the ADE model [17,18]. Note that one usually replaces the term ∝ (A -A 0 ) 2 by a hard constraint A = A 0 [START_REF] Seifert | The Structure and Conformation of Amphiphilic Membranes[END_REF], while it is important to keep the area-difference elasticity term ∝ (∆A -∆A 0 ) 2 .

Projected stress tensor 4.1 Stress tensor of a monolayer

In order to obtain the projected stress tensor of our monolayer, we need to determine the projected energy density f (r, ψ, h i , h ij ) associated with f (r, ψ, c 1 , c 2 ). We now assume that the monolayer exhibits only weak deviations from the plane (x, y). Then, if h i = O( ) and h ij = O( ), Eqs. ( 23)-( 24) are satisfied because 4). The projected free-energy density f and the projected mass density ρ are given by

c 1 + c 2 = c = ∇ 2 h + O ( 3 ) and c 1 c 2 = det(h ij ) + O(
f = f √ 1 + (∇h) 2 = f + σ 0 2 (∇h) 2 + O ( 3 ) , ( 46
) ρ = ρ √ 1 + (∇h) 2 = ρ + ρ 0 2 (∇h) 2 + O ( 3 ) , ( 47
)
Hence, defining the scaled projected density as

r = ρ -ρ 0 ρ 0 = r + 1 2 (∇h) 2 + O ( 3 ) , ( 48
)
we obtain from Eq. ( 33)

f = σ 0 + σ 1 ψ + σ 2 2 ψ 2 + σ 0 2 (∇h) 2 + k 2 ( r -e∇ 2 h ) 2 + κ 4 ( ∇ 2 h ) 2 - κ 2 (c 0 + c0 ψ) ∇ 2 h + κ 2 det(h ij ) + σ (1 + r) φ ln φ + O ( 3 ) . ( 49
)
We can now calculate the components of the projected stress tensor from ( 15) and ( 16), noting that ρ ∂ f /∂ ρ = (1 + r) ∂ f /∂ r. We will restrict ourselves to the first order in here, but the tangential components of the stress tensor are calculated at second order in Appendix A. We obtain

Σ ij = [ σ 0 + σ 1 ψ -k ( r -e∇ 2 h ) - κc 0 2 ∇ 2 h ] δ ij + κc 0 2 h ij + O ( 2 ) , ( 50
)
Σ zj = σ 0 h j + ke ∂ j ( r -e∇ 2 h ) - κ 2 ∂ j ∇ 2 h + κc 0 2 ∂ j ψ + O ( 2 ) . ( 51
)
Note that the non-analytic term coming from the entropy of a mixture, that we have put between double square brackets in (49), does not contribute to the stress tensor. This is because this term is proportional to ρ, so it disappears when one computes f -ρ ∂ f /∂ ρ. The expression of the stress tensor is thus the same whether φ is small or whether it is close to a finite value φ 0 . For h = 0, r = 0 and ψ = 0, (50) gives Σ ij = σ 0 δ ij . Hence σ 0 can be interpreted as the tension of a flat membrane with uniform density ρ 0 and uniform lipid composition φ 0 . It is consistent with the fact that σ 0 vanishes for ρ 0 = ρ eq . Besides, (30) may now be interpreted as a Hookean law for the tension of a flat membrane with no inhomogeneities [39].

Comparison with the stress tensor associated with the Helfrich model

Let us compare our results with those coming from the Helfrich model, in which the tension σ is a phenomenological parameter, namely the Lagrange multiplier implementing the area constraint. For a monolayer with elastic constants 1 2 κ and 1 2 κ, tension σ and spontaneous curvature c 0 , the Helfrich free energy density is f = σ + κc 2 /4κc 0 c/2 + κc 1 c 2 /2, so that its projected version reads at second order in :

f = σ + σ 2 (∇h) 2 + κ 4 ( ∇ 2 h ) 2 - κc 0 2 ∇ 2 h+ κ 2 det(h ij ). ( 52
)
The corresponding stress tensor takes the form [32]:

Σ H ij = ( σ - κc 0 2 ∇ 2 h ) δ ij + κc 0 2 h ij + O ( 2 ) , ( 53
)
Σ H zj = σh j - κ 2 ∂ j ∇ 2 h + O ( 2 ) . ( 54
)
Comparing Eqs. ( 53)-(54) with Eqs. ( 50)-( 51), we find that we may write

Σ ij = Σ H ij + O ( 2 ) , ( 55
)
Σ zj = Σ H zj -e ∂ j σ + ( κc 0 2 + eσ 1 ) ∂ j ψ + O ( 2 ) , (56)
if we define

σ = σ 0 + σ 1 ψ -k ( r -e∇ 2 h ) + O ( 2 ) . ( 57
)
Thus, if the scaled lipid composition ψ and the scaled density on the monolayer neutral surface rn = r -e∇ 2 h + O

( 2 ) are both homogeneous, σ is a constant. The stress tensor then has the same form in our model as in the Helfrich model. But contrary to the Helfrich tension, our σ, which may be viewed as a dynamical surface tension, can feature inhomogenities. In the inhomogeneous case, new terms appear in Σ zj . Thus, our stress tensor extends the one associated with the Helfrich model [32] to the case where the lipid density and composition are not homogeneous.

Stress tensor in the ADE model

For a one-component monolayer, the components of the stress tensor at first order in (50)-( 51) are explicitly given by

Σ xx = σ 0 -k [r -e (h xx + h yy )] - κc 0 2 h yy , ( 58
)
Σ xy = κc 0 2 h xy , ( 59
)
Σ zx = σ 0 h x + ke rx - κ 2 (h xxx + h xyy ) , ( 60
)
where ri ≡ ∂ i r and κ = κ + 2ke 2 . The other three components follow from exchanging x and y.

Even when the membrane exhibits large-scale deformations, it is possible to express the stress tensor at a given point M in the local tangent frame (X, Y ) diagonalizing the curvature tensor. Calling X (resp. Y ) the principal direction associated with the principal curvature c 1 (resp. c 2 ), we have

h X = h Y = h XY = 0, h XX = c 1 and h Y Y = c 2 at point M. Hence r = r and c = ∇ 2 h = c 1 + c 2 .
The components of the projected stress tensor read at first order:

Σ XX = σ 0 -k (r -ec) - κc 0 2 c 2 , ( 61
)
Σ Y Y = σ 0 -k (r -ec) - κc 0 2 c 1 , ( 62
)
Σ XY = Σ Y X = 0 , ( 63
)
Σ ZX = - κ 2 ∂ X c . ( 64
)
The tangential stress tensor is thus diagonal.

Choosing ρ 0 = ρ eq and using (39), which comes from partial minimization of the monolayer free energies with respect to r ± , we obtain at first order for each monolayer

Σ ± XX = k ± A ± -A ± 0 A ± κ ± c ± 0 2 c 2 , ( 65
)
Σ ± Y Y = k ± A ± -A ± 0 A ± κ ± c ± 0 2 c 1 . ( 66
)
Summing the contributions from the two monolayers, we obtain the stress tensor of a bilayer in the ADE model, still at first order:

Σ b XX = 2k b A -A 0 A -κ b c b 0 c 2 , ( 67
)
Σ b Y Y = 2k b A -A 0 A -κ b c b 0 c 1 , ( 68
)
Σ b ZX = -κ b ∂ X c . ( 69
)
where k b , A 0 , κ b and c b 0 are defined in Sec. 3.2. In the Helfrich model, the stress tensor of a bilayer with elastic constant κ b and spontaneous curvature c b 0 can be written in the principal tangent frame from (53) and (54). It reads at first order in :

Σ H XX = σ -κ b c b 0 c 2 , ( 70
)
Σ H Y Y = σ -κ b c b 0 c 1 , ( 71
)
Σ H XY = Σ H Y X = 0 , ( 72
)
Σ H ZX = -κ b ∂ X c . ( 73
)
Thus, the stress tensor in the ADE model has the same form as the one in the Helfrich model, with σ = 2k b (A -A 0 ) /A. In light of the previous section, this equivalence is not surprising since in the ADE model, φ = 0, and rn is homogeneous as shown by (39).

Force density in a monolayer

Calculation from the projected stress tensor

Now that we have obtained the stress tensor for our monolayer model, we can calculate the corresponding force per unit area p by taking the divergence of ( 50) and ( 51): p i = ∂ j Σ ij and p z = ∂ j Σ zj . This force per unit area coming from the rest of the monolayer plays an important part in a dynamical description of a membrane. Indeed, its tangential component is a term of the generalized Navier-Stokes equation describing the monolayer, while its normal component is involved in the normal force balance with the external fluid (see, e.g., Seifert and Langer [25]). We obtain at first order in :

p i = -k ∂ i ( r -e ∇ 2 h ) + σ 1 ∂ i ψ , ( 74
)
p z = σ 0 ∇ 2 h - κ 2 ∇ 4 h + ke ∇ 2 r + κc 0 2 ∇ 2 ψ . ( 75
)
Both of these results give back those of Ref. [25] in the particular case of a bilayer constituted of two identical one-component monolayers, if the reference density is ρ 0 = ρ eq . We have thus justified the expression of these force densities from the membrane stress tensor and generalized them. This force density can also be derived from the general expressions ( 18) and [START_REF] Seifert | The Structure and Conformation of Amphiphilic Membranes[END_REF]. Note that (18) indicates that p z = -δF/δh, which is indeed the force taken into account in Eq. (3) of Ref. [25]. Besides, applying [START_REF] Seifert | The Structure and Conformation of Amphiphilic Membranes[END_REF] to the freeenergy density (49) with φ = 0 shows that, in this specific case,

p i = -∂ i (δF/δr)+O ( 2 )
. This justifies the "gradient of the surface pressure" term -∇(δF/δr) used in Eq. ( 4) of Ref. [25].

Direct covariant calculation

Definitions and notations

In this section, we will not restrict ourselves to membranes undergoing small deformations around the flat shape. In general, a membrane can be considered, in a coarse-grained description, as a two-dimensional surface embedded in the three-dimensional space. The position of a fluid element in the membrane can be described by a three-dimensional vector R(u 1 , u 2 ), where u 1 and u 2 are two parameters labelling each fluid element. Mathematically, these parameters are internal coordinates in the two-dimensional surface, and physically they correspond to Lagrangian coordinates.

We are now going to introduce some basic definitions and notations used to describe the shape of a surface in differential geometry [START_REF] Aris | Vectors, tensors, and the basic equations of fluid dynamics[END_REF]34,35]. At each point R of the surface, it is possible to define two vectors tangent to the surface through

t α = ∂R ∂u α ≡ ∂ α R , ( 76
)
where α ∈ {1, 2}. These two vectors are supposed to be linearly independent. Thus,

n = t 1 × t 2 |t 1 × t 2 | (77)
is a unit normal to the surface at point R. The metric tensor of the surface can be expressed as

a αβ = t α • t β , ( 78
)
so that the area element of the surface reads

dA = √ a d 2 u , ( 79
)
where a is the determinant of a αβ , and d 2 u = du 1 du 2 . The inverse metric tensor a αβ is defined by the relation

a αβ a βγ = δ α γ , ( 80
)
where δ α γ is the Kronecker symbol. In the last relation, as well as in the following, the Einstein summation convention is used. We may now define the contravariant tangent vectors as

t α = a αβ t β . ( 81
)
A complete description of a surface is given by its metric tensor (or first fundamental form) and its curvature tensor (or second fundamental form)

b αβ = n • ∂ α t β = n • ∂ α ∂ β R . ( 82
)
The principal curvatures c 1 and c 2 of the surface are the eigenvalues of b α β = a αγ b γβ , which enables to express the total curvature and the Gaussian curvature from the curvature tensor:

c = c 1 + c 2 = b α α , (83) c 1 c 2 = det b α β . (84)

Force density in a monolayer

The surface density of internal forces q in a two-component monolayer with free energy F = ∫ dA f can be expressed as the functional derivative

q(u 1 , u 2 ) = - 1 √ a δF δR(u 1 , u 2 ) ρ √ a, ρφ √ a , ( 85
)
where, as in the previous sections, ρ is the total mass density of lipids, and φ the mass fraction of one lipid species [34,35,[START_REF] Lomholt | [END_REF]. This expression is a consequence of the principle of virtual work: for a small deformation δR of the membrane at equilibrium, the membrane free-energy variation reads

δF = - ∫ dA q • δR = - ∫ d 2 u √ a q • δR . ( 86
)
For the underlying force balance on each fluid element to be valid, the virtual deformation δR must be performed at constant total mass dm = ρ dA = ρ √ a d 2 u and composition in each fluid element. Hence, the functional derivative in (85) must be taken at constant ρ √ a and ρφ √ a. Let us calculate the force density (85) in a monolayer with free-energy density

f = σ 0 + σ 1 ψ + σ 2 2 ψ 2 + σ (1 + r) φ ln φ + k 2 (r -ec) 2 - κ 2 (c 0 + c0 ψ) c + κ 4 c 2 + κ 2 c 1 c 2 . ( 87
)
This free-energy density corresponds to (33) truncated at second order in . The two constraints on the deformation δR, δ(ρ √ a) = 0 and δ(ρφ √ a) = 0, are equivalent to

√ a δr + (1 + r)δ √ a = 0 (88) √ a δψ = 0 , ( 89
)
where we have used the fact that 1 + r > 0. To enforce these two independent local constraints, we use two local Lagrange multipliers, λ and µ. The principle of virtual work (86) then reads

δF - ∫ d 2 u { λ [√ a δr + (1 + r)δ √ a ] + µ √ a δψ } = δW , ( 90
) where δW = - ∫ dA q • δR . ( 91
)
We assume that the topology of the membrane is not affected by the virtual deformation. The Gauss-Bonnet theorem then ensures that δ (∫ dA c 1 c 2 ) = 0. The transformation of the left-hand side of (90) can be performed along the same lines as in Ref. [42]. These calculations, which are presented in Appendix B, yield:

q • t α = -k (1 + r) ∂ α (r -ec) + ( σ 1 + σ 2 ψ - κc 0 2 c ) ∂ α ψ , ( 92
)
q • n = ( σ 0 + σ 1 ψ + σ 2 2 ψ 2 - k 2 r 2 -kr ) c - κ 4 c 3 + [κc -κ (c 0 + c0 ψ) -2ke r] c 1 c 2 +ke (1 + r) c 2 + ke ∆r - κ 2 ∆c + κc 0 2 ∆ψ , ( 93
)
where ∆ is a shorthand for the Laplace-Beltrami operator (1/ √ a)∂ α (a αβ √ a∂ β ). The force density q can be expressed from its tangential component (92) and normal component (93) as

q = (q • t α )t α + (q • n)n . ( 94
)
We have thus obtained the general expression of the force density in a two-component monolayer with free-energy density (87). This expression gives back the one in Ref. [35] in the particular case of a one-component monolayer with c 0 = 0. Note that, in this section, we have used the free-energy density truncated at second order (87) as if it were exact. The force density q expressed in (92)-( 94) is the one corresponding to this model, and it contains second and thirdorder terms. This approach is consistent with the one of Refs. [42,34,35]. However, in the present paper, we have constructed the free-energy density as a general expansion around a reference state, controlling the order in of this expansion. In our approach, if the free-energy density f is kept at second order, the force density can be known only at first order.

Comparison between the two results

In the present paper, except in the previous section 5.2, we have described membranes in the Monge gauge, i.e., by their height with respect to a reference plane. Such a description is very convenient to study the membrane small deformations around the flat shape. In the Monge gauge, the position of a fluid element in the membrane is given by R(x, y) = (x, y, h(x, y)), where x and y are Cartesian coordinates in the reference plane and z = h(x, y) describes the height of the membrane with respect to the reference plane. Obviously, the coordinates (x, y) of a fluid element depend on its position in the membrane (in other words, they are Eulerian coordinates).

In Sec. 5.2, we have found the force density q in a monolayer, whatever its shape. In our derivation of q, the parameters (u 1 , u 2 ) describing the surface were Lagrangian coordinates, labelling each fluid element. However, as the force density is a physical quantity, it does not depend on the parametrization of the surface [34]. Thus, the expression we have found for q is valid (for each given membrane shape) in the Monge gauge.

We may now compare the force density obtained in Sec. 5.2 with the one obtained from the projected stress tensor in Sec. 5.1. For this, we shall write explicitly in the Monge gauge the general result obtained in Sec. 5.2. With u 1 = x and u 2 = y, the tangent vectors read in the Monge gauge t 1 = (1, 0, ∂ x h) and t 2 = (0, 1, ∂ y h). It is then straightforward to find the expression of n, a αβ and a αβ in the Monge gauge (see, e.g., Ref. [34]). Using these explicit expressions, and keeping only first order terms in , ( 92) and (93) can be written as:

q i = -k ∂ i ( r -e ∇ 2 h ) + σ 1 ∂ i ψ + O ( 2 ) , ( 95
)
q z = σ 0 ∇ 2 h - κ 2 ∇ 4 h + ke ∇ 2 r + κc 0 2 ∇ 2 ψ + O ( 2 ) , (96)
where i ∈ {x, y}. We notice that, at this order, (95) is identical to (74) and ( 96) is identical to (75). Note that q is a force density per actual unit area of the monolayer while p is a force density per projected unit area. However, this difference is irrelevant at first order. We have just shown that the force density obtained from the divergence of the projected stress tensor is consistent with the one calculated directly by using the principle of virtual work in covariant formalism. The projected stress tensor thus allows to calculate easily both the normal and the tangential components of the force density in a membrane in the Monge gauge, without having to resort to a covariant formulation.

Applications

The force density in a membrane with lipid density and composition inhomogeneities can be used to understand qualitatively and quantitatively the dynamics of a membrane submitted to a local perturbation. In order to illustrate this, we are going to focus on the local injection of a reagent close to a membrane (see Fig. 2), which modifies locally the properties of the membrane.

Forces arising from a modification of composition

Let us consider initially a one-component flat membrane at equilibrium with uniform lipid density. Let us assume that, at time t = 0, some lipids in the external monolayer of this membrane (monolayer +, as on Fig. 1), are suddenly chemically modified due to a microinjection of a reagent close to the membrane (see Fig. 2). Then there is a fraction φ(x, y), assumed to be small, of modified lipids in this monolayer. The force density in the external monolayer at time t = 0 + , just after the injection, when the shape and the density have not changed yet, is given by:

p + = σ 1 ∇φ + κc 0 2 ∇ 2 φ e z , ( 97
)
where e z is a unit vector in the z direction. This force density corresponds to Eqs. ( 74) and (75) in the case of a flat "+" monolayer with uniform density. Hence, modifying locally the lipids of a monolayer will generically induce a shape instability of the membrane. The internal monolayer is not affected by the chemical modification, so the force density remains zero in it.

Let us take the position of the micropipette injecting the reagent as the origin of our (x, y) frame. Then, φ is a decreasing function of the radial coordinate r. Let us study the case where φ is a Gaussian:

φ(r) = φ 0 exp ( - r 2 2 R 2 ) . ( 98
)
This can represent the field of modified lipids resulting from a diffusion of the reagent in the solution surrounding the vesicle before it hits the membrane. It is straightforward to calculate the corresponding force density. Its normal and radial components, nondimensionalized by their maximal values, are plotted as a function of r/R in Fig. 3.

The constants σ 1 and c0 that appear in the force density arise from the φ-dependence of the free energy per unit area of a monolayer (see Eq. 33). Physically, they describe the change in the equilibrium density and the spontaneous curvature of the membrane due to a generic modification of the lipids [28]. Their signs, and thus, those of the force density components, depend on the nature of the modification. In the case where σ 1 < 0 and c0 < 0, which corresponds to modified lipids favoring a smaller density and a larger curvature (in absolute value), the lipids in front of the pipette are submitted to a normal force going towards the exterior of the vesicle, which will yield a local deformation of the membrane in this direction. Meanwhile, a radial force drives the lipids of the external monolayer to flow in the membrane towards larger values of r, due to the fact that the modified lipids favor a smaller density.

This situation describes well the onset of the shape instability studied in Ref. [28], which is induced by microinjecting a basic solution close to a giant unilamellar vesicle. This instability can yield the formation of a membrane tubule [22].

Forces arising from a local deformation at uniform density

Besides the composition, another important effect captured by our study is the coupling between the membrane shape and the density. To shed light onto this effect, let us consider a locally deformed membrane with uniform density (on the bilayer midsurface). This can correspond to a membrane which has deformed very rapidly from a flat shape, before the density adjusts to the new deformed shape. Indeed, the symmetric density (i.e., the sum of the densities in the two monolayers) is not coupled to the deformation, while the antisymmetric density is [25,24,22,28]. Thus, intermonolayer friction is involved when the density adjusts to the deformation. The associated timescale is a few seconds, so there is indeed a lapse when a deformed membrane with non-adjusted density exists in the experiments described in Ref. [28].

To isolate the effect of the shape and density, we will focus on a one-component membrane here (note that the external monolayer is a two-component one in Ref. [28]). Let us take ρ 0 = ρ eq as our reference density in both (identical) monolayers. Then, Eq. ( 30) ensures that σ 0 = 0. The force density in monolayers "±" caused by the deforma-tion is given by:

p ± = ∓ke ∇(∇ 2 h) - κ 2 ∇ 4 h e z . ( 99
)
This force density corresponds to Eqs. ( 74) and (75) in the case of one-component monolayers with uniform density. Let us consider for instance a Gaussian-shaped deformation towards the exterior, centered on the origin:

h(r) = h 0 exp ( - r 2 2 R 2 ) . ( 100
)
The normal and radial components of the corresponding force density, nondimensionalized by the absolute value of their maxima, are plotted as a function of r/R in Fig. 4. Since ke > 0 and κ > 0, at small r, the lipids are submitted to a normal force going towards the interior of the vesicle: this will lead to a relaxation of the deformation. This is due to the fact that the membrane considered here is symmetric, so its equilibrium shape is flat. Meanwhile, a radial force drives the lipids of the external monolayer to flow in the membrane, for the density to adjust to the shape (see Fig. 5). Indeed, in a curved membrane at equilibrium, the density is uniform on the neutral surface of each monolayer, and not on the membrane midsurface. The radial forces are opposite in the external and in the internal monolayer, because the orientations of these monolayers are opposite while they share the same curvature. Hence, for the density to adjust to the deformed shape, it is necessary that the lipids in one monolayer slide with respect to the ones in the other monolayer, which involves intermonolayer friction [28].

Dynamics

The force densities in the membrane are the basis of a hydrodynamic description of the membrane (for a review, see Ref. [43]). Our work enables to take into account lipid density and composition instabilities in such dynamical Fig. 5. Membrane with Gaussian h, that has just deformed from a flat shape. The density has not adjusted to the new deformed shape yet. The lipids are thus at equal distance on the midlayer (plain line). The equilibrium density for this deformed shape would correspond to lipids at equal distance on the neutral surface of each monolayer (dashed lines). The arrows indicate the direction of the tangential force density in the membrane, consistent with Fig. 4.

studies. In Ref. [28], we have used the force densities derived in the present work to describe the dynamics of a membrane curvature instability induced by a local injection of a basic solution close to a giant unilamellar vesicle.

More precisely, we have written down generalized Stokes equations describing the balance of the forces per unit area acting tangentially in each monolayer. These equations include the tangential density of forces Eq. ( 74), the viscous force density due to the two-dimensional flow of the lipids, the tangential viscous stress exerted by the flow of the surrounding fluid, and the intermonolayer friction. We have also written down the balance of the forces per unit area acting normally to the membrane, including the normal force density given by Eq. (75), and the normal viscous stress exerted by the flow of the surrounding fluid. Our last fundamental dynamical equation is the conservation of mass.

Thanks to these equations, we have described theoretically the evolution of the deformation resulting from the local chemical modification of the external monolayer by the basic solution. This description is developed in detail, and compared to experimental results in Ref. [28].

Conclusion

We have derived a general formula expressing the projected stress tensor in a monolayer as a function of the monolayer free-energy density, taking into account inhomogeneities in the lipid density and composition. This general formula has been applied to a generic monolayer model constructed from basic principles. Our model being consistent with the ADE model, we have found in particular the stress tensor associated with the ADE model.

In the Monge gauge, the projected stress tensor provides a convenient way of deriving the force density in a monolayer, which is the basis of a hydrodynamic description of a membrane. The result is consistent with a direct calculation of the force density from the principle of virtual work in covariant formalism. We have shown an example of application to the calculation of force den-sity in a locally perturbed membrane. These force density constitute the basis of the dynamical study in Ref. [28].

Moreover, the stress tensor contains more information than the force density, since it provides the actual force exerted by a piece of membrane along its edge. Indeed, the stress tensor associated with the Helfrich model [31,32] has already been used to study various situations, such as the boundary conditions on a membrane with a free exposed edge [44], the adhesion of a fluid membrane [45,46], the surface tension of fluctuating membranes [47,48] and the force exerted by a fluctuating membrane tubule [49]. The stress tensor is also a useful tool in the study of membrane-mediated interactions [50,51,30]. Since the stress tensor studied in the present paper generalizes the one associated with the Helfrich model to monolayers with density and composition inhomogeneities, it may enable to extend such applications.

A Tangential stress tensor at second order

For the Helfrich model, the tangential components of the stress tensor can be obtained at second order from the free energy at second order [32]. On the contrary, here, because of the term ∂ f /∂ r, one cannot obtain Σ ij at order 2 without taking into account in f the terms of order 3 which depend on r. When such third order terms are included, f becomes

f = σ 0 + σ 1 ψ + σ 2 2 ψ 2 + σ (1 + r) φ ln φ + k 2 (r -ec) 2 - κ 2 (c 0 + c0 ψ) c + κ 4 c 2 + κ 2 c 1 c 2 + σ 1 2 r 2 ψ + σ 2 2 rψ 2 + k 3 r 3 + κc 0 2 r 2 c - κc 0 2 ψrc + κ 4 rc 2 + κ 2 rc 1 c 2 + O ( 4 ) , (101)
where O ( 4 ) stands for terms of order 4 , or terms of order 3 independent of r. Indeed, the latter may be discarded because they will not contribute to the stress tensor at order 2 . The new terms in f may be considered as originating from a density-dependence of the constitutive constants k, κ, κ, c 0 , c0 , σ 1 and σ 2 .

Using the relations Eqs. ( 15) and (106) yield

f = f √ 1 + (∇h) 2 = f + σ 0 2 (∇h) 2 + O ( 4 ) , ( 102
) c = ∇ 2 h + O ( 4 ) , ( 103
)
c 1 c 2 = det(h ij ) + O ( 4 ) , (104)
Σ ij = { σ 0 + σ 1 ψ + σ 2 -σ 2 2 ψ 2 + σ 0 + k 2 (∇h) 2 -(k + σ 1 ψ) r -(k + 2k ) r2 2 + κ -κ 4 ( ∇ 2 h ) 2 + [ ke - κ 2 (c 0 + c0 ψ) + κc 0 r + κc 0 2 ψ ] ∇ 2 h - κ 2 det (h ij ) } δ ij -σ 0 h i h j + κ 2 (c 0 + c0 ψ) h ij - κ 2
( h ij ∇ 2 h -h i ∂ j ∇ 2 h ) + ke (h ij r -h i ∂ j r) - κc 0 2 h i ∂ j ψ + O ( 3 ) . ( 107
)
where κ = κ + 2ke 2 as before. The non-analytic term present in f at small φ does not contribute to the stress tensor, for the same reason as before.

In the principal tangent frame, the tangential components of the stress tensor at second order are given by

Σ XX = σ 0 + σ 1 ψ + σ 2 -
Σ XY = Σ Y X = 0 . ( 109
)
Σ Y Y can be obtained by exchanging c 1 and c 2 in Σ XX . This tangential stress tensor thus remains diagonal at second order in the principal tangent frame, like the one associated with the Helfrich model [32].

B Covariant calculation of the force density

In this section, we are going to present the main steps of our covariant calculation of the force density in a monolayer, which leads to (92)-(93). This calculation follows the same lines as the one in Ref. [42].

As the Gaussian curvature term contained in the total monolayer free energy

F = ∫ dA f = ∫ d 2 u √ a f ( 110
)
does not vary during the virtual deformation δR, we may replace the free-energy density (87) by

f = f - κ 2 c 1 c 2 (111)
in our calculations. The free-energy variation during the deformation then reads

δF = ∫ d 2 u f δ √ a + ∫ d 2 u √ a δ f , ( 112
)
with

δ f = [ σ 1 + σ 2 ψ - κc 0 2 c + σ (1 + r) (1 + ln φ)
] δψ

+ [ k(r -ec) + σφ ln φ ] δr + I 1 δc , ( 113
)
where we have defined

I 1 = κ 2 (c -c 0 -c0 ψ) -ek (r -ec) . ( 114
)
The principle of virtual work (86) may now be written as

δW = ∫ dA { I 1 δc + [ k(r -ec) + σφ ln φ -λ ] δr + [ σ 1 + σ 2 ψ - κc 0 2 c + σ (1 + r) (1 + ln φ) -µ ] δψ } + ∫ d 2 u [ f -λ (1 + r) ] δ √ a , ( 115
)
where δW is given by (91). The variations δc and δ √ a only come from the variation δR of the shape of the monolayer. The coupling between δR and δr and δψ, which comes from the constraints (88)-(89), has been accounted for by introducing the Lagrange multipliers λ and µ. Hence, δR, δr and δψ should now be considered as independent variations, and the terms in δr and δψ must vanish for (115) to be valid for any virtual deformation, yielding λ = k(r -ec) + σφ ln φ , (116)

µ = σ 1 + σ 2 ψ - κc 0 2 c + σ (1 + r) (1 + ln φ) . ( 117
)
We thus obtain

δW = ∫ d 2 u (√ a I 1 δc + I 2 δ √ a ) , ( 118
)
where we have defined

I 2 = σ 0 + σ 1 ψ + σ 2 2 ψ 2 - κ 2 (c 0 + c0 ψ) c + κ 4 c 2 + k 2 (ec -r) (r + ec + 2) . ( 119
)
Note that the contribution of the term between double square brackets has vanished, as in our calculation of Σ.

Thanks to the relations

δ √ a = √ a t α • δt α , ( 120
) δc = a αβ (∂ α n) • δt β -t α • δ (∂ α n) , ( 121
) δn = -(n • δt β ) t β , ( 122
)
δW can be expressed only in terms of δt α . Performing two integrations by parts and using the relations [START_REF] Aris | Vectors, tensors, and the basic equations of fluid dynamics[END_REF] t α|β = b αβ n , (123)

n |α = ∂ α n = -b αβ t β , ( 124
)
where g |α denotes the covariant derivative (associated with the metric a αβ ) with respect to u α of a function g defined on the surface [START_REF] Aris | Vectors, tensors, and the basic equations of fluid dynamics[END_REF], we obtain

δW = - ∫ dA [ ( a αβ I 1 -b αβ I 2 ) t α + a αβ (∂ α I 2 ) n ] |β • δR . ( 125
)
Identifying ( 125) with (91) for any infinitesimal virtual deformation δR, we obtain the sought force density:

q = [( a αβ I 1 -b αβ I 2 ) t α + a αβ (∂ α I 2 ) n ] |β . ( 126
)
Performing the covariant derivative with respect to u β in (126), using ( 123)-( 124) and the relations [START_REF] Aris | Vectors, tensors, and the basic equations of fluid dynamics[END_REF] b α β|α = ∂ β c , (127)

b αβ b βα = c 2 -c 1 c 2 , ( 128
)
and taking the scalar product of q with t α (resp. n) finally leads to (92) (resp. (93)).

Fig. 2 .

 2 Fig. 2. Local injection of a reagent from a micropipette close to a vesicle. The lipids in the external monolayer of the vesicle will be affected. The injection is sufficiently local for us to focus on a small, nearly-plane zone of the membrane.

Fig. 3 .

 3 Fig. 3. Nondimensionalized force density in the external monolayer of a membrane just after a local chemical modification, with Gaussian φ. Plain line: normal component pz/p max z . Dashed line: radial component pr/p max r .

Fig. 4 .

 4 Fig. 4. Nondimensionalized force density in monolayer "+" just after a local deformation, with Gaussian h. Plain line: normal component pz/|p max z |. Dashed line: radial component pr/|p max r |. In monolayer "-", pz is identical and pr is opposite.

Anne-Florence Bitbol et al.: Membrane stress tensor in the presence of lipid density and composition inhomogeneities


Yasmina Abdeddaïm
email: y.abdeddaim@esiee.fr

Damien Masson
email: d.masson@esiee.fr

A Generic Task Model Using Timed Automata to Design and Analyze Real-Time Applications

Keywords: Timed Automata, Real-Time Scheduling, Self-Suspending Tasks

In this work, we propose an approach for automated analysis of real-time scheduling problems based on timed automata. Tasks are modeled using timed automata while schedulers correspond to feasible runs computed using model checking algorithms. Using this model, we present a method to prove on one hand the feasibility of a scheduling problem and to compute, on the other hand, an hybrid off-line/on-line scheduling policy. We expose how to use it to re-demonstrate well known results, to produce efficient schedule regarding of various properties such as minimal preemption number, and to address the problem of scheduling real-time task set composed by possibly self-suspending tasks which is known to be NP-hard.

Introduction

In real-time systems, the scheduling problem consists in producing a feasible schedule for a real-time task set whenever one exists. A slightly different problem is, given a task set and a scheduling policy, to answer the question: did all the tasks will respect all their deadlines? This is the feasibility problem. The aim of this paper is to answer to these two problems: to provide a framework to both analyze the properties of a given schedule and to generate a feasible schedule (optionally with some parameterizable properties).

The scheduling theory often concentrates on the search for an optimal scheduler. An optimal scheduling policy is a policy that, given a task set, produces a non-feasible schedule if and only if it does not exist any other scheduling policy in the same class of algorithm that can produce a feasible schedule. The precision in the same class of algorithm is of importance. Indeed, scheduling algorithms can be characterized in many ways: off-line or on-line (priority driven), workconserving or idle, preemptive or not. To illustrate the ambiguity of the word optimal, we can cite the RM, DM and EDF algorithms [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real time environment[END_REF], which are all three said optimal. In fact, RM is optimal among fixed priority driven preemptive schedulers, and only for traffic composed by periodic tasks with implicit deadlines (i.e. D i = T i ). DM is optimal among the same kind of schedulers, but for traffic composed by periodic tasks with constrained deadlines (i.e. D i ≤ T i ). Finally EDF is optimal among (not necessarily fixed) priority driven schedulers for traffic composed by periodic tasks with constrained deadlines. Any optimalityrelated result on a scheduling policy (or comparison between several algorithms) so should be considered in a well defined referential, composed by the kind of scheduler considered for comparison, and the kind of scheduling problem considered (hypotheses on the traffic model).

Our approach is to clearly separate the behavior of a real-time task from specific constraints of a scheduling policy. Using a timed automaton [START_REF] Alur | Automata for modeling real-time systems[END_REF], we model a task as a transition system between possible configurations of a task (activated, preempted, running, blocked). Feasible schedules then correspond to specific runs on the model. Different scheduling and feasibility problems can be solved using CTL [START_REF] Kozen | Logics of Programs[END_REF] model checking on the generic model without constraining the model. While changing the property to verify, we can consider feasible runs only the ones which follow specific rules (e.g. the ones leading to an EDF schedule) or just ask if there exists a way to respect real-time constraints without restriction of the scheduling policy. This design permits to use the model either to test properties of a task set (e.g. schedulable, unschedulable, ...), of a specific scheduling policy on a task set (e.g. schedulable with RM, ...) or even to generate a feasible schedule.

One important characterization of a scheduling policy is the separation between off-line schedulers in one side, and on-line schedulers in the other side. A scheduling policy is said on-line if the scheduling decisions (task start, preemption, resume...) are taken at run time. It is said off-line if the actions at runtime are limited to read a table where all scheduling decisions are store. On-line scheduling is mostly represented by a sub-category: the priority-driven schedulers. Such a scheduler maintains up-to-date a sorted-by-priority queue of active tasks, and the processor is allocated to the task in head of this queue. A Distinction can also be made between fixed priority schedulers for which a priority is assigned to each task once and for all, and dynamic priority schedulers for which the priority of a task can change at any time. Off-line scheduling policies are easy to implement and have a very low overhead cost. However, they suffer from a possibly high space complexity issue, and are not suitable on systems that have to react dynamically to external and/or unpredictable behavior such as aperiodic event happening, cost under-run or over-run of a task or deadline missed. Moreover, the generation of the table can have a high time complexity, depending of the scheduling problem considered.

In the case of scheduling problem without uncertainties, the timed automaton model permits to produce off-line a scheduling table that can be read online. In scheduling under temporal uncertainties, the model produces a kind of hybrid scheduling policy. Decisions are taken on-line according to data computed off-line. It is suitable to schedule highly dynamics systems even if the space complexity should be examined in details. Such a scheduling policy can be classified among clairvoyant algorithms. This is an interesting property since a lot of negative results for numbers of scheduling-related problems, such as NP-completeness, are applicable only to non-clairvoyant algorithms.

To illustrate the interest of our approach, we focus on one open problem in real-time systems theory: the scheduling of possibly self-suspending real-time tasks. This problem has recently been shown to be NP-hard in the strong sense [START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF], i.e. there is no optimal on-line algorithm to schedule it which takes its decisions in a polynomial time (unless P = N P ). Using our approach to model self-suspending tasks, we can both determine if a schedule exists, and provide the scheduler with extra data to drive it to an existing schedule. Moreover, the use of the model permits us to conjecture that the schedulability of such task set is sustainable regarding the variations either on execution duration and suspension duration for work-conserving scheduler (this property can be trivially proved true with an idle scheduler). The sustainability is an important property of a scheduling problem solution: it says that if a task set is feasible, then the same task set with positive modifications (such as reduced execution time or augmented period) is still guaranteed feasible. This property permits to answer the feasibility problem on the based of the worst case execution times for tasks with uncertain durations.

After a review of related works in the remainder of this Section, we introduce timed automata in Section 2. Section 3 presents our model. Section 4 exposes how to use it to test systems schedulability. In Section 5 we present the problem of schedule self-suspending tasks and how to handle it with our model. Finally we conclude in Section 6.

Task Model. In this paper, we consider real-time systems built from a set of n periodic real-time tasks τ 1 , τ 2 , ..., τ n . Each task τ i is characterized by a period T i , a worst-case execution time C i and a relative deadline D i .

Related Work. In [START_REF] Abdeddaim | Job-shop scheduling using timed automata[END_REF][START_REF] Fehnker | Scheduling a steel plant with timed automata[END_REF] timed automata are used to solve mainly non preemptive job shop scheduling problem. In these approaches, reachability analysis algorithms are used to construct optimal schedules in the sense of minimal total execution time. Handling preemptible tasks using timed automata is not possible mainly because clock variables cannot be stopped and modeling such tasks can only be done using stopwatch automata. Unfortunately, reachability analysis of stopwatch automata is an undecidable problem in the general case [START_REF] Kesten | Decidable integration graphs[END_REF]. However, there exits methods [START_REF] Waszniowski | Over-approximate model of multitasking application based on timed automata using only one clock[END_REF][START_REF] Gabor | A conservative approximation method for the verification of preemptive scheduling using timed automata[END_REF] by which we can model preemptible tasks in timed automata with over-approximation and modeling with stopwatches is recently possible using the tool Uppaal 4.1 [START_REF] David | Model-Based Framework for Schedulability Analysis Using UPPAAL 4.1[END_REF] with the support of a zone based over-approximation state exploration. In [START_REF] Fersman | Timed automata with asynchronous processes: schedulability and decidability[END_REF], authors proposed a schedulability analysis algorithm where the scheduler is modeled using an extension of timed automata called suspension automata [START_REF] Mcmanis | Suspension automata: A decidable class of hybrid automata[END_REF]. Unfortunately the number of clocks needed is proportional to the maximal number of schedulable task instances associated with the model. The authors then presents an efficient algorithm in [START_REF] Fersman | Schedulability analysis of fixed-priority systems using timed automata[END_REF] for the problem of fixed priority scheduling, as the problem is undecidable when the execution times are given as intervals, they propose an over-approximation technique. The control synthesis approach has been used to solve the job scheduling problem in the case of uncertain execution times [START_REF] Abdeddaim | On optimal scheduling under uncertainty[END_REF], and in [START_REF] Altisen | A methodology for the construction of scheduled systems[END_REF][START_REF] Altisen | Scheduler modeling based on the controller synthesis paradigm[END_REF] control synthesis is used to construct a scheduler by adding control invariant modeling the schedulability constraints.

Preliminaries

Timed Automata A timed automaton [START_REF] Alur | Automata for modeling real-time systems[END_REF] is a finite automaton augmented with a finite set of real-valued variables evolving continuously and synchronously with absolute time. Formally, let X be a set of real variables called clocks and C(X ) the set of clock constraints φ over X generated by the following grammar:

φ ::= x♯c | x -y♯c | φ ∧ φ where c ∈ IN , x,
y ∈ X , and ♯ ∈ {<, ≤, ≥, >}. A clock valuation is a function v : X → R + ∪ {0} which associates with every clock x its value v(x). Given a real d ∈ R we write v + d for the clock valuation associating with clock x the value v(x)

+ d. If r is a subset of X , [r = 0]x is the valuation v ′ such that v ′ (x) = 0 if x ∈ r, and v ′ (x) = v(x) otherwise. Definition 1 (Timed Automaton). A timed automaton (TA) is a tuple A = (Q, q 0 , X , I, ∆, Σ) where, Q is a finite set of states, q 0 is the initial state, X is a finite set of clocks, I : Q → C(X ) if the invariant function, ∆ ⊆ Q × C(X ) × Σ × 2 X × Q is a finite set of

transitions and Σ is an alphabet of actions augmented with the action ⊥ that represents the empty action

A configuration of a timed automaton is a pair (q, v) consisting of a state q and a dim(X ) vector v of clock valuations. The semantic of a timed automaton is given as a timed transition system with two types of transition between configurations defined by the following rules:

a discrete transition (q, v) a -→ (q ′ , v ′ ) where there exists δ = (q, φ, a, r, q ′ ) ∈ ∆ such that v satisfies φ and

v ′ = [r = 0]x, -a timed transition (q, v) d -→ (q, v + d1) d ∈ R + with v and v + d satisfying I(q) the invariant of state q.
A run of the timed automaton starting from a configuration (q 0 , v 0 ) is a sequence of time and discrete transitions ξ : (q 0 , v 0 )

t1 -→ (q 1 , v 1 ) t2 -→ • • • .
The reachability problem for timed automata is decidable and PSPACEcomplete [START_REF] Alur | Automata for modeling real-time systems[END_REF][START_REF] Alur | A theory of timed automata[END_REF]. Since this result, many model checking algorithms [START_REF] Alur | Model-checking in dense real-time[END_REF][START_REF] Henzinger | Symbolic model checking for real-time systems[END_REF][START_REF] Larsen | Model-checking for real-time systems[END_REF] and timed model checkers [START_REF] Yovine | Kronos: A verification tool for real-time systems[END_REF][START_REF] Larsen | Uppaal in a nutshell[END_REF] have been developed.

This basic model of timed automata can be extended to allow the use of integer variables. Let V be a set of integer variables and B(V) a set of integer constraints over V generated by the following grammar:

φ V ::= x < c | x ≤ c | x > c | x ≥ c | φ ∧ φ where c ∈ IN and x ∈ V.
We define a timed automaton using integer variables as a tuple A = (Q, q 0 , X , V, I, ∆, Σ) where a transition (q, φ x , φ v , a, r x , r v , q ′ ) ∈ ∆ is defined, in addition to the basic transitions of a timed automaton, by the integer guard φ v ∈ B(V) and the set r v of linear update functions over integer variables.

A network of timed automata is the parallel composition A 1 || . . . ||A n of a set of timed automata A 1 , . . . , A n . Parallel composition used an interleaving semantic and synchronous communication can be done using input actions denoted a? and output actions denoted a! (this hand shake synchronization is used in the tool Uppaal [START_REF] Larsen | Uppaal in a nutshell[END_REF]) while asynchronous communication is done using shared variables.

Timed Game Automata A timed game automaton model [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF] is an extension of the timed automaton model, which has been introduced for the synthesis of timed controllers.

Definition 2 (Timed Game Automaton). A timed game automaton (TGA) is a timed automaton where the set of transitions ∆ is split into controllable (∆ c ) and uncontrollable (∆ u ) transitions.

The timed game automaton model defines the rules of a two players game with on one side the controller (mastering the controllable transitions) and on the other side the environment (mastering the uncontrollable transitions). Given a timed game automaton and a logic formula, solving a timed game consists in finding a strategy f s.t. the automaton starting at the initial configuration and supervised by f always satisfies the formula whatever actions are chosen by the environment. More precisely, a strategy is a partial mapping f from the set of runs of the TGA to the set ∆ c ∪ {λ} s.t. for a finite run ξ, if f (ξ) = e ∈ ∆ c then execute the controllable transition e from the last configuration of the run ξ and if f (ξ) = λ then wait in the last configuration of the run ξ. A strategy is memory-less whenever its result depends only on the last configuration of the run. It has been shown that solving a timed game is a decidable problem [START_REF] Maler | On the synthesis of discrete controllers for timed systems (an extended abstract)[END_REF][START_REF] De Alfaro | Symbolic algorithms for infinite-state games[END_REF].

3 Real-Time Task Automata Model We present our method by introducing a generic model for scheduling using timed automata. For the sake of clarity, we assume a single processor scheduling s e p problem and all tasks periodic with their deadlines equal to their periods. However, the proposed approach can be easily generalized to less specific scheduling problems.

s e f p cτ = C -1 ∧ xτ = 1 xτ := 0 ∧ t := 1 cτ < C ∧ xτ = 0 ∧ t = 0 xτ := 0 xτ := 0 ∧ cτ := 0 ∧ t := 1 xτ ≤ 1 cτ < C -1 ∧ xτ = 1 cτ := cτ + 1 ∧ xτ := 0 ∧ t := 0
xτ := 0 ∧ t := 1 cτ < 6 ∧ xτ = 0 ∧ t = 0 cτ = 4 ∧ xτ = 1 f xτ = 1 xτ ≤ 2 cτ < 4 ∧ xτ = 1 xτ := 0 ∧ cτ := 0 ∧ t := 1 xτ ≤ 1 cτ := cτ + 2 ∧ xτ := 0 ∧ t := 0
A real-time task can be among one of the possible configurations: active, running, waiting for the next activation, or preempted. The timed automaton model is used in our approach to model a real-time task as a transition system where the set of states represents the possible configurations of a task and the set of transitions the possible moves between different configurations of a task.

For every real-time task τ = (C, D, T ), we associate a 4-states timed automaton A τ with one clock x τ , one integer variable c τ and a set of states Q = {s, e, f, p}. State s is the waiting state where the task is active and not yet executed, e is the active state where the task is running, p is the state where the task is preempted and f is the one where the task is waiting for the next activation. An example is given by Figure 1. To compute the execution time of a task, we use the clock x τ and the variable c τ as follows: the automaton can stay in the execution state exactly one time unit and the variable c τ keeps track of how many time unit has been performed. This is represented by an invariant x τ ≤ 1 on state e and a loop transition that increments c τ . The transition from e to f is enabled when the total time spent in the active state is equal to C (the completion time of the task τ ). The preemption of the task is modeled using transitions from e to p and from p to e. Note it is supposes in this model that a task can be preempted only at integer times, moreover, if we consider that preemptions can occur only at certain known integer points we can arise the staying condition of state e as shown in Figure 2.

To make this task periodic of period T , we use a second timed automaton A T with one state and a loop transition enabled every T time unit. This transition is labeled with an output action T ! and synchronize with the transition from waiting to active state of the automaton A τ . Definition 3 (A real-time task automata model). Let τ (C, T, D) be a realtime periodic task with T = D. A real-time task automata model for τ is a tuple (A τ , A T ) where A τ and A T are two timed automata defined as follows:

1. A τ (Q τ , s, x τ , c τ , I τ , ∆ τ , Σ τ ):
-Q τ = {s, e, f, p}, -I τ (e) = x τ ≤ 1 and ∀q ∈ Q τ , if q = e I τ (q) = true, -∆ τ is composed of the following transitions • (s, true, true, ⊥, x τ , c τ = 0, e),

• (e,

x τ = 1, c τ = C -1, ⊥, ∅, ∅, f ),
• (e, x τ = 0, c τ < C, ⊥, ∅, ∅, p),

• (p, true, true, ⊥, x τ , ∅, e),

• (e,

x τ = 1, c τ < C -1, ⊥, x τ , c τ = c τ + 1, e),
• (f, true, true, T ?, ∅, ∅, s). In the treated problem, we consider that all the task parameters (execution time, period and deadline) are known precisely in advance. In real-time scheduling, task processing can take more or less time than expected, thus the system is evaluated according to its worst case behavior. Let us consider a non deterministic scheduling problem where the execution time of a task is not given in advance, but restricted to be bounded within an interval. Let us call a real-time task τ an uncertain real-time task if at every activation of the task, the execution time can vary within an interval [C l , C u ]. An uncertain real-time task is modeled using a timed game automata noted T GA τ . The modeling of an uncertain task τ is similar to the model of definition 3, with a distinction between controllable and uncontrollable transitions. The start and preemption transition are controlled by the scheduler, while the end transition, from state e to f , is controlled by the environment and can be taken within c ∈ [C l , C u ]. This is modeled using a guard c τ ≥ C l -1 on the end transition and a guard c τ < C u -1 on the loop execution transition. Figure 3 represents an uncertain real-time task automaton where the uncontrollable transition is represented using a dashed line.

-Σ τ = {T ?}. 2. A T (s T , s T , x T , ∅, I T , ∆ T , Σ T ): -I T (s T ) = x T ≤ T , -∆ T = {(s T , x T = T, true, T !, x T , ∅, s T )}, -Σ τ = {T !} . xτ := 0 ∧ t := 1 cτ < Cu ∧ xτ = 0 ∧ t = 0 cτ ≥ C l -1 ∧ xτ = 1 xτ := 0 xτ ≤ 1 p e s f xτ := 0 ∧ cτ := 0 ∧ t := 1 cτ < Cu -1 ∧ xτ = 1 cτ := cτ + 1 ∧ xτ := 0 ∧ t := 0

Schedulability Using Real-Time Task Automata

This section introduces some of the properties we are interested to verify from a scheduling point of view on real-time systems and exposes how to check these properties on our model using CTL [START_REF] Kozen | Logics of Programs[END_REF] properties.

The first and more natural property we want to verify on the model is the feasibility of the task set. Second, one can be interested to know if the system can be scheduled with a specific scheduling policy (e.g. RM, DM, EDF). However, the feasibility of the system is not the only property of importance in the choice of the scheduling policy. One can be interested in a lot of other characterizations for the obtained schedule: minimizing average response times, the number of preemptions, the maximum jitter or ensure the sustainability of the schedule. A feasible schedule is said sustainable regarding to a parameter if it remains feasible when this parameter is changed in a positive way. For example, if a sustainable scheduling policy regarding to tasks' cost succeed to schedule a task set Σ, it should feasibly-schedule any task set where all the tasks has the same parameters as the ones of Σ except their cost which can be lesser or equal 3 . Sustainability of scheduling policies can be studied on a lot of parameters, such as the periods, the number of tasks or the system load. We will expose how these properties can be checked on the model.

A scheduling algorithm defines the rules which control when and how transitions between different task states occur. In our approach, the scheduling algorithm corresponds to a run in the product automaton of real-time task automata.

e p

x T ≤ T T ?

T !

x T = T T ?

T ? ST OP T ? T ? xτ := 0 s xτ < 1 xτ ≤ 1 s T f Fig. 4.
Real-time task automata model. To be more readable, we omit the guards and updates represented in Figure 1.

Let Σ = {τ 1 (C 1 , D 1 , T 1 ) . . . τ n (C n , D n , T n )}
be a finite set of real-time tasks. We associate to every task τ i a real-time task automata model

(A i τ , A i T ) with A i τ (Q i τ , s i , x i τ , c i τ , I i τ , ∆ i τ , Σ i τ ) and Q i τ = {s i , e i , f i , p i }.
We use a global variable proc indicating if the processor is idle (proc = 1) or not (proc = 0). This variable is tested for every transition leading to the execution state and reset to one for every transition starting from it. In the case of multi-processor scheduling, one can use several global variables, each of them corresponding to a processor.

We introduce a new state ST OP i for each automaton A i τ , which is reached if a task misses its deadline. For every state of A i τ , a transition labeled by an input action T i ? leads to state ST OP i . Thus, if a new instance of the task is active before the execution of the previous one, the automata goes to the state ST OP i as shown by Figure 4.

Let us note (q i , v i ) a configuration of the parallel composition of A 1 τ || . . . ||A n τ s.t. q i = (q 1 i , q 2 i , . . . , q n i ) where q j i is a state of the automaton A j τ . We note q 0 the tuple where ∀j ∈ 1 . . . n q j 0 = s j . The set of clocks is augmented with a global clock t which is never reset. We note t i the valuation of the clock t in the configuration (q i , v i ) and v 0 the vector of valuations where all the clocks are set to 0. Definition 4 (Scheduling run). Let Σ((A τ1 , A T1 ), . . . , (A τn , A Tn )) be a set of real-time task automata. A scheduling run ξ :

(q 0 , v 0 ) -→ (q 1 , v 1 ) -→ • • • is a run of the parallel composition of A τ1 || . . . ||A τn . A scheduling run is feasible iff ∀(q i , v i ) ∈ ξ ∀j ∈ {1 . . . n} q j i = ST OP j .
Thus, a set τ = {τ 1 (C 1 , D 1 , T 1 ) . . . τ n (C n , D n , T n )} is schedulable if and only if there exists a feasible scheduling run in the corresponding real-time task automata model. Then, schedulability can be verified using the CTL formula 1.

φ Sched : EG¬( i=1...n ST OP i ) (1)
A scheduling algorithm can be computed using a feasible trace satisfying this formula if it exists using Algorithm 1. This algorithm takes as input a feasible scheduling run and defines a scheduling policy. The proof of termination is based on the fact that if ξ is a feasible schedule, then it exists a configuration (q i , v i ) with q i = q 0 and t i > 0 where all active tasks have terminated their execution without missing their deadline. In the case where D ≤ T , this configuration is reached at least at the hyper-period, the least common multiple of the periods of all tasks. We then just have to repeat this algorithm to obtain an infinite schedule.

Work-Conserving Scheduler Our model generates possibly idle schedules, i.e. schedules where the processor can be idle at any time. Scheduling theory often implicitly addresses problems for work-conserving schedulers (that produce schedule where the processor can be idle only when there is no ready task) because leaving the processor idle when tasks are ready seems to result in a resource wasting. But it is important to not restrict the model to work-conserving schedules since it exists task sets only feasible with idle schedulers. Moreover, the scheduling problem in idle-scheduling environment was proved undecidable without a clairvoyant algorithm.

Algorithm 1 Scheduling Run Algorithm

SchedRun[ξ]

i ← 0, ti ← 0 while qi = q0 or ti == 0 do if ∃q j i = q j i+1 then if q j i+1 = e j then execute task τj at time t = ti end if if q j i+1 = p j then preempt task τj at time t = ti end if end if i ← i + 1 end while Definition 5 (Work-Conserving Scheduling Run). A scheduling run ξ is work conserving iff ∀(q i , v) ∈ ξ if ∃q j i ∈ {s i , p i } and q j i+1 = q j i then ∃k = j s.t. q k i ∈ {e k } Work conserving runs corresponds to work-conserving scheduling algorithms. Let P, be the set of n permutations of the set {p, s, f } where the tuple (f, f, . . . , f ) is excluded. We note k i ∈ P the tuple (k 1 i , k 2 i . . . k n i ) where ∀j ∈ 1 . . . n k j i ∈ {p , s j , f j } the preemptive, start and final states of the task τ j . To compute work-conserving runs, we use the CTL formula 2

φ W C : EG¬( ki∈P (k 1 i ∧ k 2 i ∧ . . . ∧ k n i ∧ P OS)) (2)
where P OS is a state of an observer automaton, which is in state P OS if ∃j ∈ 1 . . . n s.t. x j τ > 0 and in state N P OS if ∀j ∈ 1 . . . n x j τ = 0. This formula states that it exists a run where in all the configurations containing an active task, it exists a task in its execution state.

Fixed-Priority Scheduler

We consider that the tasks τ i are sorted according to a priority function, such that τ 1 ≤ τ 2 . . . ≤ τ n . A system of task is schedulable according to a fixed priority algorithm if and only if the CTL formula 3 is satisfied.

φ F P : EG¬( i=1...n-1 (s i j=i...n e j ) ∨ i=1...n-1 (p i j=i...n e j ) ∨ ( i=1...n ST OP i )) (3)
This formula states that it exists a run where in all the configurations, a task cannot be in its execution state if a less priority task is active.

Earliest-Deadline-First Scheduler Earliest deadline first (EDF) is a dynamic priority scheduler where, if two tasks are active, the scheduler assigns the processor to the task the closest to its deadline. A system of task is schedulable according to EDF if and only if the CTL formula 4 is satisfied.

φ edf : EG¬( i=1..n j =i=1..n (s i ∧e j ∧P ij ) i=1..n j =i=1..n (p i ∧e j ∧P ij )∨( i=1..n ST OP i )) (4)
In this formula, P ij and N P ij are the states of an observer automaton, the automaton is in state

P ij if x Ti -x Tj > T i -T j and in N P ij if x i T -x i T ≤ T i -T j where x i
T and x j T are the clocks of the period automata A i T and A j T . This formula states that it exists a run where in all the configurations, a task cannot be in its execution state if a task with a closer deadline is active.

Minimal Preemption Scheduler Our approach can be used to find scheduling algorithms that minimizes the number of preemptions. This can be done simply by adding a variable that counts the number of preemptions. A task cannot be preempted no more than its duration, so a guard is used to limit the number of preemptions. A global variable preemp is used to count the total number of preemptions. To verify if a task set is schedulable with less then c preemptions, we use the CTL formula:

φ P reemp : EG¬( i=1...n ST OP i ∨ preemp > c) (5)
Sustainable Scheduler Let τ = {τ 1 ([C 1 l , C 1 u ], D 1 , T 1 ) . . . τ n ([C n l , C n u ], D n , T n
)} be a finite set of uncertain real-time tasks. We associate to every task τ i an uncertain real-time task automaton T GA i τ . A scheduling algorithm is sustainable in [C i l , C i u ]∀i ∈ 1 . . . n if all tasks do not miss there deadline whatever are the execution times. Definition 6 (Scheduling strategy). Let A be a network of n uncertain realtime task automata. A scheduling strategy is a strategy f from the set of configurations of A to the set ∆ 1 c ∪ . . . ∆ n c ∪ {λ}. A scheduling strategy is feasible iff the automata A starting at the initial configuration and supervised by f never reach a state ST OP i ∀i ∈ 1 . . . n.

In CTL, this means that the network supervised by f satisfies the safety formula:

φ sust : AG¬( i=1...n ST OP i ) (6)
Given a network of uncertain real-timed automata and a formula φ sust , finding a sustainable scheduling algorithm (in [C l , C u ]) consists in the construction of a feasible strategy if one exists. The scheduler can be computed using Algorithm 2.

Algorithm 2 Scheduling Strategy Algorithm

SchedStrategy[f ]

(q, v) ← (q0, v0), t ← 0 while q = q0 or t == 0 do while f ((q, v)) = λ or no task has finished do v = v + 1 end while if f ((q, v)) = tr ∈ ∆ j c then (q k , v k ) is the successor of (q, v) while taking the transition tr if ∃q j k = q j and q j k = e j then execute task τj at time t = t k end if if ∃q j k = q j and q j k = p j then preempt task τj at time

t = t k end if end if if a task τj has finished then let be (q k , v k ) s.t.: ∀i = 1 . . . n i = j q i k = q i , v k (x i τ ) = v k (x i τ ) and q j k = f j , v k (x j τ ) = 0 end if (q, v) = (q k , v k ) end while

Self-Suspending Tasks

Traditional real-time scheduling theory is build on the hypothesis that the tasks cannot suspend themselves [START_REF] Liu | Scheduling algorithms for multiprogramming in a hard real time environment[END_REF]. However, real-world tasks have to communicate, to synchronize and to perform external input/output operations. Each of these actions can result in a significant suspension interval for a task. It is especially true when modern architectures are considered: a task executing on a multi-core processor can have to synchronize with another one which is executing on a different core ; a task executing on an heterogeneous ISA multiprocessor system can have to migrate temporarily in order to execute specific instructions ; in the context of embedded systems where some accelerated co-processors such as Digital Signal Processors (DSPs) or Graphics Processing Units (GPUs) can be available, executing a part of a task on one of them results in a self-suspension of the task from the main processor point of view. Moreover, the durations of the suspension intervals are mostly independent of the processor (CPU) speed and so of tasks' worst case execution times. The easiest way to deal with such tasks is to account these suspension intervals as a part of the task computation time itself, and to enforce busy-waiting, wasting the gain-time offered by the hardware parallelism.

Task Model with Self-Suspensions Two task models are considered in the literature. In the simplest one, each task is characterized by the tuple τ i = (C i , E i , T i , D i ) where C i is its worst case execution time and E i its worst case suspension time (T i and D i being as usual the period and the deadline). There is no assumption on the number and time arrival of suspensions. The only information is the sum of their durations. This model was used in [START_REF] Liu | Task scheduling with self-suspensions in soft real-time multiprocessor systems[END_REF][START_REF] Liu | Real-Time Systems[END_REF]. The second is an extension where a task is characterized by the tuple τ i = (P i , T i , D i ) with P i its execution pattern. An execution pattern is a tuple

P i = (C 1 i , E 2 i , C 2 i , E 2 i , ..., C m i )
Where each C j i is a worst case execution time and each E j i a worst case suspension time. This model was used in [START_REF] Lakshmanan | Scheduling self-suspending real-time tasks with rate-monotonic priorities[END_REF][START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF].

In [START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF], authors expose three negative results on these kind of systems scheduled on-line with a non-clairvoyant algorithm. These three points can be solved using our framework:

The scheduling problem for self-suspending tasks is NP-hard in the strong sense. However, with our model we can say if a task set is feasible. T AS Fig. 5. Feasible schedule exists but neither pfp or EDF is able to find it Classical algorithms do not maximize tasks completed by their deadlines. Figure 5 shows that it exists systems for which feasible schedules exist whereas neither PFP (with all priority assignments possibilities) nor EDF are able to produce one of them. However our timed-automata-assisted-scheduler (TAS on the Figure) permits to find a feasible schedule whenever one exists.

Scheduling Anomalies can occur at run-time. suspend exactly for their worst case execution time. We can see that the system is feasible all over its hyper period. Sub-figure 6(b) exposes the schedule where the third instance of τ 1 executes with the pattern P 1 = (

C 1 1 2 , E 1 1 2 , C 2 1
). It results in a deadline missed for the τ 3 at time 49. Note that it exists a simple way to enforce the sustainability: forcing the system to simulate activity when a task completes earlier than it was supposed to. Figure 6(c) exposes the resulting schedule of this strategy on the previous example. However, this results in resource wasting. The gain time resulting from the earlier completion of a task could better be used for increase system dynamism or robustness, e.g. to compensate a cost overrun. Figure 6(d) produces a feasible schedule for the previous example where the system is never idle when work is pending. It is possible at the expense of a priority inversion at time 21. The time automata assisted scheduling can easily find such feasible schedule, where the problem is probably hard for a pure non-clairvoyant scheduler. To conclude with this example, note that an EDF schedule suffers from the exactly same anomaly. The system is feasible when the costs are fixed, but deadlines are missed (by τ 2 at time 40 and τ 1 at time 50) if the third instance of τ 1 executes according to the modified pattern.

Self-Suspending Automata To handle self-suspending scheduling problem with timed automata, we had to modify our model by adding a suspended state for tasks, however, all results presented in the previous section will remain valid with the modified model. Let τ = (P, T, D) be a self suspending task with an execution pattern P = (C 1 , E 1 , C 2 , E 2 , ..., C m ) composed of m steps. To model a self suspending task using timed automata, we add to the real-time task automaton a new state sp where the execution of a task can be suspended (see Figure 7(a)). The automaton can stay in the state exactly E time units, this is modeled using an invariant x τ ≤ E and a guard x τ = E. At each step i of the execution pattern, the variables E and C are settled to E i and C i respectively and the task finishes when all the steps has been computed (see the guard i = m on the end transition of Figure 7

Conclusion

We presented a real-time application model using timed automata and exposed how to use it to handle classical and known-unfeasible scheduling problems. The approach has been tested using the tools Uppaal [START_REF] Larsen | Uppaal in a nutshell[END_REF] and Uppaal-Tiga [START_REF] Behrmann | Uppaaltiga: Time for playing games![END_REF]. In our experiments, we focused on the efficiency of the model to construct different schedules by changing the properties to verify. We managed to use our model to test the feasibility of scheduling problem according to fixed priority and dynamic priority schedulers and to produce feasible schedules when these methods are not optimal. We tested also the model for minimal preemption scheduling and to generate sustainable schedules when they exists. All the tests have been done in both cases of work-conserving and idle scheduling. Actually, the schedules are generated by hand using the feasible trace and the winning strategy generated by the tools. Future works must propose an automatic generation of the scheduler and must examine other solutions to model the preemption or focus on optimization algorithm for state space exploration. Indeed, we used a discretization of the execution time of a task to handle preemption. This can leads to a combinatoric explosion of the model when the execution time of tasks grows.

Fig. 1 .

 1 Fig. 1. Real-time task automaton

Fig. 2 .

 2 Fig. 2. Real-time task automaton, C = 7 and preemption every 2 time units

Fig. 3 .

 3 Fig. 3. Uncertain real-time task automaton

Figure 6 Fig. 6 .

 66 Fig.6. PFP is not sustainable for self-suspending tasks: counter example on system {τ1, τ2, τ3} = {((2, 2, 4), 10, 10, High) , ((2, 8, 2), 20, 20, M edium) , (2, 12, 12, Low)}

(cτ ≥ C l - 1 ∧Fig. 7 .

 17 Fig. 7. Self suspending automata

  (a)). The timed game automaton of Figure 7(b) represents the model for an uncertain self suspending task. In this model, the execution pattern P is a set of intervals ([C 1l , C 1 u ], [E 1 l , E 1 u ], [C 2 u C 2 u ], [E 2 l , E 2 u ], ..., [C m l C m u ]) representing uncertain execution time and suspension time.

e.g. this property is not verified for most of the traditional multiprocessor scheduling algorithms