François Dubois
email: francois.dubois@math.u-psud.fr

Pierre Lallemand
email: pierre.lallemand1@free.fr

Mohamed Mahdi
email: mohamedmahdi.tekitek@fst.rnu.tn

Taylor expansion method for linear lattice Boltzmann schemes with an external force Application to boundary conditions

Keywords: bounce-back, adiabatic wall, acoustics, Poiseuille flow AMS classification: 6505, 76N15, 80A20, 82C20

In this contribution we show that it is possible to get the macroscopic fluid equations of a lattice Boltzmann scheme with an external force using the Taylor expansion method. We validate this general expansion by a detailed analysis of the boundary conditions. We derive "quartic parameters" that enforce the precision of the boundary scheme. We explicit and validate the corresponding relations for a Poiseuille flow computed with the D2Q13 lattice Boltzmann scheme.

Introduction

• The lattice Boltzmann scheme is a mesoscopic method. It deals with a small number of functions that can be interpreted as populations of fictitious particles. It has been developed more than 20 years ago by Qian et al. [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF], Succi et al. [START_REF] Succi | Lattice Gas Methods for Homogeneus and Inhomegeneus Hydrodynamics, Discrete kinetic theory, lattice gas dynamics and foundations of hydrodynamics[END_REF] among others. This method can simulate various fluid dynamics problems. The dynamics of those particles is such that time, space and momentum are discretized. They move at successive discrete times between the nodes of a regular square lattice. The velocity space is discretized by a reduced set of discrete velocities. The unknown is the distribution { f i } which is function of velocities, space location and discrete time. In each time step there are two fundamental steps : advection and collision. The advection step is linear and allows the connection of a given vertex with its neighbors. It corresponds to the characteristic method with a Courant number exactly equal to one. The collision step is nonlinear and local in space. The discrete dynamics solve a system of differential equations that converge to an equilibrium state for large time (see the details e.g. in [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF]).

• A theoretical analysis of the lattice Boltzmann scheme was proposed by Qian et al. [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF] and d'Humières [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] with a Chapmann Enskog expansion coming from statistical physics. We refer to Asinari and Ohwada [START_REF] Asinari | Connection between kinetic methods for fluid-dynamic equations and macroscopic finite-difference schemes[END_REF] for a method of analysis based on the Grad moment system. A fruitful idea developed by Junk et al [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF] and our team [START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] is to use the so-called equivalent equation method derived by Lerat-Peyret [START_REF] Lerat | Noncentered Schemes and Shock Propagation Problems[END_REF] and Warming-Hayett [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF] in the context of classical finite difference schemes. The time step is considered as an infinitesimal parameter and the finite differences are expanded into a family of equivalent partial differential operators. We observe that the emerging viscosity with the lattice Boltzmann method is in some sense a discrepancy of the scheme that can explicited.

• The main goal of this study is to extend the linearized analysis of a lattice Boltzmann scheme using the Taylor expansion method [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF][START_REF] Dubois | Equivalent partial differential equations of a lattice Boltzmann scheme[END_REF] to the case where there is an external body force. In fact few methods [START_REF] Buick | Gravity in a lattice Boltzmann model[END_REF][START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF][START_REF] Luo | Unified theory of lattice Boltzmann models for non-ideal gases[END_REF][START_REF] Martys | Evaluation of the external force term in the discrete Boltzmann equation[END_REF][START_REF] Mohamad | A critical evaluation of force term in lattice Boltzmann method, natural convection problem[END_REF] have been used to incorporate external forces in lattice Boltzmann scheme. In this study, a fraction θ (with 0 ≤ θ ≤ 1) of the force term is added before the collision step and a fraction (1 -θ ) just after. Usually θ is equal to 1 2 . In this contribution we consider a given external force variable in time and space.

• Boundary conditions are performed to impose a given velocity (or a given pressure) on a part of the physical boundary. These conditions have to be translated in terms of the lattice Boltzmann scheme. When the scheme is completely defined for internal vertices a numerical boundary conditions of bounce-back or anti-bounce-back type (see e.g. Bouzidi et al. [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF]) is equivalent to impose a zero value for some combination of the particle distributions. With the help of the Taylor expansion, it is possible to make in evidence the errors associated to several methods for one dimensional thermal problems or two dimensional fluid problem. This defines so-called "quartic parameters" (see [START_REF] Dubois | On a superconvergent lattice Boltzmann boundary scheme[END_REF]). For the general framework of this kind of methodology, we refer to Ginzburg and Adler [START_REF] Ginzburg | Boundary flow condition analysis for the three-dimensional lattice Boltzmann model[END_REF], Ginzburg et al. [START_REF] Ginzburg | Two-relaxation-time lattice boltzmann scheme: about parametrization, velocity, pressure and mixed boundary conditions[END_REF] and Ginzburg and d'Humières [START_REF] Humières | Multi-reflection boundary conditions for lattice Boltzmann models[END_REF]. Our approach is to extend our previous work [START_REF] Dubois | On a superconvergent lattice Boltzmann boundary scheme[END_REF] to three representative examples. In particular a precise choice of relaxation rates allows to spectacularly decrease the error associated to the  scheme. In the case of the Poiseuille flow the numerical evaluation of the analytical solution can even be possible with this methodology.

• In the first section, we precisely define our numerical scheme. In section 2, we extend the linearized analysis of a lattice Boltzmann scheme using the Taylor expansion method [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF]. In section 3, we study a very simple one dimensional thermal problem with a given temperature on the boundary. The existence of the source term allows an experimental validation of our methodology applied to the boundary conditions. In section 4, we study a two dimensional scheme with the D2Q9 scheme in the context of the Poiseuille flow. This flow is considered with an external force term instead of physical pressure term. Theoretical analysis allows to propose quartic values for the relaxation coefficients and for the parameter θ . These parameters are confronted to a numerical test for this model with 13 velocities in section 5. To our knowledge this kind of determination and validation has not been proposed by other teams previously.

1) Lattice Boltzmann scheme with a forcing drift

• In the following, the notation "DdQq" denotes a lattice Boltzmann scheme with d space dimensions and q velocities. Space is discretized by a regular lattice L parametrized by a spatial scale ∆x. The time step is denoted by ∆t. A numerical celerity λ is naturally defined by λ ≡ ∆x ∆t . The q discrete velocities can be written as v j = λ e j for 0 ≤ j ≤ q -1. A classical example is given by the D2Q9 scheme illustrated on the Figure 1. • The lattice Boltzmann scheme is formulated in both spaces of particles and moments m k , as proposed by d'Humières [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF]. We introduce an invertible matrix

M = (M k j ), 0 ≤ k, j ≤ q -1.
Then the moments m k for 0 ≤ k ≤ q -1 are defined by linearity:

m k ≡ q-1 ∑ j=0 M k j f j .
The moment vector m ∈ R q can be written as

(1) m = V Y , 
with two kinds of quantities: conserved moments V ∈ R N and slave moments Y ∈ R q-N . The conserved moments V are not affected by the collision step when there is no forcing term. If N = 1 there is exactly one scalar partial differential equation and when N > 1 we obtain an underlying system of N partial differential equations. The family of slave moments Y relax during the collision step towards an equilibrium distribution. An elementary time evolution of a lattice Boltzmann scheme with a force is composed by four steps (see Guo et al. [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]) : computation of the dynamical variables, collision (through simple relaxation), forcing of the conserved moments and advection.

• The computation of the dynamical variables W from the first component V of the moments (1), is given by:

(2)

W = V + θ ∆t F ,
where θ is a fixed scalar in [0, 1] and F ∈ R N is a given force term. Remark that the postprocessing is obtained from the knowledge of these dynamical variables.

• The collision and relaxation step is a redistribution of the populations { f i } at each node x. The N first moments W are not affected by the relaxation.The other "non-conserved" moments relax with a time constant τ k towards the equilibrium values m eq k . In this contribution the equilibrium values are a linear function of the dynamical variables W . We suppose that the equilibrium values Y eq of the non-conserved moments are linear function of dynamical variables. We have

(3) Y eq = E W ,
where E is a fixed matrix with q -N lines and N columns. Thus the relaxation is simply described by an ordinary differential equation:

d dt (m k -m eq k ) + 1 τ k (m k -m eq k ) = 0, for k ≥ N.
Using an explicit first order Euler scheme, we obtain the algorithm:

(4) m * k = (1 -s k ) m k + s k m eq k
, where the superscript * indicates "post-collision" state. The relaxation rates s k ≡ ∆t τ k have to satisfy 0 < s k < 2 in order to maintain stability of the scheme. With the notation (1) the relaxation step is performed in the moment space as follows:

(5)

Y * = (Id -S)Y + SY eq ,
where S is the diagonal matrix of the relaxation times s k for k ≥ N. We remark that this collision step is local in space.

• Due to the force term, the conserved moments during the collision step evolve according to (6)

V * = V + ∆t F = W + (1 -θ )∆t F.



If ( ) t denotes the transpose of the matrix ( ), the moments m * ≡ (V * , Y * ) t after collision are determined by the relations ( 5) and ( 6). After the forcing and collision steps, the particle distribution f * j is recovered by the inversion of a linear system of small size:

(7) f * j = q-1 ∑ k=0 M -1 jk m * k , 0 ≤ j ≤ q -1.
• The advection step corresponds to a method of characteristics with a Courant number equal to 1 for the advection with velocity v j . We assume here that for each node x and each velocity v j , the vertex xv j ∆t is also a node of the lattice. The particles moves from a lattice node x to either itself (for v 0 = 0) or one of the q -1 neighbors x j = x + v j ∆t for velocities v j = 0 as presented in the Figure 1. So a time step of a lattice Boltzmann scheme can be written as:

(8) f i (x, t + ∆t) = f * i (x -v i ∆t, t) , 0 ≤ i ≤ q -1 , x ∈ L , where f *
i denotes the post-collision population of particles evaluated at the end ( 7) of the previous step.

2) Taylor expansion method for the equivalent equations

• In this contribution, we extend the "Berliner version" [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF] of the Taylor expansion method in order to derive macroscopic equivalent equations when an external forcing term is present. We suppose on one hand that the dynamical variables W satisfy a partial differential equation:

(9) ∂ t W = α 1 W + ∆t α 2 W + ∆t 2 α 3 W + • • • + γ 0 F + ∆t γ 1 F + . . .
and on the other hand that the non conserved moments Y follow a dynamical expansion of the type

(10) Y = E W + ∆t β 1 W + ∆t 2 β 2 W + • • • + ∆t ρ 0 F + ∆t 2 ρ 1 F + . . .
In the equations ( 9) and ( 10), α j , γ j , β j and ρ j are space derivative operators of order j. More- over α j and γ j are N × N matrices and β j and ρ j are (q -N) × N matrices. Note here that the equation ( 9) is a system of N partial differential equations that represent the evolution of the conserved variables. It gives the macroscopic (equivalent) equation of the physical problem. The equation [START_REF] Higuera | Lattice Gas Dynamics with Enhanced Collisions[END_REF] describes the expansion of the "slave variables" Y in terms of the conserved quantities. Due to (5), we first write the collision step as follows:

Y * = (Id -S)Y + SY eq = (Id -S)Y + S E W
due to (3). Then, according to the relation (2), we have

Y * = (Id -S)Y + S E (V + ∆t θ F) = S E V + (Id -S)Y + ∆t θ S E F.
• Then the moments m * after forcing and relaxation can be presented as follows:

(11) m * = V * Y * = I 0 SE I -S V Y + ∆t I θ SE F .

 FRANÇOIS DUBOIS, PIERRE LALLEMAND AND MOHAMED-MAHDI TEKITEK

Secondly, we rewrite the scheme (8) in moment space:

V Y k (x,t + ∆t) = m k (x,t + ∆t) = ∑ j,l M k, j M -1 j,l m * l (x -v j ∆t,t) then (12) V Y k (x,t + ∆t) = ∑ j,l,µ M k, j M -1 j,l (-v j ) µ ∆t |µ| µ! ∂ µ m * l (x,t).
In the previous relation ∂ µ is a space derivation of order µ = (µ 1 , µ 2 , . . ., µ d ) ∈ N d and µ! = µ 1 ! . . . µ N !. The length |µ| of the multi-index µ is equal to ∑ j µ j . We inject the value [START_REF] Humières | Generalized Lattice-Boltzmann Equations[END_REF] of m * in the previous expression. We arrange the previous formal series as increasing powers of ∆t. We obtain the following equation:

(13) V Y (t + ∆t) = ∑ n≥0 ∆t n A n B n C n D n V Y (t) + ∑ n≥0 ∆t n+1 G n H n F.
In [START_REF] Junk | Asymptotic analysis of the lattice Boltzmann equation[END_REF], the matrices A n , B n , C n , D n , G n and H n are space derivation operators of order n. They are easy to identify with the help of the right hand side of the equation [START_REF] Humières | Multi-reflection boundary conditions for lattice Boltzmann models[END_REF]. The order zero of the development ( 13) relative to ∆t is deduced from (11) by identification:

(14) A 0 B 0 C 0 D 0 = I 0 S E I -S , G 0 H 0 = I θ S E .
• Proposition 1. Taylor expansion at order 1.

The equivalent equations and the non conserved moments for linear lattice Boltzmann scheme with external force up to first order are given by the relations

(15) ∂ t W = α 1 W + γ 0 F + O(∆t) , Y = E W + ∆t β 1 W + ρ 0 F + O(∆t 2 )
with ( 16)

α 1 = A 1 + B 1 E , β 1 = S -1 (C 1 + D 1 E -E α 1 ) , γ 0 = G 0 = I , ρ 0 = S -1 (H 0 -E γ 0 -θ S E) .
Proof of Proposition 1.

We have from (13) :

V (t + ∆t) = V + ∆t ∂ t V + O(∆t 2 ) = V + ∆t (A 1 V + B 1 Y ) + ∆t G 0 F + O(∆t 2 )
and due to [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF],

Y (t + ∆t) = Y + ∆t ∂ t Y + O(∆t 2 ) = (Id -S)Y + S E V + ∆t (C 1 V + D 1 Y ) + ∆t H 0 F + O(∆t 2 ).
Using the relation [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF] we have

V = W -θ ∆tF 
and the development at equilibrium (order zero)

Y = E W + O(∆t) ,
we obtain on one hand

∂ t W = (A 1 + B 1 E)W + G 0 F + O(∆t)
and on the other hand

Y = EW + ∆t S -1 [C 1 + D 1 E -E (A 1 + B 1 E)W ] + (H 0 -Eγ 0 -θ S E) F + O(∆t 2 ) .
By identification of the previous relations with the ansatz (9) [START_REF] Higuera | Lattice Gas Dynamics with Enhanced Collisions[END_REF], we obtain the relations (15) at first order with the coefficients α 1 , β 1 , γ 0 and ρ 0 explicited in the relations [START_REF] Luo | Unified theory of lattice Boltzmann models for non-ideal gases[END_REF].

• Proposition 2. Taylor expansion at order 2.

The equivalent equations for linear lattice Boltzmann scheme with external force, as described in first section, up to second order, are:

   ∂ t W = α 1 W + γ 0 F + ∆t α 2 W + γ 1 F + O(∆t 2 ) Y = E W + ∆t β 1 W + ρ 0 F + ∆t 2 β 2 W + ρ 1 F + O(∆t 3 ) (17)
with coefficients α 1 , β 1 , γ 0 and ρ 0 proposed in ( 16) and ( 18)

                           α 2 = B 1 β 1 + A 2 + B 2 E -1 2 α 2 1 γ 1 = B 1 ρ 0 + G 1 -θ A 1 -1 2 α 1 γ 0 + (γ 0 -2θ ) ∂ t β 2 = S -1 D 1 β 1 -E α 2 -β 1 α 1 -E α 2 1 2 +C 2 + D 2 E ρ 1 = S -1 D 1 ρ 0 + H 1 -Eγ 1 -β 1 γ 0 -1 2 Eα 1 γ 0 -1 2 E γ 0 ∂ t -ρ 0 ∂ t -θ C 1 .
Proof of Proposition 2. We first differentiate relative to time the first relation of ( 15):

∂ 2 t W = ∂ t α 1 W + F + O(∆t) = α 1 α 1 W + F + ∂ t F + O(∆t) = α 2 1 W + (α 1 + ∂ t ) F + O(∆t).
We can therefore develop the dynamical variables W at time t + ∆t up to second order. Due to the definition (2), we have:

W (t + ∆t) = V (t + ∆t) + θ ∆t F(t + ∆t) = V (t) + ∆t A 1 V (t) + B 1 Y (t) + ∆t 2 A 2 V (t) + B 2 Y (t) +θ ∆t F(t + ∆t) + O(∆t 3 ) = W (t) -θ ∆t F(t) + θ ∆t F(t + ∆t) + ∆t A 1 (W -θ ∆t F) + B 1 (E W + ∆t β 1 W + ∆t ρ 0 F) + ∆t 2 A 2 W + B 2 E W + O(∆t 3 ) .

 FRANÇOIS DUBOIS, PIERRE LALLEMAND AND MOHAMED-MAHDI TEKITEK

Then taking into account the Taylor formula at second order accuracy and the expression of ∂ 2 t W :

W + ∆t ∂ t W + ∆t 2 2 α 2 1 W + (α 1 + ∂ t ) F + O(∆t 3 ) = W + θ ∆t 2 ∂ t F + ∆t A 1 W + B 1 E W +∆t 2 -∆t F + B 1 β 1 W + B 1 ρ 0 F + A 2 W + B 2 E W + O(∆t 3 ).
After dividing by ∆t:

∂ t W = (A 1 + B 1 E)W + F + ∆t A 2 + B 2 E + B 1 β 1 -1 2 α 2 1 W + ∆t G 1 + B 1 ρ 0 -θ A 1 -1 2 α 1 + (θ -1 2 ) ∂ t F + O(∆t 2
) and the expressions of α 2 and γ 1 in ( 18) are established. To get the second relation of ( 17), we first look to the derivation relative to time of the nonconserved moments Y . We have

∂ t Y = ∂ t E W + ∆t β 1 W + ρ 0 F + O(∆t 2 ) = E α 1 W + F + ∆t (α 2 W + γ 1 F) + O(∆t 2 ) +∆t β 1 α 1 W + F + ∆t ρ 0 ∂ t F + O(∆t 2 ) = E α 1 W + F + ∆t (E α 2 + β 1 α 1 )W + (E γ 1 + β 1 + ρ 1 ∂ t ) F + O(∆t 2 ) .
We differentiate relative to time at second order:

∂ 2 t Y = ∂ t E α 1 W + F + O(∆t) = E α 1 α 1 W + F) + E ∂ t F + O(∆t) = E α 2 1 W + (E α 1 + E ∂ t ) F + O(∆t).
Then we have

SY (t) = S E V -Y (t + ∆t) -Y (t) + ∆t (C 1 V + D 1 Y + θ S E F) +∆t 2 (C 2 V + D 2 Y + H 1 F) + O(∆t 3 ) = S E (W -θ ∆t F) -∆t E α 1 W + F + ∆t (E α 2 + β 1 α 1 )W + (E γ 1 + β 1 + ρ 1 ∂ t ) F -∆t 2 2 E α 2 1 W + (E α 1 + E ∂ t ) F +∆t C 1 (W -θ ∆t F) + D 1 (E W + ∆t β 1 W + ∆t ρ 0 F) + θ S E F +∆t 2 (C 2 W + D 2 E W + H 1 F) + O(∆t 3 ) = S E W + ∆t β 1 W + ∆t ρ 0 F +∆t 2 C 2 + D 2 E + D 1 β 1 -1 2 E α 2 1 -(E α 2 + β 1 α 1 ) W + H 1 + D 1 ρ 0 -θ C 1 -1 2 (E α 1 + E ∂ t ) -(E γ 1 + β 1 + ρ 0 ∂ t ) F + O(∆t 3 )
and the expressions of β 2 and ρ 1 proposed at the second line of the relations (18) are established because γ 0 is the identity matrix.

 3) Boundary conditions for D1Q3 thermal problems

• In general there is a difference between the assumed wall position and the "measured" one, where the numerical solution vanishes. For some particular lattice Boltzmann parameters, this difference is null up to order two on ∆x. In this section, we study the impact of the knowledge given in the relations (17) [START_REF] Mohamad | A critical evaluation of force term in lattice Boltzmann method, natural convection problem[END_REF] on the analysis of the boundary conditions for the lattice Boltzmann schemes. In general, the incoming particles are given as a simple function of the outgoing ones through specular reflection (bounce-back) [START_REF] Bouzidi | Momentum transfer of a Boltzmann-lattice fluid with boundaries[END_REF] or specular "anti-reflection" (anti-bounceback) that are detailed in the corpus of the text. We first focus on the D1Q3 scheme for a scalar problem.

We introduce the following moment matrix M for the D1Q3 (see Figure 2) lattice Boltzmann scheme:

(19) M =   1 1 1 0 λ -λ -2 λ 2 λ 2 λ 2   .
We consider only one conservation law (N = 1) and we set ρ = W = f 0 + f 1 + f 2 . The equilibrium values Y eq of the two non conserved moments are given according to ( 20) With the previous lattice Boltzmann D1Q3 scheme, described by relations (2), ( 6), ( 8), ( 19) and ( 20), we have the equivalent equation

Y eq =   0 λ 2 ζ   ρ ≡ E ρ.
∂ ρ ∂t -κ ∂ 2 ρ ∂ x 2 -F = O(∆t 2 )
at second order for a stationary force F, with the thermal diffusivity

κ = 2 + ζ 3 ∆t λ 2 σ 1 with σ 1 ≡ 1 s 1 - 1 2 .

Proof of Proposition 3.

There is only one scalar conservation and W ≡ ρ. According to the proposition 2, we have the equivalent equation up to order two :

∂ t W = α 1 W + F + ∆t α 2 W + γ 1 F + O(∆t 2 ).

 FRANÇOIS DUBOIS, PIERRE LALLEMAND AND MOHAMED-MAHDI TEKITEK

The associated coefficients can be computed easily:

α 1 = 0, α 2 = λ 2 σ 1 ζ +2 3 , γ 0 = 1 and γ 1 = (θ -1
2 ) ∂ t . We have in consequence γ 1 F = 0 if ∂ t F = 0 and the equivalent equation is established.

• We consider the elementary one-dimensional Poisson problem with homogeneous Dirichlet boundary conditions:

(21) -△ u(x) = c for 0 < x < 1 , u(0) = u(1) = 0.
The solution of the problem (21) is quadratic. A uniform body force F is considered to take into account the source term c in the right hand side of the Poisson equation. We denote by x b the first node at the left in the domain ]0, 1[ and by x e a fictitious outside node (see Figure 3). The node x b is supposed to be located at a distance ∆q of the solid wall position x i (see Figure 3). We implement a so-called "anti-bounce-back" boundary condition to take into account the Dirichlet boundary conditions:

(22) f 1 (x b ,t + ∆t) = -f * 2 (x b ,t) + Φ
where Φ is a given boundary term directly related to the boundary data. Observe that this term is null if we consider e.g. homogenous Dirichlet boundary conditions. • Proposition 4. Quartic treatment of a D1Q3 boundary condition. For the above D1Q3 lattice Boltzmann scheme, we define

ξ ≡ -3 + 8σ 1 σ 2 -8σ 1 (1 -2θ ).
The numerical solution of the Poisson problem [START_REF] Warming | The modified equation approach to the stability and accuracy analysis of finite difference methods[END_REF] vanishes at third order formal accuracy for the position ∆q = ∆x 2 if the condition ξ = 0 is satisfied. Proof of Proposition 4. Using the equation [START_REF] Higuera | Lattice Gas Dynamics with Enhanced Collisions[END_REF] on the non conserved moment, we have

Y = E W + ∆t β 1 W + ∆t 2 β 2 W + ∆t ρ 0 F + ∆t 2 ρ 1 F.
Then we perform a collision step to get

Y * = (I -S)Y + S E Y + ∆t θ S E F.
With the help of the inverse moments matrix M -1 , we deduce:

f * 1 = 1 2λ 2 [2m * 2 + λ m * 1 ] and f * 2 = 1 2λ 2 [2m * 2 -λ m * 1 ].
We impose the anti-bounce-back condition on the boundary :

f 1 (x b , t + ∆t) = -f 2 (x e , t + ∆t) + Φ.
Due to the fundamental time iteration of a lattice Boltzmann scheme (8), we can write this relation under the form

f * 1 (x e ,t) = -f * 2 (x b ,t) + Φ.
Using a Taylor development of the moments, we obtain :

f * 1 (x e ) + f * 2 (x b ) = (2+ζ ) 3 ρ(x i )+ + ∆t 2 λ 2 (2+ζ ) 72 -3 + 8σ 1 σ 2 -8σ 1 (1 -2θ ) ∂ 2 ρ ∂ x 2 (x i ) + O(∆t 3 ).
We can precise that Φ = (2 + ζ ) 3 ρ(x i ) in the right hand side of (22). Then the conclusion holds if we give a null value for the second order term in the previous relation.

• Remark We have extended our result [START_REF] Dubois | On a superconvergent lattice Boltzmann boundary scheme[END_REF] for a general values of the parameter θ .

-0.004 -0.002 0 0.002

0.004 -2 -1 0 1 2 ∆ q - ∆x/2 ξ θ=0.2 θ=0.5 θ=0.7

4) D2Q9 for a linear acoustic type fluid

• This scheme is described and analyzed in [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]. The matrix M is now 9×9. For a fluid model, we have N = 3 conserved quantities: density m 0 ≡ ρ and the two components m 1 ≡ J x and m 2 ≡ J y of the momentum. The non conserved moments m k for k ≥ 3 are detailed in [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]. Their equilibrium values are classical: The relaxation rates are labelled from s 3 to s 8 with the two constraints s 5 = s 6 and s 7 = s 8 to recover an isotropic model compatible with the physics. The particle directions are represented on the figure 1. With a forcing term of the type F = (0, F x , 0), the partial equivalent equations at the order 2 are given by ( 23)

             ∂ ρ ∂t + ∂ J x ∂ x + ∂ J y ∂ y = O(∆t 2 ) , ∂ J x ∂t + c 2 s ∂ ρ ∂ x -ζ ∂ ∂ x ∂ J x ∂ x + ∂ J y ∂ y -ν ∆J x -F x = O(∆t 2 ) , ∂ J y ∂t + c 2 s ∂ ρ ∂ y -ζ ∂ ∂ y ∂ J x ∂ x + ∂ J y ∂ y -ν ∆J y = O(∆t 2 ),
where the sound speed c s , the shear viscosity ν and the bulk viscosity ζ satisfy respectively

c 2 s = λ 3 , ν = ∆t σ 7 λ 2
3 and ζ = ∆t σ 3 λ 2 3 for the D2Q9 scheme. • We consider the two-dimensional Poiseuille flow with the D2Q9 lattice Boltzmann scheme. We model this classical problem with the help of an imposed volumic force. Let Ω = [1, N x ] × [1, N y ]. We enforce periodic condition at the inlet (i = 1) and at the outlet (i = N x ). A "bounceback" boundary condition is imposed on the walls to take into account the null velocity J x = 0 on the upper and lower boundary. A uniform body force F x is applied to model the pressure gradient. We observe that for stationary problems, a lot of time iterations are necessary in order to obtain correct converged results. In that case the lattice Boltzmann method is not the most efficient taking into account the fact that due the exact advection of the lattice the resolution of the stationary problem is not a simple task. Poiseuille flow for the D2Q9 scheme. Solid wall position ∆q vs. the parameter ξ ≡ -3 + 8σ 8 σ 5 -8σ 8 (1 -2θ ). Variations of the parameter θ of the relation (2) for the forcing term (θ = 0.2 with bullets, θ = 0.5 with circles and θ = 0.7 with boxes) and of the parameter σ 5 .

• The proof of the above Proposition is elementary with traditional algebraic methods (see e.g. [START_REF] Dubois | On a superconvergent lattice Boltzmann boundary scheme[END_REF]) but quite long to develop. We have implemented it without difficulty with a formal computer software. In order to convince the reader, we test the resulting condition ξ = 0 with appropriate numerical experiments presented in Figure 6 and7. These experiments extend the ones presented in [START_REF] Dubois | On a superconvergent lattice Boltzmann boundary scheme[END_REF].

5) D2Q13 for a linearized fluid

• To construct the scheme D2Q13, we add four more velocities to the previous scheme D2Q9 (see the figure 8). The details can be found e.g. in Higuera et al. [START_REF] Higuera | Lattice Gas Dynamics with Enhanced Collisions[END_REF]. Nine moments are analogous to those proposed previously. The equilibrium values are parameterized according to m eq 3 = a 1 ρ, m eq 4 = m eq 5 = 0, m eq 6 = c 1 λ J x , m eq 7 = c 1 λ J y , m eq 8 = -65+63 c 1 24 λ J x , m eq 9 = -65+63 c 1 24 λ J y , m eq 10 = a 2 ρ, m eq 11 = a 3 ρ, m eq 12 = 0.

In order to simulate isotropic models, the relaxation rates s 3 to s 12 satisfy the following relations: s 4 = s 5 , s 6 = s 7 and s 8 = s 9 as proposed e.g. in [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF]. With the previous choices, with a forcing term F = (0, F x , 0) analogous to the one proposed previously for the D2Q9 scheme, the partial equivalent equations take the form (23). In this case the sound speed c s , the shear viscosity ν and the bulk viscosity ζ satisfy respectively c 2 s = 28+a 1 26 λ , ν = ∆t σ 5 λ 2

2 and ζ = -∆t σ 3 2+a 1 26 λ 2 . In our compuation we set a 1 = -12. Then ζ > 0. We study a Poiseuille flow with the D2Q13 lattice Boltzmann scheme in the framework of an imposed source term. The general choices are analogous to the ones we did with the D2Q9 lattice Boltzmann scheme. is null and with the above choices for the lattice Boltzmann scheme, the solid wall is located at ∆q = ∆x 2 at third order accuracy if and only if ξ = 0. This proposition has been solved with the help of formal calculus. We did various experiments analogous to those presented for the D2Q9 scheme. They are summarized in the captions of Figures 9 and10 (1 -2θ ). Variations of the parameter θ of the relation (2) for the forcing term (θ = 0.2 with bullets, θ = 0.5 with circles and θ = 0.7 with boxes) and of the parameter σ 5 .

D2Q9 Conclusion

• In this contribution, we have extended the "Berliner version" [START_REF] Augier | On rotational invariance of lattice Boltzmann schemes[END_REF] of the formal Taylor expansion method for linear lattice Boltzmann scheme to the case of a possible external force.

We have established that the quartic parameters to enforce the precision of the lattice Boltzmann scheme depend not only on the relaxation rates but also on the choice of the parameter θ for the time integration (2) of the source term.

• We have derived quartic parameters for the D2Q13 model for a Poiseuille flow. Our numerical results validate the formal Taylor development. Nevertheless, a rigorous numerical analysis is still an open question. In future work, we wish to extend the quest of quartic parameters for more general physical problems, study time dependent boundary conditions and give general boundary conditions for any position of the wall relative to the lattice.

Figure 1 .

 1 Figure 1. Stencil for the D2Q9 lattice Boltzmann scheme.

Figure 2 .

 2 Figure 2. Stencil of the lattice Boltzmann scheme D1Q3.

Figure 3 .

 3 Figure 3. "Anti-bounce-back" numerical boundary condition for the D1Q3 lattice Boltzmann scheme. The boundary node x b is located at the distance ∆q of the theoretical boundary. The incoming particles f 1 (x b ) are determined from the knowledge of the outgoing density of particles f 2 (x b ).

Figure 4 . 1 =1/ 3 ,

 413 Figure 4. Results of simulations for the solid wall position ∆q vs. the parameter ξ ≡ -3 + 8σ 1 σ 2 -8σ 1 (1 -2θ ) for the D1Q3 scheme for thermal model. We have considered variations of the parameter θ of the relation (2) for the forcing term (θ = 0.2 with bullets, θ = 0.5 with circles and θ = 0.7 with boxes) and of the parameter σ 2 relative to the last moment.

Figure 5 .

 5 Figure 5. Analogous to the figure 4. On the top: regular variation of the parameter θ for θ ∈ [0.1, 0.9]. Bottom: zoom of the top figure for θ ∈ [0.31, 0.32].

m eq 3 =

 3 -2ρ, m eq 4 = ρ, m eq 5 = -

Figure 6 .

 6 Figure 6.Poiseuille flow for the D2Q9 scheme. Solid wall position ∆q vs. the parameter ξ ≡ -3 + 8σ 8 σ 5 -8σ 8 (1 -2θ ). Variations of the parameter θ of the relation (2) for the forcing term (θ = 0.2 with bullets, θ = 0.5 with circles and θ = 0.7 with boxes) and of the parameter σ 5 .

Figure 8 .

 8 Figure 8. Stencil of the D2Q13 lattice Boltzmann scheme.

  .

Figure 9 .

 9 Figure 9. Poiseuille flow for the D2Q13 scheme. Solid wall position ∆q vs. the parameter ξ ≡ -7 + 40σ 5 σ 7 -8σ 5(1 -2θ ). Variations of the parameter θ of the relation (2) for the forcing term (θ = 0.2 with bullets, θ = 0.5 with circles and θ = 0.7 with boxes) and of the parameter σ 5 .

• Proposition 5. Precise position of the Poiseuille boundary with D2Q9. With the above choices for the lattice Boltzmann scheme, the solid wall is located at ∆q = ∆x 2 up to third order when the parameter

Aknowledgements

The authors thanks the referee for helpful questions after the first draft of this contribution.


domain: info.info-au spi sdv

text:
F Grognard

Y Chitour
email: yacine.chitour@math.u-psud.fr

G Bastin
email: bastin@auto.ucl.ac.be

EQUILIBRIA AND STABILITY ANALYSIS OF A BRANCHED METABOLIC NETWORK WITH FEEDBACK INHIBITION

Keywords: Metabolic network, small gain, stability, limit cycle

This paper deals with the analysis of a metabolic network with feedback inhibition. The considered system is an acyclic network of mono-molecular enzymatic reactions in which metabolites can act as feedback regulators on "preceding" enzymes of their own pathway, and in which one metabolite is the root of the whole network. We show, under mild assumptions, the uniqueness of the equilibrium. In the simplified case where inhibition only acts on the reactions having the root of the network as substrate, we show that the equilibrium is globally attractive. This requires that we impose conditions on the kinetic parameters of the metabolic reactions.

INTRODUCTION

The cellular metabolism is defined as the (huge) set of biochemical reactions that occur inside a living cell for growth and reproduction. It is usually represented by an intricate network connecting the involved biochemical species (called "metabolites"). The pathways of the network are called "metabolic pathways". In the metabolic engineering literature, it is widely accepted that "despite their immense complexity, metabolic systems are characterized by their ability to reach stable steady states" (quoted from [START_REF] Stephanopoulos | Metabolic Engineering: Principles and Methodologies[END_REF], Chapter 4). It should however be fair to recognize that a mathematical analysis of this fundamental stability property is a difficult question which was not much investigated. We shall limit ourselves to simple metabolic pathways which are made up of sequences of mono-molecular enzyme-catalysed reactions in the form X s → X p .

Those reactions can be inhibited by the presence of other metabolites in the network. Without this inhibition, the stability analysis is straightforward; we will then concentrate on networks with inhibition, like the one given by the aspartate amino-acid pathways [START_REF] Umbarger | Amino acid biosynthesis and regulation[END_REF], see Figure 1. In this network, each produced amino-acid inhibits an enzyme of its own pathway. This action can be seen as a negative feedback, that regulates the behavior of the network. Indeed, if we, for example, consider a large excess of isoleucine (X 20 ), the reaction X 16 -→ X 17 is shut down, so that the concentration of isoleucine is progressively reduced.

In Section 2, a model of metabolic networks such as the one of Figure 1 will be presented. The equilibria of these models will then be studied in Section 3, followed by a stability analysis in Section 4, where global attractivity of a single equilibrium is shown. x 20 Fig. 1. Metabolic network representing the aspartate amino-acid pathways: the solid lines represent the reactions and the dash-dotted lines the inhibition produced by the state at the start of the arrow onto the reaction that lies at the end of the arrow. The root of the metabolic pathway is x 1 (aspartate), and the products are the corresponding amino acids: lysine (x 9 ), methionine (x 14 ), threonine(x 16 ), and isoleucine (x 20 ).

The non-genericity of the stability of the equilibrium is then illustrated in Section 5.

MODEL OF A METABOLIC NETWORK

In our model of a metabolic network made of enzymecatalysed reactions in the form X s → X p which are inhibited by other metabolites of the network, we will denote the inhibition factor by x [p] , a n p dimensional vector containing the molar fractions of all the metabolites inhibiting the reaction X s → X p . This reaction is characterized by a velocity ϕ sp (x s , x [p] ).

On the other hand, some metabolites of the considered network are used in stages of the metabolism that are not modeled in the considered network. Therefore, we must include consumption terms in the model for those reactions in the form X s → ... We will generically denote those terms ϕ s0 (x s ).

We impose that these reaction velocities satisfy the following assumption:

Assumption 1.

• For all s, p such that the reaction X s → X p belongs to the metabolic network, the function ϕ sp (x s , x [p] ) satisfies ϕ sp (0) = 0, is non-decreasing in x s for x s ≥ 0 and non-increasing in x

[p] j . Note that this must also be valid when the value of p is 0.

In order to define the class of metabolic networks that we consider, we need the following definition from graph theory: Definition 1. A directed graph is called an arborescence if, from a given node x, known as the root node, there is exactly one elementary path from this node to any other node y.

The metabolic networks that we consider in this paper then satisfy the following assumption Assumption 2.

• the involved species are denoted

X 1 , X 2 , • • • , X n
• the graphic representation of the network (with the different metabolites as nodes and the different reactions as oriented edges) is an arborescence with X 1 as root • the inhibition acting on a reaction X s → X p only results from the action of metabolites from the (sub)-arborescence rooted in X p

Stemming from the definition and known properties of arborescences, Assumption 2 has the following consequence on the class of metabolic network that we consider:

(i) Each metabolite is produced by a single other metabolite; (ii) There is no cycle of reactions; With these definitions and notations, we shall now define a mass-balance dynamical model in the form

ẋ = Φ(x) -µx + ce 1 where x = (x 1 , • • • , x n ) T ∈ IR n ,
and x i denotes the molar fraction of the metabolite X i inside the cell. The factor µ ≥ 0 represents the specific growth rate of the cell: we assume that the cell metabolism is analyzed during a period of exponential cell growth with a constant specific growth rate µ. The vector e 1 = (1, 0, • • • , 0) T and the scalar c denote the constant supply rate of the metabolite X 1 at the root of the network. The function Φ includes all the reaction velocities, whether they correspond to reactions inside the network or reactions consuming the metabolites of the network for use in subsequent stages of the metabolism.

In order to specify Φ(x), we introduce the following notations:

Notation 1.

• P(j) = {k| the reaction X j → X k belongs to the network }. P(j) defines the set of all metabolites that are produced by reactions having X j as substrate. If there is a consumption term in the form ϕ j0 (x j ) in the derivative of x j , the index 0 is included in P(j).

• A(j) = {k|X k belongs to the arborescence with its root in X j }. "0" is not included in A(j).

It can easily be seen that P(j) \ {0} is a subsets of A(j) under Assumption 2.

From the arborescence structure, it is clear that we can separate the metabolites into three different families:

• the root X 1 : we suppose that there is a constant supply rate c of X 1 so that the corresponding mass-balance equation is the following:

ẋ1 = c - k ∈P(1) ϕ 1k (x 1 , x [k] ) -µx 1 (1)
• the intermediate metabolite X j , which is the result of the reaction X i → X j :

ẋj = ϕ ij (x i , x [j] ) - k ∈P(j) ϕ jk (x j , x [k] ) -µx j (2)
• the boundary metabolite X j (such that P(j) = {0}), which is the product of a reaction X i → X j :

ẋj = ϕ ij (x i , x [j] ) -ϕ j0 (x j ) -µx j (3)
Under Assumption 2, we can only have

x [j] = x j or x [j] = ∅.
A particular case of this network is the metabolic chain with sequential feedback inhibition (cf. Chitour, [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF])

X 1 → X 2 → • • • → X n .
with the last metabolite X n acting as an inhibitor of the first reaction X 1 -→ X 2 . In this case, all x [j] are empty, except x [2] = x n .

EQUILIBRIUM OF A METABOLIC NETWORK

Based on equations ( 1)-( 2)-(3), we can now compute the mass-balance of the whole arborescence:

d dt ( n l=1 x l ) = c -µ n l=1 x l - {k|0 ∈ P(k)} ϕ k0 (x k ) (4)
and of the arborescence that has its root in

X j d dt ( l ∈A(j) x l ) = ϕ ij (x i , x [j] ) -µ l ∈A(j) x l - {k|0 ∈ P(k) and k ∈A(j)} ϕ k0 (x k ) (5)
Those expressions will be critical in the proof of the following proposition:

Proposition 1. If Assumptions 1 and 2 are satisfied, then:

(A) the system (1)-( 2)-( 3) is positive;

(B) if all ϕ s,p (x s , x [p] ) (where p can be 0) are increasing in x s then there is at most one equilibrium

x = (x 1 , • • • , xn ) of (1)-(2)-(3) in IR n + ; (C) if µ > 0, then system (1)-(2)-(3) has a unique equilibrium x = (x 1 , • • • , xn ) in IR n + .
Moreover, the solutions of ( 1)-( 2)-( 3) are bounded for any initial condition in IR n + .

Proof: (A) is easily seen by considering the system on the boundaries of the positive orthant.

The proofs of (B) and (C) are very similar. We write the proof for (B) and highlight the differences that arise for the proof of (C) .

We will first consider system (2)-( 3) with x 1 = x1 as constant input. For any value of x1 , we will denote by (x 2 , • • • , xn ) the equilibrium of ( 2)-( 3); this equilibrium is a function of x1 , so that we will state that it is

(x 2 , • • • , xn )(x 1
). We will now show, by induction, that every element xi is an increasing function of x1 (resp. non-decreasing in case (C)).

The initial step of the proof considers the equilibrium of (3) with x i = xi as constant input

ϕ ij (x i , x[j] ) -ϕ j0 (x j ) -µx j = 0
where x[j] = xj or x [j] = ∅. When xi = 0, xj = 0 is the only solution. Also, the left-hand side of this equation is an increasing function of xi (resp. nondecreasing in case (C)) and a decreasing function of xj . It is then easily seen that, if we increase xi , xj needs also to be increased (resp.increased or kept constant) to keep this equality satisfied. We then have that, in this case, xj (x i ) is an increasing function such that xj (0) = 0 (resp. non-decreasing function such that xj (0) = 0). When µ = 0 (which can only happen in case (B)), the definition of xj could be limited to an interval [0, xm i ). Let us now make the following induction hypothesis: for a given j, the functions xk (x j ) are increasing (resp. non-decreasing) functions for all k ∈ A(j) with xk (0) = 0. We then study the equilibrium of the massbalance of the arborescence that has its root in X j . From (5):

ϕ ij (x i , x[j] (x j )) -µ l ∈A(j) xl (x j ) - {k|0 ∈ P(k) and k ∈A(j)} ϕ k0 (x k (x j )) = 0
With a similar argument to that of the initial step, we see that xj is an increasing (resp. non-decreasing) function of xi and that xj (0) = 0. The same can be said for all xk with k ∈ A j because they are already increasing (resp. non-decreasing) functions of x j .

By induction, we then have that every xk (x 1 ) is an increasing function of x1 defined on the interval [0, xm 1 ) (resp. non decreasing function defined for all x1 ). An equilibrium of the whole system then has to satisfy the equilibrium of the total mass-balance. From (4), this comes to:

µ n l=1 xl (x 1 ) + {k|0 ∈ P(k)} ϕ k0 (x k (x 1 )) = c
The system admits as many equilibria as this equation has roots. In case (B), the left-hand side is increasing from 0 when x1 increases from 0 to xm 1 (because of the second term). Therefore, if there exists an equilibrium, it is unique. In case (C), the left-hand side is strictly increasing from 0 to +∞ when x1 increases from 0 to +∞ (because of the µx 1 term), so that the equilibrium exists and it is unique.

The final point of (C) is a direct consequence of (4); this implies

d dt ( n l=1 x l ) ≤ c -µ n l=1
x l which clearly implies boundedness of the solutions when µ > 0. 2

Uniqueness of the equilibrium, especially when it is coupled with boundedness of solutions, gives hope of the possibility of having some general result about the structural global asymptotic stability of the equilibrium. In the next section, we study a special case of feedback inhibition, in which we will be able to prove global attractivity.

STABILITY OF A CLASS OF METABOLIC NETWORKS

In this section, we will study the stability of the equilibrium of a class of metabolic networks belonging to the family that was described in Section 2. In order to do that, we will first present a technical lemma.

Technical lemma

In order to analyze the stability of the considered metabolic networks, we will use the small-gain theorem of De Leenheer, Angeli & Sontag (2003) for the interconnection of monotone systems. Due to space limitation, the details of this theory are omitted here, and we only give two lemmas that we will use in the following sections (and which are proven in [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]. Let us consider the stability of a system

ẋ = f (x) x ∈ IR n (6)
which is not cooperative, but which satisfies the following assumption:

Assumption 3. Each off-diagonal element of the Jacobian J of f (x) is sign-definite (independent of x).

We then build a new system ẋ = F (x, -v) according to the following construction: for all i, j with j = i, if J ij < 0, then replace x j in f i (x) with the constant -v j . Define u 1 as the vector that contains all the constants v j that are necessary for the preceding construction (not necessarily all j have been required); we denote those constants by v s1 , • • • , v s l . From this construction, it directly appears that the system

ẋ = F (x, -u 1 ) (7)
is cooperative (in the sense defined in De Leenheer et al. (2003)) We then define the output of ( 7) as

y 1 = (x s1 , • • • , x s l )
T . The interconnection of system ( 7) with the trivial system

y 2 = u 2 (8)
through the feedback connections u 1 = -y 2 and u 2 = y 1 results in the original system (6). We can then show the following lemma Lemma 1. Suppose that the solutions of ( 6) are bounded and that system ( 7) has a unique globally attractive equilibrium for any constant input u 1 (this equilibrium defines an input output characteristic ȳ1 = k y1 (u 1 ) for system [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]). Suppose that the expression of y 1 = k y1 (u 1 ) is unknown, but that it is known that it satisfies equation

y 1 = G(y 1 , u 1 )
then system (6) has a unique globally attractive equilibrium if there exists a norm . such that

sup u1, y1 ∈ I R m ∂G ∂y 1 + ∂G ∂u 1 < 1 (9)
and if the unique equilibrium of system [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF] with u 1 constant is globally attractive for any u 1 .

Proof: see [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF] 4

.2 Feedback inhibition

The metabolic networks that we will consider for our stability study are made up of an arborescence of mono-molecular enzyme-catalysed reactions which satisfy the assumption that only one type of inhibition is present: the inhibition of the reactions that use X 1 as substrate. Also, a metabolite X m can only inhibit the reaction X 1 -→ X p , which is the first reaction of the unique elementary path linking X 1 to X m (this is a consequence of the third point of Assumption 2). This translates into Assumption 4.

• For all k ∈ P(1),

x [k] = ∅.
For simplification of notations, we will denote the elements of P(1)

as {k 1 , • • • , k r }.
A mass-balance model for such a network fits in the structure that was described in the previous section, so that we know that there is a single equilibrium. The model can be written as ( 1)-( 2)-( 3), but is particularized due to the presence of Assumption 4:

ẋ1 = c - r i=1 ϕ 1ki (x 1 , x [ki] ) -µx 1 (10) ẋki = ϕ 1ki (x 1 , x [ki] ) - j ∈P(ki) ϕ kij (x ki ) -µx ki (11) ẋk = ϕ lk (x l ) - j ∈P(k) ϕ kj (x k ) -µx k (12)
where x ki is a product of x 1 and x k of x l = x 1 . We impose the boundedness of the partial derivatives of ϕ ij in the following assumption:

Assumption 5. There exist d ij ≥ 0, α [kj ] b ≥ 0 such that 0 ≤ ∂ϕ ij ∂x i ≤ d ij for all i, j -α [kj ] b ≤ ∂ϕ 1kj ∂x [kj ] b ≤ 0 for all j ≤ r, b ≤ n j
We now define a new notation Notation 2. From the arborescence structure, we know that there exists a unique path from X 1 to any metabolite X s . This path takes the form

X 1 → X kj → • • • X k → X l → • • • → X w → X s ; if X s is
an arbitrary metabolite, we will store the indices of this path (without 1 and s) in C s ; alternatively, if x s corresponds to some x

[kj ] b , we will also denote this path C

[kj ] b . Similarly, we denote g s (k) or g [kj ] b (k) the index of the metabolite that follows X k in the path that connects X 1 to X s . This allows for the following theorem Theorem 2. If Assumptions 1, 2, 4 and 5 are satisfied and

r j=1 n k j b=1    d 1,kj µ + 1 k ∈ C [k j ] b d kg [k j ] b (k) µ + d kg [k j ] b (k)    < µ max ks,c α [ks] c (13)
and ϕ 1ki is bounded for all i ∈ {1, • • • , r} (0 ≤ ϕ 1ki ≤ B ki ), then the equilibrium of system ( 10)-( 11)-( 12) is globally attractive in the positive orthant.

Proof: The proof of this Theorem is an application of Lemma 1 (see [START_REF] References Chitour | Stability analysis of a metabolic model with sequential feedback inhibition[END_REF]). The decomposition of ( 10)-( 11)-( 12) as in ( 7) is done by replacing all the inhibiting terms x [p] with the constant vector -u 1 .

We see that, despite the fact that the inhibition is classically presented as a negative feedback that regulates the system, we have only been able to prove global attractivity under the restrictive condition (13), while stability is easily seen in the absence of inhibition. This condition is very strong, especially if the specific growth rate is small. We will analyze this condition further in the next section.

LIMIT CYCLES IN METABOLIC NETWORKS

Having obtained the sufficient result for global attractivity of Theorem 2, it is relevant to ask two questions: is condition (13) necessary and sufficient, on the one hand, and is Theorem 2 still valid without condition (13) on the other hand? The answer to the first question is easily seen to be "no": the stability of the equilibrium is retained even if condition ( 13) is slightly violated (a few simulations of simple systems is already convincing). In this section, we will show that the answer to the second question is also "no": without condition (13), the stability can be lost, so that we see that the stability of the metabolic networks is not a simple consequence of the structure of the models. Indeed, in this section, we shall exhibit an example where the equilibrium becomes unstable with a limit cycle (Hopf bifurcation) when condition ( 13) is not satisfied. We will concentrate on the stability of the equilibrium of a simple sequential pathway of four metabolites without branching and with sequential feedback inhibition (that was presented at the end of Section 2). Each metabolite produces a single other metabolite, and X 1 -→ X 2 is inhibited by the last metabolite, X 4 . We can directly apply Theorem 2 to this system:

                   ẋ1 = 1.71 - 1 1 + (x 4 /19) p 3.2x 1 1 + x 1 -0.01x 1 ẋ2 = 1 1 + (x 4 /19) p 3.2x 1 1 + x 1 - 1.4x 2 1 + x 2 -0.01x 2 ẋ3 = 1.4x 2 1 + x 2 - 1.2x 3 1 + x 3 -0.01x 3 ẋ4 = 1.2x 3 1 + x 3 - x 4 1 + x 4 -0.01x 4 (14)
where we take ϕ 1 (x 1 , x 4 ) = which is Hurwitz for p < 56.4519, and is not Hurwitz for p larger than that value. This transition from a stable to an unstable equilibrium is illustrated on Figure 2, where the time responses of the four states is illustrated for the value of p = 0 (no inhibition), p = 10 (weak inhibition), and p = 60 (strong inhibition). In the latter case, oscillations appear. This corresponds to a limit cycle in the state-space. A Hopf bifurcation has taken place in p = 56.4519. Despite this oscillation, the reaction rates for the three reactions X 2 -→ X 3 , X 3 -→ X 4 , and X 4 -→ ... are close to their maximum after the transient. The only limiting reaction is X 1 -→ X 2 which, due to the inhibition is far from its maximum reaction rate.

CONCLUSION

In this paper, we have shown that a large class of models of metabolic system only has a single equilibrium. We then have proven that, under a small gain condition, this equilibrium is globally attractive in the particular case where inhibition only acts on reaction having the root as substrate. Finally, we have shown that stability of this equilibrium is not a generic property of the metabolic systems: a condition needs to be imposed on the parameters to have stability (similar to the small gain condition that we found).

Fig. 2 .

 2 Fig.2. Evolution of the states of system (14) for p = 0 (dotted line), p=10 (dash-dotted line) and p = 60 (solid line).

  At the equilibrium, we can see that the characteristics polynomial of the Jacobian matrix has the form (s + 0.0140)(s + 0.0135)(s + 0.0130)(s + 0.0125) +4.2 10 -7 p(s + 0.01)

					1 1+(x4/19) p	3.2x1 1+x1 . The
	parameter p mainly influences the maximal slope of
	the inhibiting factor	1 1+(x4/19) p (which takes place in
	x 4 = 19). Condition (13) becomes
	α µ	d 1 µ	+ 1	d 2 µ + d 2	d 3 µ + d 3

< 1 ≡ p < 3.16 10 -5 which is a condition very similar to what was obtained in Chitour et al. (2003). This condition is very strong because it is based on a small-gain analysis. System (14) has a single equilibrium in x = (19, 19, 19, 19) T .


Alexis Heloir
email: alexis.heloir@univ-valenciennes.fr

Fabrizio Nunnari
email: fabrizio.nunnari@dfki.de

Christophe Kolski
email: christophe.kolski@univ-valenciennes.fr

Adaptive Hand-Tracked System for 3D Authoring

Keywords: Gestural input, 3D authoring ACM Classification H, 5, m, Information Interfaces and Presentation (e, g, HCI): Miscellaneous

We present the interaction design and the component architecture of an adaptive authoring system based on a consumer-range 3D input device. We claim that this system can help both novice and experienced users performing authoring tasks in a 3D authoring environment. The system uses a keyboardless self-adaptive interaction controller built upon a rule-based system that learns and infers the user's behavior/condition on the fly according to her actions; rearranging rules when necessary and suggesting breaks to avoid performance drops caused by fatigue or the so-called gorilla-arm effect.

INTRODUCTION

Animated characters are paramount in interactive media, where movement quality ad gesture repertoire has now become an intrinsic characteristic of the essence of a character, at the same level as appearance, story background and skills. Animated characters might be a privileged support of the narrative and they also have the capability to support the expression of deaf individuals on the internet, by allowing them to express opinions anonymously in their primary language: Sign Language. All these applications -customization of character's movement, animated narrative and sign language generation -share a common requirement: they all need a convenient and straightforward 3D animation authoring system.

Recent advances in consumer-range interaction devices like the Kinect 1 or the Leap Motion 2 has opened the door to an unpreceded range of new user interfaces, interaction modalities and metaphors where gesture and bodily interaction are the cornerstones [START_REF] Wigdor | Brave NUI world: designing natural user interfaces for touch and gesture[END_REF]. Taking inspiration 1. http://www.xbox.com/kinect (25 Feb 2014) 2. http://www.leapmotion.com (25 Feb 2014) from these trends, we propose a self adaptive architecture that has the potential to assist users in 3D authoring task. We currently focus on position objects in 3D space along 6 degrees of freedom. We believe that the interaction metaphors we propose in this paper can help user perform 3D editing in the following applications: character customization, animated narrative, and sign language synthesis.

Our major contribution consists of a keyboardless selfadaptive interaction controller built upon a rule-based system that learns and infers the user's behavior/condition on the fly according to her actions; rearranging rules when necessary and suggesting breaks to avoid performance drops caused by fatigue or the so-called gorilla-arm effect.

RELATED WORK

Animation is a highly cross-disciplinary domain which spans across acting, psychology, movie-making, computer graphics, and programming. Therefore, learning the art of traditional animation requires time, effort, and dedication.However, recent consumer-range technology has proved to be capable of enabling inexperienced users to either author animation of human-like bodies or to interactively control physical and digital puppets. [START_REF] Sturman | Computer puppetry[END_REF][START_REF] Sanna | A kinect-based interface to animate virtual characters[END_REF] When an animation is realized using a specialized input device or motion capture system, is often referred to as performance-driven animation or puppeteering [START_REF] Sturman | Computer puppetry[END_REF]. In performance-driven animation, the animator interactively operates a dedicated interface that allows him or her to control a synthetic character with a reinterpretation of the animator's original performance. For instance, Sanna et al. [START_REF] Sanna | A kinect-based interface to animate virtual characters[END_REF] used the Microsoft Kinect as an input-device for driving the locomotion of an animated character. Other work proved that specific input devices and interactive user interfaces can improve the authoring and animation process. On the one hand, some devices are better suited to record straight-ahead animation by capturing the motion or the dynamics of the user [START_REF] Chai | Performance animation from low-dimensional control signals[END_REF], on the other hand, other interaction devices are more suited to locate and rotate objects in a 3D space [START_REF] Lin | A sketching interface for sitting pose design in the virtual environment[END_REF]. Only a few systems are actually suited to both methods [START_REF] Nguyen | Multitouch puppetry: creating coordinated 3D motion for an articulated arm[END_REF] and to our knowledge, no proposed architecture accounts for the possibility to switch seamlessly between the two modes during an editing session. We believe that the system we present has the potential to fit both modes during an edit session. Even if, in the current paper, we focus on object displacement, the hand tracking device we currently use could also be used to record the velocity profile of the user's hand. lacking behind those for 2D graphical user interfaces. as stated in Bowman et al in 2013 [START_REF] Bowman | The Encyclopedia of Human-Computer Interaction: 3D User Interfaces[END_REF]: "There is currently no standard 3D UI (and it's not clear that there could be, given the diversity of input devices, displays, and interaction techniques), and few well-established guidelines for 3D UI design. While general HCI principles such as Nielsen's heuristics [START_REF] Nielsen | Heuristic evaluation of user interfaces[END_REF] still apply, they are not sufficient for understanding how to design a usable 3D UI." For many novice users, effective operation of 3D UIs requires going through a learning phase, even if these UIs are well designed. Even if it might seems unrealistic to completely get rid of such a learning phase, we would like to test the relevance of a self adaptive architecture for 3D UI. This architecture would be, at the beginning, very errortolerant and easy to interact with, yet rather inefficient in terms of efficiency and accuracy. The UI would however accompany the user from apprentice to mastery, reconfiguring it's interaction mechanics and thresholds, achieving better efficiency as the user gains confidence and proficiency while interacting with the UI. Previous work on adaptive UIs indeed shown that increasing accuracy led to strongly improved satisfaction. Increasing accuracy also resulted in improved performance and higher utilization of the adaptive interface. [START_REF] Gajos | Predictability and accuracy in adaptive user interfaces[END_REF] INTERACTION DESIGN With a mouse-based interface or a multitouch screen, users can control at most three to four degrees of freedom at the same time ([X, Y, scroll] or [X, Y, pinch, rotate]). In contrast, a 3D input device such as the Leap Motion gives simultaneous access to 6 degrees of freedom: translation and rotation along the three axis. We thus expect users to be faster when directly manipulating the 3D scene rather than when using the mouse and the keyboard. Indeed, for single target selection, Sears and Shneiderman [START_REF] Sears | High precision touchscreens: design strategies and comparisons with a mouse[END_REF] have shown that direct-touch outperforms the mouse.

We design our prototype using the Blender 3 3D editing tool. The standard editing mechanism of Blender uses a combination of mouse and key strokes to select the mapping between the 2 DOFs of the mouse and the spatial properties of a virtual object. The editing system is albeit quite complex, and works as follows 4 . By pressing the key G (Grab) the object location in space and the mouse position are stored. From this moment the editor enters a modal state where the 2D offset of the mouse on the real desktop is mapped to a 3D translation of the virtual object. The translation occurs on a virtual plane parallel to the screen projection plane. During the interaction, the use of keys X, Y and Z can constraint the movement of the object to a different reference axes (x,y,z) or planes (xy, xz, yz). Hitting the left mouse button terminates the modal control, sticking the object to its current location. Hitting the right mouse button, or the ESC key, cancel the modal control, bringing the object back to its initial location. This last case does not appear in the figure for clarity reasons. A similar behavior is achieved by pressing R (rotate), where the default rotation is around an axis perpendicular to the projection plane, and x,y,z constraints work similarly to the grab operation.

The interaction mechanism we implemented for the Leap Motion is much simpler: After selecting a virtual object, the user presses T to begin a simultaneous modal control for both rotation and translation, i.e 6 DOFs. We call the phase within the begin of the modal control and its confirmation/cancellation an editing action. During early trials, we observed a number of phenomena suggesting that performance editing with the HT&M approach can further improve together with the reliability of tracking technologies. For instance, when performing docking tasks, the user has a tendency to grab the virtual object as if he/she was holding a physical object.

This observations suggest that we should design an interaction mode based solely on the hand tracking device, aiming at increasing the performance and eliminating the need to learn keyboard commands. Finally, we believe that self-adaptation can be useful among different users as well as during long working sessions of a single user. In the latter case, the self-adaptation will tune the interaction mechanic of the system according to the fatigue condition of the user. The following section describes the design we propose endow self-adaptation capabilities to our system.

SELF-ADAPTIVE ARCHITECTURE

The overall architecture of our interactive system is summarized in Figure 1. It follows a feedback-controlled loop model where the user input (hand motion) is filtered out and analyzed by a component called the Motion Analyzer. This component infers a set of mid-level motion primitives and sends them to the interaction manager: a reconfigurable rule based system in charge of triggering the right interaction mode according to its input motion primitives. The interaction manager delivers a flow of edit actions. This flow is continuously analyzed by the Status and performance assessor.

Mo2on) Analyzer) • filter)bank) Interac2on) Manager) • Interac2on)state) • State)transi2on)rules)

Status)and)Performance)

Assessor)

• Ac2ons)analysis) User)

Visual)Feedback)(e.g.,)object)grabbed))

Sugges2ons) (e.g.,)"have)a)break")))

Hand)Mo2on)+) midHlevel)mo2on) primi2ves) Updated)parameter)values) Edit) Ac2ons) Figure 1.
Our architecture is inspired by a feedback controller pattern

Our system should adapt to the user, whether the user is experienced or novice. Furthermore, if the user is novice, the system should accompany the user from beginner to master level. Doing so, it should enforce best-practices and monitor user's progress while adapting to the user's performance improvement. The system should also be capable of recognizing performance level deterioration and take measures against it. Again, by enforcing good practices or by suggesting pauses. Indeed gestural interaction, when used extensively, might induce fatigue, sometimes called "gorilla arm" 5 . To avoid this, the assessor maintains a vector of descriptors that could be viewed as a simplistic user model accounting for proficiency and fatigue level. When the fatigue level increases too much, the system suggests a short break to the user. Also, the Status and Performance Assessor is continuously tuning the rules of the interaction manager to improve the user's comfort and level of performance.

Motion Analyzer

The activity of the Motion Analyzer is based on the information stream received from the hand tracking device. Each data frame contains, among others, information about the number of detected hands and a set of datastructure instances storing the position and the orientation of each palm and each detected fingertip.

The Motion Analyzer filters out and analyzes the flow of frames streamed from the Leap. It consists of a set of primitive functions performing the analysis on a timesliding window buffer of the Leap frames received in the last 2 seconds. These functions are predicates and a selection of them can be arranged according to the instructions delivered by the Status and Performance Assessor. These predicates indicate if the user's hand appeared in the current frame (newHand()), of if the hand disappeared (handLost(). Some predicates have parameters, for instace the function telling if the hand has been stable for a certain amount of time supports two parameters: time span and distance threshold (isHand-Stable[lookback time, stability radius]()). Furthermore, handFastMovement(threshold speed) detects if the user made a fast movement, which should be interpreted as moving away the hand from the editing task.

Interaction Manager

The Interaction Manager is an online reconfigurable rulebased/production system in charge of switching to the right interaction state according to the values computed by the mid-level motion primitives introduced in the previous section. The system handles three interaction states, as depicted in the right side of Fig. 1: -HOVER: this state is active when a hand has been detected by the tracking system, but no editing action is actually carried on, -GRAB: active when the user has selected an object and is moving it (visual hint), -IDLE: active when no hand is visible by the tracking device. When defined, new motion analysis primitives are bound to a selection of dynamic variables. They mostly represent time and distance thresholds. The value of these variables is updated on the fly according to how the user performs.

A selection of rules govern the basic state transition triggered by the presence/absence of the hand. When a new hand is detected, the system switches to state HOVER, regardless of its previous state. Similarly, when a hand tracking is lost, the system switches to the IDLE state. The rest of the interaction is based on the principle that when the user wants to start an editing action he needs to stabilize her hand into the sensor action space. When the hand is stable for enough time, the selected 3D object will start following the hand. This rule requires the user hand to be somehow far from the position where an object was dropped (GRAB state exited) in order to avoid undesired re-grabbing. This is done with a call to the isDistant-FromLastDrop() primitive and was inserted according to preliminary observation of users' behaviour. Stabilizing the hand again terminates the editing action.

Self-adaptivity -Status and Performance Assessor

In the previous section, we saw how the state-transition was governed by a set of rules in the Interaction Manager. These rules can be modified by reparameterizing the primitive functions composing the rules and their arrangement. Such changes have an influence on the interaction dynamics. The goal of the Status and Performance Assessor is to guarantee that the current rule arrangement and parameterization maximize the user's comfort and efficiency.

We monitor aspects of the user interaction, such as her ability to keep her hand still, to control the velocity profile of the movement, as well as the reaction time to visual cues. This architecture is user-agnostic, which means that it does not build and track a model of the user, it rather tracks the short-term evolution of the user interaction and apply adaptation strategies in order to either increase the user's efficiency or to limit the decrease of performance level.

We now present how the Status and Performance Assessor monitors the user's performance level. Table 1 lists the variables we use to measure user performances. The values of the variables are calculated through the observation of the last N o = 25 recognized actions. In the following, we describe each assessment variable and its influence on the Interaction Manager's ruleset parameterization.

Here, we apply the principle behind closed feedback loops systems, where the output (user performance) is constantly monitored the amplification (or attenuation) of the input signal (user actions) until, with a certain delay and oscillation, the system becomes stable. So, as the user performance is observed to increase, the thresholds will be decreased, to accommodate user need of system reactivity. On the contrary, if the user performance is decreasing, thresholds will be relaxed to make the system less prone to timeouts and movement range exceeding.

The value e is the exponential moving average of the last 25 edit actions. If Y i is the duration of the last edit action, e i is computed as follows: The coefficient α represents the degree of weighting decrease, a constant smoothing factor between 0 and 1. A higher α discounts older observations faster. Since α may be expressed in terms of N time periods, where α =2/(N + 1). And we want N to be around 25, we use a value of 0.077 for α.

The evolution of the resulting e value gives us a hint about the user's ability to perform faster or slower edit actions. If e decreases, we assume that the user needs a more responsive system. Therefore, we decrease the value of the GRAB START/STOP STABILITY TIME parameters (Rules 6 and 9). On the opposite, if e increases, we increase the value of the GRAB START/STOP STABILITY TIME parameters. Practically, e is used as an adjustment for a feedback gain which multiplies GRAB START/STOP STABILITY TIME by (1 + e i ) at each iteration.

m is computed exactly the same way as e on the cumulative distance which is traveled by the hand during an action. Since an experienced user is capable of accomplishing an edit with only a few large actions (GRABS) of the hand in space, a decreasing value of m indicates the user tendency to perform longer movements, including large positioning (hand is fast) and fine positioning (hand is slow) of the manipulated object. This suggests that the user is gaining in efficiency and that the interactions rules may be accommodated by decreasing the GRAB START/STOP STABILITY RADIUS required to enter and exit the GRAB status. In other words GRAB START/STOP STABILITY RADIUS are multiplied at each iteration by the following feedback gain:

(1 + m i ).

The value of c depends on how many of the last N o operations have been canceled. An operation is canceled when the user presses the ESC key in order to restore the 3D object in its initial position. For each operation with 1 <i<N o , we consider c i =0if the operation has not been canceled and c i =1id it did. We calculate the tendency by interpolating a line among the sampled results. The tendency c 0 is calculated as the tangent of the interpolated line. We consider a positive tendency as the fact that that too many operations started when the user didn't really mean to do. This will be used to increase both the GRAB START/STOP STABILITY ... TIME and RADIUS thresholds.

The value of l is calculated, similarly to c, by counting how many times the hand tracking has been lost while performing the last N o operations. A positive tendency tells us that the user hand is exiting too often from the editing space. This means that the user hardly feel uncomfortable in extreme hand positions. We use this as hint that we can increase the sensibility of the overall system, i.e., increase the ration between the quantity of motion performed by the dragged 3D object with respect to the same quantity of motion performed by the hand in real world. Consequently, also the value of DROP OFF DISTANCE (rule 6) and DIRECTION CHANGE SIZE (rule 8) are modulated.

The value of f is calculated similarly to c and l, by counting how many times the GRAB state has exited because a fast hand movement has been detected. If this occurs too frequently, the system modulates the value of the FAST MOVEMENT SPEED variable (rule 7), and may suggest the user to have a short break so that she could recover from the fatigue that might be induced by the gorilla ram effect. We are conducting further tests involving long edit session to correctly adjust this variable.

CONCLUSION AND PERSPECTIVES

To sum up, we have presented an animation system that has the potential, on the one hand, to enable novice users to author complex 3D edits and/or animations of humanoid characters and on the other hand, to increase the productivity of experienced users. This system is built upon both the Leap Motion and provides an intuitive authoring interface that has the potential to support two different edit schemes: performance capture and pose-topose animation. We further proposed a self-adaptive architecture inspired by a closed-loop controller that monitors the user's performance and tune on the fly the interaction controller in order to either increase the user's efficiency or to limit the decrease of performance. Since the device we used is commercially available at a consumerrange price and since we used open source software to build this system, we are publishing on-line 6 all the sources that are necessary to build and reproduce the described architecture. We will assess the benefit brought by such a system in further work.

  5. http://en.wikipedia.org/wiki/Touchscreen #.22Gorilla arm.22 (20 Jan 2014)

  3. http://www.blender.org, 6 Aug 2014 4. A comprehensive description of the editing basics can be found in the early tutorials of the Blender software.

Table 1 .

 1 List of variables used to assess user performance

	Name meaning
	e	action edit time
	m	action linear hand movement
	c	cancelled actions
	l	lost hand tracking
	f	fast hand movement detected

Session 4 : Techniques d'interaction: dimensions > 2 IHM'14, Villeneuve d'Ascq, France