Aymen Laadhari

Yves Barral

Gábor Székely

A data-driven optimal control method for endoplasmic reticulum membrane compartmentalization in budding yeast cells

In this paper, we propose a variational data assimilation approach for including data measurements in the simulation of the mobility of fluorescently labeled molecules in the yeast endoplasmic reticulum. The modeling framework aims to provide numerical evidence for compartmentalization in the endoplasmic reticulum. Experimental data is collected and an optimal control problem is formulated as a regularized inverse problem. To our knowledge, this is the first attempt to introduce an optimization formulation constrained by partial differential equations to study the kinetics of fluorescently labeled molecules in budding yeast. We derive the optimality conditions and use an Optimize-then-Discretize approach. A gradient descent algorithm allows accurate estimation of unknown key parameters in different cellular compartments. The numerical results support the empirical barrier index theory suggesting the presence of a physical diffusion barrier that compartmentalizes the endoplasmic reticulum membrane by limiting the exchange of proteins between the mother and its growing bud. We report several numerical experiments on real data and geometry, with the aim of illustrating the accuracy and efficiency of the method. Furthermore, a relationship between the size ratio of mother and bud compartments and the barrier index ratio is provided.

Introduction

In budding yeast cells, a small daughter cell emerges from the mother cell, which can produce nearly 20 -50 rejuvenated daughters before it dies. The endoplasmic reticulum, referred to as ER, of the budding yeast cell is made up of membranous tubules and sheet-like cisterns [Moor, 1967]. The principles underlying the anisotropic protein exchange in budding yeast have been the subject of several experimental studies showing that this is not due to the geometric organization of the ER, see [START_REF] Puhka | Endoplasmic reticulum remains continuous and undergoes sheet-to-tubule transformation during cell division in mammalian cells[END_REF], Dayel et al., 1999[START_REF] Lippincott-Schwartz | Studying protein dynamics in living cells[END_REF], Lippencott-Schwartz et al., 2003]. The endoplasmic reticulum membrane is composed of morphologically and functionally diverse domains, such as the rough ER, the smooth ER and the nuclear envelope, whose distinction is based on the contribution of several proteins [START_REF] English | Endoplasmic reticulum structure and interconnections with other organelles[END_REF]. The membrane can be divided into large compartments roughly corresponding to future daughter cells, covering the mother compartment Γ m which includes the nuclear envelope and is larger than the bud compartment Γ b ; The bud neck area Γ r , however, has a sheet morphology instead of tubules [START_REF] English | Endoplasmic reticulum structure and interconnections with other organelles[END_REF], see Fig. 1. In budding yeast cells as well as in a broad spectrum of eukaryotes [START_REF] Caudron | Septins and the lateral compartmentalization of eukaryotic membranes[END_REF], experimental studies have established that the ER membrane is physically continuous throughout the cell but that a lateral diffusion barrier may exist at the level of the bud neck and compartmentalize the membrane into an anterior and a posterior domain by limiting the exchange of ER proteins between these two domains. See e.g. [START_REF] Puhka | Endoplasmic reticulum remains continuous and undergoes sheet-to-tubule transformation during cell division in mammalian cells[END_REF], West et al., 2011]. From a biological point of view, the biologists suggested that the barrier could represent a specialized lipid domain in the cortical ER membrane at the bud neck with a different composition from the rest of the ER.

Studying the mobility of fluorescently tagged molecules in the ER membrane can provide biologists with important insights into the cellular function and organization. The Fluorescence Loss In Photobleaching technique, referred to as FLIP, is a fluorescence microscopy technique used in the laboratory to measure the mobility and molecular dynamics of proteins in living cells. It consists of repeatedly exposing mobile fluorescent molecules in a defined area of interest to an intense light pulse over time by a high-intensity laser beam, thus inducing irreversible photochemical bleaching of this area, see [START_REF] Clay | A sphingolipid-dependent diffusion barrier confines er stress to the yeast mother cell[END_REF], Lippencott-Schwartz et al., 2003]. Due to molecular mobility, an almost uniform decrease in fluorescence signal is observed over time throughout the membrane and is proportional to the ER protein concentration. Areas disconnected from the photobleached area should however continue to fluoresce. The rate at which fluorescence intensity changes after bleaching can provide experimenters with information about the movement of bleached molecules and the properties of surface diffusion in a particular region of a living cell. We notice that the amount of fluorescence measured after bleaching in each cell compartment is always normalized to the amount before bleaching.

Intriguingly, several experimental studies have shown that the kinetics of membrane proteins are very slow between the mother and bud compartments compared to the kinetics within each region where the rates of fluorescence loss are comparable. Indeed, applying FLIP in Γ m results in rapid depletion of fluorescence signal in Γ m , while fluorescence is lost only slowly in Γ b ; Conversely, applying FLIP in Γ b leads to rapid photobleaching of the marker in Γ b but not in Γ m [START_REF] Clay | A sphingolipid-dependent diffusion barrier confines er stress to the yeast mother cell[END_REF]. Accordingly, biologists rationalize that although the ER is physically continuous, its membrane should be compartmentalized, suggesting the presence of a barrier between the different membrane compartments, where protein diffusion is restricted somewhere [START_REF] Luedeke | Septindependent compartmentalization of the endoplasmic reticulum during yeast polarized growth[END_REF].

If confirmed, this theory will have a major role in maintaining and possibly optimizing the life expectancy of nascent cells. In fact, daughter cells are born with a full life expectancy, while parent cells age after each cell division. It has been shown in the state of the art that ER stress triggers mechanisms that age yeast cells. Interestingly, compartmentalization may play a key role in retaining damage and aging factors between mothers and renewing cells by preventing the aging-causing stresses from entering the daughter cells [START_REF] Tabas | Integrating the mechanisms of apoptosis induced by endoplasmic reticulum stress[END_REF], Shcheprova et al., 2008, Caudron and Barral, 2009]. This remains a very thorny area of research where several elements remain unclear. To date, it remains challenging to prove and elucidate the biological relevance of such barriers. To numerically study the properties of molecular diffusion and compartmentalization on the ER membrane, we first perform several FLIP experiments in the laboratory using confocal microscopy [START_REF] Clay | A sphingolipid-dependent diffusion barrier confines er stress to the yeast mother cell[END_REF]. In each FLIP experiment, we collected pointwise measurement data of the total amount of fluorescence separately in each cell compartment after photobleaching and at specific time instants. The experimental data will be used in the mathematical modeling framework detailed later. The data are also made available as supplementary material.

Numerical simulation tools based on mathematical modeling using optimal control have the potential to further investigate the aforementioned diffusion barrier theory. Optimal control problems governed by convection or convection-diffusion equations play an important role in many biomedical and engineering applications, see e.g [START_REF] Douglas | Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures[END_REF], Martínez et al., 2000, Zhou and Yan, 2010, Zhu and Zeng, 2003, Kim and Park, 2008, Hogea et al., 2008, Jiang and Zhang, 2000, Lenhart and Workman, 2007, Garvie and Trenchea, 2007]. Data assim-ilation formalism represents an appropriate framework for estimating model parameters, where the problem is formulated as a PDE-constrained optimization so that the key parameters maximize a performance criterion subject to an appropriate set of constraints. An appropriate cost functional measuring the discrepancy between the numerical solutions and the experimental data is minimized, while the model parameters represent the control variables that allow the observations to be "best approximated" in some sense [Gunzburger, 2002, Liao et al., 2008]. The first-order necessary conditions for optimality, called Karush-Kuhn-Tucker conditions, enable to derive an optimality system composed by "the state problem", "the adjoint problem" and "the inversion equations". In the existing literature, Lagrange finite element discretizations of optimal control problems have been widely used in such numerical frameworks. There have been several theoretical analyzes and contributions to numerical algorithms for various applications, see e.g [START_REF] Becker | Optimal control of the convection-diffusion equation using stabilized finite element methods[END_REF], Ruo et al., 2002, Yan and Zhou, 2008, Yan and Zhou, 2009, Khaksar-E Oshagh and Shamsi, 2017, Hollis et al., 1987, Garvie and Trenchea, 2014, Casas et al., 2008, Hoppe and Neitzel, 2022, Zhou and Yan, 2010, Yücel et al., 2015]. Moreover, an extensive literature on the topic of numerical methods for the optimal control of PDEs can be found in [Betts, 2020, Mills, 1983, Kelley, 1999].

In this paper, we present a deterministic mathematical modeling framework to estimate the membrane diffusion parameters using optimal control theory and numerically explore the barrier theory in yeast cells. This is part of an ongoing work to better understand the mobility of ER proteins in different biological cells, while considering both more accurate mathematical models and uncertainty quantification in solutions. The outline of this article is as follows. In section 2, we provide the mathematical framework of the transient optimal control problem and derive first-order optimality conditions. We present in the section 3 the numerical method and the algorithm. Section 4 is devoted to the numerical experiments and exploration of the diffusion barrier theory.

Fitting in Γ b

Mathematical setting

In this section, we present the mathematical framework of the data-driven model describing the dynamics of fluorescently labeled molecules on the yeast ER membrane. The acquisition of experimental data is briefly described.

Data acquisition

A series of FLIP experiments, using time-lapse laser-scanning confocal microscopy expressing the membrane marker Sec61-GFP, allows the collection of experimental data. The experiments were performed at an ambient temperature of 23 degrees Celsius as described in [START_REF] Clay | A sphingolipid-dependent diffusion barrier confines er stress to the yeast mother cell[END_REF], with some modifications as detailed in the following. Cells are imaged on the LSM 780 (Carl Zeiss, Jena, Germany) confocal microscope, while ZEN 2011 software (Carl Zeiss) is used to control the microscope. Photobleaching is applied with 100 iterations at a laser power of 100%. All intensity values were normalized to the total cell intensity and set to the value of 100% at the start of each experiment. In each FLIP experiment, pointwise measurement data of the total amount of fluorescence are collected separately in each cell compartment after photobleaching and at specific time instants. Mean values and standard deviations (SD) are calculated from repeated experiments. Figure 1 (left) shows a budding yeast cell expressing the GFP-HDEL marker through serial optical sections, in which images were processed using the software ImageJ 1.49g (National Institutes of Health).

To investigate the diffusion barrier, fluorescence signal measurement data is collected from two different sets of FLIP experiments. In the first series of experiments, photobleaching is applied in the mother cortex on Γ 1 ⊂ Γ m , while it is applied to the bud compartment on Γ 2 ⊂ Γ b in the second series of experiments, see For the numerical computations, a realistic ER geometry, referred to as Γ, is segmented based on serial optical sections through a yeast cell. Let T h be partition of Γ consisting of geometrically conforming open simplicial elements K (triangles), such that Γ = ∪ K∈T h

. Several mesh tools are used to build, improve mesh quality, and avoid badly stretched mesh elements. The initially segmented surface mesh in STL format is repaired and smoothed locally using MeshLab1 and the modules for mesh optimization in Netfabb2 . We also rely on Autodesk Maya3 for mesh quality improvement. Finally, the unstructured triangular mesh is remeshed and optimized using the Frontal algorithm of Gmsh4 , see more details in [Laadhari andSzékely, 2017, Laadhari andSzékely, 2017]. Fig. 1 shows a segmented mesh of the ER membrane with 110 760 triangular elements.

Direct problem

FLIP experiments have a timescale of minutes, while the most relevant dynamic process is fluorescent protein drain due to photobleaching. We assume, as in so many biological situations, that the fluorescently labeled molecules are free. Therefore, the FLIP recovery of the reporter protein reflects a dominant time-resolved pure surface diffusion scenario.

Let us consider a sufficiently smooth, orientable and closed ER surface Γ ⊂ R 3 , with an outward facing normal n. Let (0, T ) represent the time period of the experiment and note Γ i,T = (0, T ) × Γ i with an index i ∈ {m, b}. Let u 1 and u 2 be the unknown protein concentration fields on Γ at time t ∈ (0, T ) respectively for the series of FLIP experiments described above, i.e. when photobleaching is applied to G 1 and Γ 2 , respectively.

The optimal control formalism allows to estimate the key parameters of the model of unknown values. Indeed, the set of control variables µ m , µ b and µ r represent the surface diffusion parameters of the reporter protein respectively in the mother, bud and bud neck zones. We denote by χ i the characteristic function which serves as an indicator of the surface Γ i , with i ∈ {1, 2, m, r, b}. The total diffusion is then expressed as

µ tot = µ m χ m + µ b χ b + µ r χ r .
Essential boundary conditions must be prescribed on ∂Γ 1 in the first set of experiments and on ∂Γ 2 when photobleaching is applied to the bud. To allow the use of the same finite element mesh in both sets of FLIP experiments and avoid further mesh manipulations on the boundaries, we consider a penalty method such as a penalty parameter 1/ε helps maintain zero concentration in bleached areas, see e.g. [START_REF] Perić | Computational model for 3-D contact problems with friction based on the penalty method[END_REF], Janela, J et al., 2005, Laadhari and Quarteroni, 2016, Astorino et al., 2012] for other applications of penalty methods. We will pay attention to the setting of the penalty parameter because a choice of too large penalty can deteriorate the conditioning of the resulting linear system and induce numerical instabilities, see for example [Laadhari and Székely, 2017].

Moreover, we assume pointwise bounds on key parameters to avoid unphysical values during the optimization process. We consider an admissible set of parameters

U ad = (µ m , µ b , µ r ) ∈ R 3 : 0 < µ m A, 0 < µ b B, 0 < µ r C ,
where A, B and C are fixed values. Let I represent the identity tensor. For an arbitrary and sufficiently regular scalar function φ and vector v, we introduce the surface gradient, surface divergence, and Laplace-Beltrami operators as follows:

∇ Γ φ = (I -n ⊗ n)∇φ = ∇φ -(n.∇φ) n, div Γ v = (I -n ⊗ n) : ∇v = div v -((∇v) n).n, ∆ Γ φ = div Γ (∇ Γ φ) ,
where the tensorial product of two vectors is denoted ⊗, and the two times contracted product between tensors is denoted by the semicolon. They correspond to the two-dimensional operators evaluated in the tangent plane of Γ. The projector tensor In ⊗ n has zero eigenvalue, i.e. (In ⊗ n) n = 0, leading to non-diffusion in the direction normal to the surface. Equipped with the initial conditions u 1,0 and u 2,0 , the forward problem reads: SP(u 1 , u 2 ) : For given µ m , µ b , µ r belonging to U ad , find u 1 and u 2 such that

∂u 1 ∂t -div Γ ((µ m χ m + µ b χ b + µ r χ r ) ∇ Γ u 1 ) + 1 ε u 1 χ 1 = 0 in (0, T ) × Γ, (2.1a) ∂u 2 ∂t -div Γ ((µ m χ m + µ b χ b + µ r χ r ) ∇ Γ u 2 ) + 1 ε u 2 χ 2 = 0 in (0, T ) × Γ, (2.1b) u 1 (t = 0, •) = u 1,0 (•) in Γ, (2.1c) u 2 (t = 0, •) = u 2,0 (•) in Γ.
(2.1d)

We use the so-called optimize-then-discretize approach rather than the discretize-then-optimize approach. Indeed, the optimality conditions will first be derived and then discretized after writing the variational formulation.

Inverse problem and parameter identification

The data collection framework provides experimental measurements of fluorescence loss kinetics separately in each cellular compartment at a discrete set of time points throughout the interval (0, T ), such that the averaged protein concentration in each compartment is normalized to the concentration before photobleaching. The spatiotemporal evolution of the fluorescence decay cannot be provided by the experimental setting. We first proceed with a nonlinear least-squares curve fitting of the experimental data of the total amount of fluorescence over time, as depicted in Fig.

2.

Let F m,1 (t) = exp -9.71 × 10 -2 t 0.526 and F b,1 (t) = exp -2.10 × 10 -3 t be the time evolution functions of fluorescence loss in Γ m and Γ b , respectively, obtained by fitting after averaging twenty FLIP experiments in which photobleaching is applied to the mother at Γ 1 . Similarly, F m,2 (t) = exp -8.50 × 10 -3 t 0.758 and F b,2 (t) = exp -2.48 × 10 -1 t 0.352 denote the time evolution functions of fluorescence decay in Γ m and Γ b , respectively, obtained by fitting after averaging twenty FLIP experiments in which photobleaching is applied at Γ 2 .

The inverse problem consists in finding the optimal parameters µ m , µ b , and µ r such that the averaging solutions of the direct problem SP(u 1 , u 2 ) match as closely as possible the target averaging concentrations in each compartment. For t ∈ (0, T ), l ∈ {m, b} and i ∈ {1, 2}, and given the curve fitting functions F l,i (t), the optimization problem is expressed as:

Find optimal parameters (µ m , µ b , µ r ) = arg inf µm,µ b ,µr∈U ad J u 1 , u 2 ; µ m , µ b , µ r ,
subject to the forward problem SP (2.1a-2.1d) as a constraint.

(2.2)

The cost functional J is constructed by matching the temporal evolution of the average concentration predicted by the model in each cell compartment with the corresponding target experimental measurements. It depends on both the state variables and the controls and is expressed as:

J (u 1 , u 2 ; µ m , µ b , µ r ) = i=1,2 α m,i 2 T 0 Γm u i (t, x) - T 0 F m,i (t) Γm u 1,0 (x) 2 + i=1,2 α b,i 2 T 0 Γ b u i (t, x) - T 0 F b,i (t) Γ b u 2,0 (x + i=1,2 γ m,i 2 Γm u i (T, x) -F m,i (T ) Γm u 1,0 (x) 2 + i=1,2 γ b,i 2 Γ b u i (T, x) -F b,i (T ) Γ b u 2,0 (x) 2 + δ 2 µ m -µ b 2 + δ m 2 µ 2 m + δ b 2 µ 2 b + δ r 2 µ 2 r .
(2.

3)

The terms weighted by α ς,i , with ς ∈ {m, b} and i ∈ {1, 2}, measure the mismatch between the solution numerical and target data throughout the time interval (0, T ). However, the terms weighted by γ ς,i , with ς ∈ {m, b} and i ∈ {1, 2}, measure the discrepancy at the final instant t = T and lead to an optimal control with the desired states only at this instant. An appropriate choice of these weights, possibly zero, allows more emphasis to be placed on the solutions corresponding to the targets either throughout the duration of the experiment or at the end time, which can also be replaced by any other particular instant in (0, T ). The δ-weighted term helps account for the expectation of biologists that similar protein diffusion rates are expected in the mother and bud compartments. Finally, the terms weighted with δ j , with j ∈ {m, b, r}, are the Tykhonov regularization used to cope with the possibly illposed character of the inverse problem in the sense of Hadamard and thwart the tendency of controls to become locally unlimited. See e.g. [START_REF] Tikhonov | Numerical Methods for the Solution of Ill-Posed Problems[END_REF][START_REF] Engl | Regularization of Inverse Problems[END_REF] for a description of the mathematical theory of regularization methods. A more in-depth discussion on the choice of the regularization parameter is available in [Vogel, 2002] but is outside the scope of this work. Remark that the weights in the objective functional helps to put more emphasis on some components with respect to the other components. In practice, we either choose α ς,i = 0 or γ ς,i = 0, with ς ∈ {m, b} and i ∈ {1, 2}. In the sequel, the dependence on x and t will be omitted to alleviate the notations.

To derive the Karush-Kuhn-Tucker optimality conditions, we introduce the adjoint variables v 1 (t, x) and v 2 (t, x) of the corresponding state variables u 1 (t, x) and u 2 (t, x), respectively. The associated Lagrangian functional is expressed as:

L u i∈{1,2} ; v i∈{1,2} ; µ j∈{m,b,r} = J u i∈{1,2} ; µ j∈{m,b,r} - i∈{1,2} Γ T ∂u i dt v i + Γ T (µ m χ m + µ b χ b + µ r χ r ) ∇ Γ u i • ∇ Γ v i + Γ T 1 ε χ i u i v i .
(2.4)

The first-order optimality conditions are found by imposing the stationarity of the Lagrangian functional with respect to the adjoint, state, and inversion variables, respectively.

Let DL [ψ, δψ] be the directional derivative of L at ψ along the direction δψ. The stationarity of the Lagrangian functional with respect to the adjoint variables DL v i∈{1,2} , δv = 0 gives the direct problem SP (2.1a-2.1b). The adjoint system is obtained by imposing the stationarity of the Lagrangian functional with respect to the state variables, i.e. DL u i∈{1,2} , δu = 0. The adjoint equations read:

AP (v 1 , v 2 ):
Given the state variables u 1 and u 2 and the control variables µ m , µ b , µ r ∈ U ad , find v 1 and v 2 such that

- ∂v 1 ∂t -div Γ ((µ m χ m + µ b χ b + µ r χ r ) ∇ Γ v 1 ) + 1 ε v 1 χ 1 = ς=m,b j=1,2 α ς,j T 0 Γς u 1 - T 0 Γς F ς,j u 1,0 χ ς , in (0, T ) × Γ, (2.5a) - ∂v 2 ∂t -div Γ ((µ m χ m + µ b χ b + µ r χ r ) ∇ Γ v 2 ) + 1 ε v 2 χ 2 = ς=m,b j=1,2 α ς,j T 0 Γς u 2 - T 0 Γς F ς,j u 2,0 χ ς , in (0, T ) × Γ, (2.5b) v 1 (t = T, •) = ς=m,b j=1,2 γ ς,j Γς u 1 (T, •) -F ς,j (T ) Γς u 1,0 χ ς , in Γ, (2.5c) v 2 (t = T, •) = ς=m,b j=1,2 γ ς,j Γς u 2 (T, •) -F ς,j (T ) Γς u 2,0 χ ς , in Γ. (2.5d)
Remark that the adjoint equations (2.5a) and (2.5b) are solved backwards in time, which therefore requires terminal conditions (2.5c) and (2.5d) instead of initial conditions.

For given a, b, x ∈ R, let π [a,b] {x} = max min (x, b) , a be the projection of x on the interval [a, b]. By imposing the stationarity of the Lagrangian with respect to the inversion parameters, i.e. DL µ j∈{m,b,r} , δµ = 0, the inversion equations provide an explicit characterization of the the optimal controls as follows:

Given the state variables u 1 and u 2 and the adjoint variables v 1 and v 2 , the optimal control variables read:

µ m = π [0,A]    δ b δδ m + δ b (δ + δ m ) i=1,2 T 0 Γm ∇ Γ u i • ∇ Γ v i + δ δδ m + δ b (δ + δ m ) ς=m,b i=1,2 T 0 Γς ∇ Γ u i • ∇ Γ v i    , (2.6a) µ b = π [0,B]    δ + δ m δ µ m - 1 δ i=1,2 T 0 Γm ∇ Γ u i • ∇ Γ v i    , (2.6b) µ r = π [0,C]    1 δ r i=1,2 T 0 Γr ∇ Γ u i • ∇ Γ v i    . (2.6c)
From a numerical point of view, we rather opted for the use of a gradient optimization algorithm instead of an exact evaluation of the control parameters. This turned out to be more stable but requires the evaluation of the Lagrangian gradient against the control variables as detailed afterwards.

Numerical approximation

In this section, we describe the solution method based on a second-order time discretization and the numerical algorithm

. Let us divide [0, T ] into N + 1 subintervals [t n , t n+1 ), with n = 0, • • • , N of constant step ∆t. For any n 1, we denote by u n 1 , u n 2 , v n 1 and v n 2 the approximations of u 1 , u 2 , v 1
and v 2 at time step n , respectively. We apply a gradient descent algorithm with an adapted step length yielding a sequence of approximations to the optimal solutions and optimal key parameters. For each iteration k 0 of the gradient method, we solve the forward and inverse problems SP and AP using a fully implicit scheme. For each iteration k of the gradient method, the second-order backward differentiation formula, referred to as BDF2, is used for the time derivative terms. For any n ∈ [0, N ] and k > 0, and given the control variables µ k j with j ∈ {m, b, r}, the state unknowns u k,n+1 i with i ∈ {1, 2} are computed such that

Γ 3u k,n+1 i -4u k,n i + u k,n-1 i 2∆t ξ + Γ µ k m χ m + µ k b χ b + µ k r χ r ∇ Γ u k,n+1 i • ∇ Γ ξ + Γ 1 ε χ i u k,n+1 i ξ = 0, ∀ξ ∈ H 1 (Γ).
The numerical solution of the adjoint problem is computed backward in time starting from the final time t = T .

By change of variable τ ≡ Tt ∈ (0, T ), we obtain ∂ ∂τ = -∂ ∂t and the terminal condition becomes an initial condition. A similar numerical scheme is used in the semi-discrete time approximation of the adjoint system. For

any n ∈ [1, N + 1] and given u k,m 1 , u k,m 2 , µ k m , µ k b and µ k r , with m ∈ [1, N + 1], the semi-discrete adjoint problem consists in finding v k,n-1 i with i ∈ {1, 2} such that Γ 3v k,n-1 i -4v k,n i + v k,n+1 i 2∆t ξ + Γ µ k m χ m +µ k b χ b + µ k r χ r ∇ Γ v k,n-1 i • ∇ Γ ξ + Γ 1 ε χ i v k,n-1 i ξ = ς=m,b j=1,2 α ς,j T 0 Γς u k i - T 0 Γς F ς,j u i,0 χ ς ξ, ∀ξ ∈ H 1 (Γ).
Regarding the space discretization by finite elements, high order piecewise polynomial approximations using P κ (K), with κ 1 and K ∈ T h , are used for both state and adjoint variables.

Algorithm 1 Strategy of the optimal control problem 0: Set k ← 0, err ← 2 J 0: Let (u 1,0 , u 2,0 ) be the known initial condition, and µ 0 m , µ 0 b , µ 0 r be the initial guess

0: for n = 0, • • • , N ≡ T ∆t do 0: Solve SP u 0,n+1 1 , u 0,n+1 2 with u 0,0 i = u i,0 (initialization)
0: end for 0: Evaluate J 0 0: while err J do 0:

λ ← 6λ/5 0:

for n = N + 1, • • • , 1 do 0: Solve AP v k,n-1 1 , v k,n-1 2 using u k,N +1 1 and u k,N +1 2
for terminal conditions 0:

end for 0:

k ← k + 1 0: Compute GU µ k m , µ k b , µ k r 0: for n = 0, • • • , N do 0: Solve SP u k,n+1 1 , u k,n+1 2 with u k,0 i = u i,0 0:
end for 0:

Evaluate J u k i∈{1,2} ; µ k j∈{m,b,r} 0: while J k > J k-1 do 0: λ ← 3λ/4 0: if λ < λ then return "Algorithm stagnated" 0: end if 0: Re-compute GU µ k m , µ k b , µ k r 0: for n = 0, • • • , N do 0: Re-solve SP u k,n+1 1 , u k,n+1 2 with u k,0 i = u i,0 0:
end for 0:

Re-evaluate J u k i∈{1,2} ; µ k j∈{m,b,r} 0:

end while

0: err ← |J k -J k-1 | |J k-1 | 0: end while=0
At the computational level, we have used a gradient-based optimization algorithm because some instability issues are encountered when using an exact evaluation of the optimal controls (2.6a-2.6b-2.6c) at each optimization iteration. That requires the evaluation of the gradient of the Lagrangian with respect to the control variables. Provided with an initial guess for the controls µ 0 m , µ 0 b , and µ 0 r and a starting value for the descent length, the control variables are updated along the gradient direction using an adapted step length λ k as follows (3.1a-3.1b-3.1c).

GU µ k m , µ k b , µ k r : For k > 0, given u k-1 i , v k-1 i with i ∈ {1, 2} and µ k-1 j with j ∈ {m, b, r}, compute: µ k m = π [0,A]    µ k-1 m -λ k   δ m µ k-1 m + δ µ k-1 m -µ k-1 b - i=1,2 T 0 Γm ∇ Γ u k-1 i • ∇ Γ v k-1 i      , (3.1a) µ k b = π [0,B]    µ k-1 b -λ k   δ b µ k-1 b + δ µ k-1 b -µ k-1 m - i=1,2 T 0 Γ b ∇ Γ u k-1 i • ∇ Γ v k-1 i      , (3.1b) µ k r = π [0,C]    µ k-1 r -λ k   δ r µ k-1 r - i=1,2 T 0 Γr ∇ Γ u k-1 i • ∇ Γ v k-1 i      . (3.1c)
We proceed with an adaptation strategy for the step length λ k to ensure the decrease of the cost functional. Indeed, the step length is rejected and therefore decreased in the case where the cost functional does not decrease.

Control variables are updated only if the step length is accepted, while we exit the algorithm if λ drops below a threshold value λ . The gradient descent iterations are repeated until the relative change in the cost functional becomes smaller than a given tolerance J . Hence, convergence is achieved at iteration k > 0 if

err ≡ J k -J k-1 J k-1 < J , with J k ≡ J u k i∈{1,2} ; µ k j∈{m,b,r} .
The pseudo-code of the overall iterative scheme is detailed in Algorithm 1.

mesh size h

Numerical results

The method has been implemented using Rheolef [Saramito, 2020], which represents a C++ library for scientific computing with special emphasis on finite elements. The parallel computing relies on MPI5 and MUMPS which is used for factorization and as direct solver on distributed-memory architectures. Rheolef relies also on the libraries Boost6 , Blas7 , Scotch8 , and UMFPACK9 for much of its functionalities. Numerical results are visualized graphically using the open source application Paraview10 and Gnuplot11 . The computations are performed on a workstation with an Intel R Core TM i7-4790 (3.6 GHz) processor. In the following, we provide a set of numerical examples to evaluate the performance of the presented numerical method, with the ultimate goal of investigating numerically the barrier index theory.

Example 1

The purpose of this example is to validate the accuracy of the resolution of the direct problem, mainly the assembly of surface operators for the numerical solutions of elliptic surface partial differential equations. We consider a test case presented in [START_REF] Deckelnick | An h-narrow band finiteelement method for elliptic equations on implicit surfaces[END_REF] and study numerically the convergence properties of the solution.The problem consists in solving:

u -∆ Γ u = α 3x 2 y -y 3 in Γ, with α = - 13 8 35 π .
The geometry is the unit sphere x ∈ R 3 : |x| = 1 and the exact solution is given by:

u(x) = α |x| 2 12 + |x| 2 3x 2 y -y 3 .
We study the spatial accuracy for high order finite element approximations by computing the error in the norms L 2 (Γ), H 1 (Γ) and L ∞ (Γ) of the computed solution u h on successively refined meshes with respect to the reference exact solution π h u, where π h is the Lagrange interpolation operator in the corresponding finite element space.

The decrease in errors with respect to mesh size and convergence rates are shown in Fig. 3. Results depict the conformity between numerical and exact solutions for several polynomial finite element approximations. For example, for the Lagrange polynomials P 3 , the experimental order of convergence is 4 in the norm L 2 (Γ), whereas it is equal to 3 in the H 1 (Γ) norm. The results matches the theoretical error estimates u hπ h u 0,2,Γ < Ch κ+1 and u hπ h u 1,2,Γ < Ch κ , where C represents a constant, h = max K∈T h diameter(K) represents the mesh size, and κ is the degree of the polynomial approximation with κ ∈ {1, 2, 3, 4}, see [START_REF] Deckelnick | An h-narrow band finiteelement method for elliptic equations on implicit surfaces[END_REF].

Example 2: Optimal control simulations in simplified geometry

In this example, we proceed with a numerical validation of the proposed method in the case of simplified geometry and data. We consider an idealized geometry encompassing the convex envelopes of the mother and bud domains, and we generate successively refined semi-regular meshes. Photobleaching is only applied to the bud domain and we consider the corresponding set of experimental measurements of fluorescence loss. The different subdomains are described as regularized Heaviside functions (see for example [START_REF] Laadhari ; Laadhari | An operator splitting strategy for fluid-structure interaction problems with thin elastic structures in an incompressible newtonian flow[END_REF], Laadhari, 2017, Laadhari, 2018a]), and are highlighted in Fig. 4. Thus, the direct problem consists in finding u 2 satisfying (2.1b-2.1d), with

α m,1 = α b,1 = γ m,1 = γ b,1 = 0 in (2.3).
To test the robustness of the method for different choices of the energy parameters, we consider two different sets of parameters so that we minimize the difference between the solution and the targets either over the whole time interval or at the final instant, respectively. In the first experiment, we consider h = 0.043, ε = 10 -10 , P 1 Lagrange polynomials, α m,2 = α b,2 = 2 × 10 4 , γ m,2 = γ b,2 = 0, δ = 10 13 , δ m = δ b = δ r = 1, and µ 0 m = µ 0 b = µ 0 r = 0.1. For the optimal control with desired states at the final time, we choose h = 0.043, ε = 10 -10 , P 1 Lagrange polynomials, α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1, and µ 0 m = µ 0 b = µ 0 r = 0.05. We run the simulations for t ∈ (0, T = 349) until the convergence of the optimal control algorithm is reached. The time evolution of fluorescence loss kinetics is displayed in Fig. 6 (top) for several values of k, showing a slightly different kinetics as the matching with the data is enforced in different ways. In Fig. 6 (top right), the iterative procedure converges and allows estimating the optimal parameters so that the solution corresponds to the target kinetics at t = T . In Fig. 6 (center), we plot using a semi-logarithmic scale the decrease of the cost functional as a function of the number of iterations for the two experiments, where J δ represents the energy term weighted by δ in (2.3) (similarly for the other energy terms). We observe a first phase characterized by a rapid decrease is followed by a plateau characterizing convergence, where the optimal objective values are respectively 2.17 × 10 9 and 8.83 × 10 6 . We can also verify that the control parameters converge in Fig. 6 (bottom), showing a rapid increase from the starting values followed by a plateau characterized by µ m ≈ µ b . The estimated optimal key parameters are µ m = 0.232, µ b = 0.229, µ r = 0.124 in the first test case, while µ m = 0.195, µ b = 0.195, µ r = 0.078 in the second case. Using the optimal parameters found, Fig. 5 provides snapshots of the numerical solution at the same times for the two aforementioned cases.

Sensitivity of the optimal control solution to the photobleaching area: We now study the sensitivity of the optimal control solution to the change in the photobleaching area in the same compartment Γ b . We consider two test cases with different photobleaching areas as shown in Fig. 7. We consider a mesh size h = 0.043 and use P 1 Lagrange polynomials. The cost function parameters are set as follows: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1. Fig. 7 (top left) shows similar kinetics of fluorescence decay at convergence. Energy is minimized and convergence is achieved. The optimal model parameters are slightly different at convergence since the photobleaching is applied to different areas in Γ b . Nevertheless, we notice that we continue to have almost identical values of µ m and µ b at the optimal state in each experiment, which remains significantly higher than the diffusion parameter at the bud neck Γ r , as will be subsequently discussed.

Sensitivity of the optimal control solution to the initial guesses: We now investigate the sensitivity of the optimal control solution to the choice of starting values of the key parameters µ m , µ b and µ r . We run the simulations

α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1 and µ 0 m = µ 0 b = µ 0 r = 0.05.
for h = 0.043 using P 1 Lagrange polynomials. The cost functional parameters are set as follows: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , delta = 10 12 , δ m = δ b = δ r = 1 so that more emphasis is placed on the solution matching targets at instant t = T . For different choices of the starting values, Fig. 8 plots the changes in µ m , µ b , the objective functional J and the loss of fluorescence in Γ m and Γ b , showing that the optimal control solution is relatively insensitive to changes in the initial guesses. As expected, choosing a relatively far initialization of the model parameters results in slower convergence to the desired state. Indeed, convergence is reached after 250 iterations for µ 0 m = µ 0 b = µ 0 r = 0.05, whereas it is reached after 1500 iterations for µ 0 m = µ 0 b = µ 0 r = 1.

Convergence study for high order piecewise finite element approximations: In this test case, we study the sensitivity of the optimal parameters to the choice of the finite element approximation. We run the simulations searching for solutions in spaces generated by Lagrange polynomials of higher order P κ , κ 1. We study the estimated optimal solutions with respect to the total number of degrees of freedom, referred to as Dof, for a fixed number of mesh elements. We consider the following parameters: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 and δ m = δ b = δ r = 1. The time step size is chosen small enough not to influence substantially the overall accuracy. The time horizon is (0, T = 349). The evolution of the model parameters and the evolution of the energy cost with respect to the number of iterations are provided in Fig. 9, showing an overall convergence for higher degrees. Table 1 reports the optimal quantities calculated at convergence, i.e. when k → ∞, for higher polynomial degrees as well as the corresponding reference values, obtained by numerical continuation.

P κ Dof µ m µ b µ r J P 1 2407
Spatial convergence study: Hereafter, we carry out a quantitative convergence study with respect to the reference solution q t,ref obtained by fitting the available measurement data. Let NTS design the number of time steps in each numerical simulation. We introduce the relative error corresponding to the temporal evolution of the quantity q t and its rate of convergence as follows:

e h Γ b ,2 = NTS t=1 |q t,ref -q t,h | 2 1/2
, with q t = Γ b u 2 (t, x) dx and ROC = log 10 ( e h Γ b ,2 ) log 10 (h) .

The parameter h represents an average mesh size within a given mesh refinement level, while q t,h is an approximation of q t using T h . For different values of NTS, we appropriately consider the standard linear interpolation applied to the solution q t . The results in Table 2 show that the aforementioned quantity converges with a more than linear order of convergence in the l 2 -norm.

h Dof µ m (k → ∞) µ b (k → ∞) µ r (k → ∞) ROC 0.
P 1 Lagrange polynomials, α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1, and µ 0 m = µ 0 b = µ 0 r = 0.05.
Computational cost: To study the computational features of the proposed optimal control strategy for different choices of energy parameters, we run two different simulations with the same spatial resolution h = 0.014 using the P 1 Lagrange polynomials. We choose the parameters α ς,2 = 2 × 10 4 , γ ς,2 = 0, δ = 10 13 , δ m = δ b = δ r = 1 in the first simulation, while we set α ς,2 = 0, γ = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1, with ς ∈ {m, b}, in the second simulation. Parallel calculations are performed until convergence on a workstation using 4 cores. The computational cost is 749.4 s in the first case, while it is 636.2 s in the second case, which shows that the prescription of the matching at the final instant presents lower but not significant computational savings compared to the case where matching is imposed in the entire time interval. However, it should be noted that the two choices of energy parameter lead to qualitatively similar results in terms of fluorescence kinetics and ratios between the optimal parameters at convergence, as show above in Fig. 6. Therefore, the choice α ς,i = 0 will be preferred in the sequel.

Example 3: Diffusion barrier in yeast endoplasmic reticulum

Hereinafter, we focus exclusively on the numerical investigation of the empirical conclusions of compartmentalization by diffusion barrier, where the optimization is performed over the time interval of each FLIP experiment. We run numerical simulations using the data-driven model to investigate key model parameters in different cellular compartments, paying particular attention to the change in the cost functional. However, the effect of compartment size ratio on the kinetics of fluorescence will be subsequently explored in Example 4.4.

In this regard, we consider the setting and measurements of the aforementioned series of FLIP experiments described in subsection 2.1, where photobleaching is applied either at the mother or bud compartments. A unstructured realistic ER mesh of nearly regular triangular elements is considered. The initial guess of the key parameters to be estimated iteratively is µ 0 m = µ 0 b = µ 0 r = 0.43, while we choose ε = 10 -10 and δ m = δ b = δ r = 1.

Iteration k

Cost functional J k 5 × 10 3 2.5 × 10 3 0 1.5 × 10 14 , 15, 20, 50, 80, 110, 250, 355}. The optimal control algorithm provided satisfactory results illustrating the convergence of the iterative procedure. Moreover, we follow the decay of the objective functional J k with respect to the number of descent gradient iterations in Fig. 11 (left). The graph shows a rapid decrease in cost in an initial phase, followed by a horizontal plateau reached when the minimum is reached; That corresponds to the optimal objective functional obtained with the estimated optimal parameters. The evolution of the adapted step length λ k with respect to the number of iterations is reported in Fig. 11 (right). It shows in a first phase several fluctuations due to the adaptation of the step length to ensure the decrease of the cost functional, featuring here a rapid decrease. The variations of J are small thereafter, leading to lower and lower values of λ k until it becomes lower than the tolerance threshold λ when convergence is reached.

The change of model parameters with respect to the number of iterations is provided in Fig. 13. At convergence, the estimated model parameters are µ m = 0.661853, µ b = 0.661895 and µ r = 0.0312565. The direct problem is then solved for both FLIP experiments using the estimated optimal parameters. We evaluate the normalized intensity values throughout the time periods in both bleached and unbleached cell compartments. Fig. 14 experimental data of the fluorescence signal at specific measurement times, along with its maximum and minimum values. It depicts a good agreement between the numerical results and the experimental measurements. In figures 10 and 12, we visualize the fluorescence intensity field for eight consecutive snapshots at the optimal state after applying a photobleaching to Γ m and Γ b , respectively. In both cases, the effect of the diffusion barrier at the bud neck region can be seen.

Remark that the estimated optimal control parameters satisfy µ m /µ b ≈ 1.0 and µ m /µ r ≈ 21.2. That is, the numerical results pinpoint that the diffusion of the reporter protein is similarly fast in the mother and in the bud compartments. However, it is 21 times slower in the barrier zone at the level of the bud neck. This is in agreement with the conclusions of many biologists that a physical diffusion barrier exists between the main compartments, whereas the protein Sec61-GF moves at similar speeds in the mother and bud compartments, see e.g. [START_REF] Luedeke | Septindependent compartmentalization of the endoplasmic reticulum during yeast polarized growth[END_REF]. Consequently, the restricted diffusion in the bud neck region helps to slow the loss of fluorescence in the unbleached compartment in each FLIP experiments. From a biological point of view, the bud neck area features a specific morphology and contributes primarily, and among other things, to the retention of damage and aging factors between the mother and daughter cell, as mentioned above. Although our numerical study is in quantitative agreement with the diffusion barrier theory, we are still far from a high-fidelity mathematical description characterizing the intrinsic processes underlying such barriers.

Example 4: Barrier index and dependence on compartment size ratio

Last but not least, we explore the effect of the compartment size on the kinetics of fluorescence. An index of great importance to biologists is the Barrier Index, referred to as BI. We first define T 1/2 as the time corresponding to a concentration equal to 50% of the initial concentration before photobleaching in a given compartment, as shown graphically in Fig 14 . The BI represents the ratio of the times required to lose 50% of the fluorescence signal in the unbleached compartment compared to the bleached compartment.

An agreement in the calculation of the barrier index between the photobleached and unbleached compartments between the numerical and experimental results is obtained. Remark that the BI reaches much higher values (18.22± 1.90 experimentally and 18.572 numerically) when photobleaching is applied to the bud than when it is applied to the mother (7.76±2.27 experimentally and 6.907 numerically), see Fig 15 . The biologists rationalized that this might simply be due to the fact that the ER volume in the bud is smaller than that in the mother cell, leading to a much more rapid depletion of the smaller compartment upon FLIP.

In order to test whether this explanation indeed explains the yeast observations, the dependence between the compartment size ratio and the barrier index is studied numerically. We calculate how the observed BI would change when the relative volume of mother and bud is allowed to vary. Using the segmented ER geometry, we generate different meshes without changing the ER organization in the mother cortex but where the bud is increasingly bigger, until reaching almost the size of the mother compartment, see Fig 16 . The surface meshes are subsequently remeshed in order to optimize the quality of the mesh for carrying out finite element calculations.

We consider the optimal parameters µ m , µ b and µ r found previously. Parallel FLIP experiments were then simulated in each of these meshes, FLIPing either the mother or the bud compartment. In each geometry, the values of the BI obtained by photobleaching in the mother and in the bud are then compared. Whereas the ratio between the BI in the bud divided by the BI in the mother is high as long as the bud domain is significantly smaller than the mother cell, this ratio becomes close to 1 as the size of the bud approaches the size of the mother, see Fig. 17. A noticeable consequence is that, at least in budding yeast cell, the difference in BI depending on the place of photobleaching reflects the fact that (i) the ER membrane is compartmentalized and (ii) one compartment is smaller than the other.

Conclusions

This contribution presents a numerical framework for the modeling of the kinetics of fluorescently tagged molecules on the endoplasmic reticulum in asymmetrically dividing yeast cell. Using the Fluorescence Loss In Photobleaching technique, experiments are carried out in laboratory and provide measurement data on the decrease in the level of fluorescence in the different cellular compartments. To study the anisotropic molecular diffusion, we present a datadriven mathematical model based on the use of partial differential equations constrained optimization. Optimality conditions are derived and a gradient descent algorithm is used to estimate the diffusion parameters in the different cellular compartments. We address the main features of the method and we provide simulations with the aim of providing numerical evidence of compartmentalization in the ER membrane. Our computational model supports some biological conclusions indicating that the exchange of membrane proteins between mother and bud compartments is very slow, almost 20 times slower, compared to the diffusion of proteins into each compartment due to some surface diffusion limitations at the bud-neck zone. This compartmentalization is also conserved in a broad range of cellular contexts, such as mouse neural stem cells and the early C. Elegans embryo, and contributes to processes as diverse as the confinement of protein aggregates during aging and the patterning of developing embryos. Further improvements are needed to better explore the compartmentalization in budding cells, while accounting for high fidelity mathematical descriptions. This is part of an ongoing work to explore concordant hypotheses for anisotropic exchange of ER proteins in yeast and and other cells [START_REF] Laadhari | Fully eulerian finite element approximation of a fluid-structure interaction problem in cardiac cells[END_REF], Laadhari, 2011], with more accurate mathematical formalism accounting for uncertainty quantification in the solution and the sensitivity of data noise on model parameters.

Figure 1 :

 1 Figure 1: Confocal images of FLIP experiments applied to mother (top) or bud (bottom) domains expressing the membrane marker Sec61-GFP. The gray bars indicate the first pre-bleached images, while the black bars indicate snapshots at specific instants of repeated bleaching in the rest of the experiment. Photobleached areas are highlighted by white rectangles. (Right) Segmented geometries showing: Γ m (green), Γ b (yellow), Γ r (red) and photobleached areas Γ 1 ∈ Γ m (top, blue) or G 2 ∈ Γ b (bottom, blue).

FittingFigure 2 :

 2 Photobleaching in mother

  Fig 1. Fig 2 reports the decay the amount of fluorescence signal over time in the mother and bud compartments. Each experiment was repeated on average twenty times and the error bars indicate the standard deviation of the data collected.

  For any n ∈ [1, N + 1] and k > 0, the approximations of the state unknowns u k,n 1 and u k,n 2 are computed by induction, using values at previous time steps. Similarly, for any n ∈ [0, N ], the adjoint unknowns v k,n 1 and v k,n 2 are computed by induction, using values at next time steps. For i ∈ {1, 2}, the scheme is bootstrapped by initial conditions (2.1c-2.1d) u -1 i = u 0 i = u i,0 and terminal conditions (2.5c-2.5d) v N +2 i = v N +1 i = v i (T ), where u -1 i and v N +2 i only stand for convenient notations.

Figure 3 :

 3 Figure 3: Example 1: Study of the spatial convergence in the norms L 2 (Γ) (top), H 1 (Γ) (bottom left) and L ∞ (Γ) (bottom right) for high order Lagrange polynomials P κ and high order curved surface meshes. The logarithmic scale is used.

Figure 4 :

 4 Figure 4: Example 2: Idealized and simplified geometry of yeast ER highlighting in red the mother, bud, bud neck and photobleaching areas, respectively.

Figure 5 :

 5 Figure 5: Example 2: Snapshots showing the computed fluorescence loss kinetics at optimal state for times t ∈ {10, 70, 170, 340} with a mesh h = 4 × 10 -3 . Top: α m,2 = α b,2 = 2 × 10 4 , γ m,2 = γ b,2 = 0. Bottom: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 .

Figure 6 :

 6 Figure 6: Example 2: Numerical results of the optimal control with the desired states enforced over the entire time period (Left: γ m,2 = γ b,2 = 0) and with the desired states at the final time T (Right: α m,2 = α b,2 = 0). Top: Loss of fluorescence over time at certain iterations towards convergence. The arrows show the convergence towards the optimal control solution for certain values of k between 0 and 3500. Center: Minimization of the functional cost. Bottom: Convergence of key model parameters.

Figure 7 :

 7 Figure 7: Example 2: Test cases, referred to as TC-1 and TC-2, with different photobleaching areas colored in red (top left). Time evolution of the fluorescence decay (top right). Change in optimal parameters (bottom left) and energy decay (bottom right).

Figure 8 :

 8 Figure 8: Example 2: Modification of estimates of control variables and optimal control solution for different choices of starting values for key parameters. Simulation parameters: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1.

Figure 9 :

 9 Figure 9: Example 2: Change in the estimates of control variables and cost functional for high order finite elements. Parameters: α m,2 = α b,2 = 0, γ m,2 = γ b,2 = 10 8 , δ = 10 12 , δ m = δ b = δ r = 1 and µ 0 m = µ 0 b = µ 0 r = 0.05.

Figure 10 :

 10 Figure 10: Example 3: Flipping in mother compartment. Snapshots showing fluorescence kinetics loss in numerical simulations with optimal model's parameters at times t ∈ {5, 15, 20, 50, 80, 110, 250, 355}.

Figure 11 :

 11 Figure 11: Example 3: (Left) Convergence of the objective functional with respect to the iterations count k. (Right) Change in the step length with respect to the iteration count. The logarithmic scale is used on the y axis. Parameters: δ = 10 6 , α ς,j = 10 6 and γ ς,j = 0 for ς ∈ {m, b} and j ∈ {1, 2}, and δ m = δ b = δ r = 1.

Figure 12 :

 12 Figure 12: Example 3: Flipping in bud compartment. Snapshots showing fluorescence kinetics loss in numerical simulations with optimal model's parameters at times t ∈ {5, 15, 20, 50, 80, 110, 250, 355}.

Figure 13 :

 13 Figure 13: Example 3: Change in the estimates of control variables versus the iteration count k. Scaling coefficients: δ = 10 6 , α ς,j = 10 6 and γ ς,j = 0 for ς ∈ {m, b} and j ∈ {1, 2}, and δ m = δ b = δ r = 1.

Figure 14 :

 14 Photobleaching in mother

Figure 15 :

 15 Figure 15: Example 4: Bar chart representing the T 1/2 of the unbleached compartment over the T 1/2 of the bleached compartment in the experimented in comparison to the simulated data. Mean ± SD.

Figure 16 :

 16 Figure 16: Example 4: (Left) ER mesh featuring a size of the bud nearly identical to the size of the mother cell. (Right) Segmented realistic realistic geometry of the yeast ER.

Figure 17 :

 17 Figure 17: Example 4: Change in the ratio BI bud /BI mother when the size of the bud increases relative to the size of the mother. It tends to 1 for similar sizes of ER compartments.

Table 1 :

 1 Example 2: Convergence history for high order finite element approximations using P κ with κ 1.

			0.195178 0.194918 0.0784062 8.83096E+6
	P 2	9622	0.201067 0.200895 0.0977545 4.98737E+6
	P 3	21647	0.202328 0.202163 0.102388 3.89932E+6
	P 4	38492	0.203342 0.203204 0.101571 2.89142E+6
	P 5	60127	0.204080 0.203959 0.103173 2.23366E+6
	Reference values 0.206645 0.206725 0.102447 1.51009E+6

Table 2 :

 2 Example 3: Convergence history of same outputs with respect to spatial resolution. Simulation parameters:

	094	1450	0.17493	0.17427	0.07649	2.77
	0.060	2634	0.18637	0.18594	0.06838	2.33
	0.043	4810	0.19538	0.19512	0.07592	2.08
	0.038	7262	0.19930	0.19912	0.07714	1.85
	0.029	12556	0.20756	0.20755	0.08490	1.65
	0.019	26692	0.21499	0.21509	0.09372	1.53
	0.014	28348	0.21669	0.21682	0.09342	1.53
	0.009	102084	0.22317	0.22343	0.10016	1.37
	0.007	150522	0.22501	0.22531	0.10168	1.31
	0.004	368908	0.22799	0.22836	0.10386	1.19
	Reference values	0.2305	0.2309	0.1084	-

MeshLabhttp://meshlab.sourceforge.net

Netfabbhttp://www.netfabb.com

Mayahttp://www.autodesk.com/products/maya

Gmshhttp://www.geuz.org/gmsh

Message Passing Interfacehttp://www.mpich.org

Boost librarieshttp://www.boost.org

Basic Linear Algebra Subprograms libraryhttp://www.netlib.org/blas

Scotchhttp://www.labri.fr/perso/pelegrin/scotch

Umfpack routineshttp://www.cise.ufl.edu/research/sparse/umfpack/

Paraviewhttp://www.paraview.org

Gnuplothttp://www.gnuplot.info

Acknowledgement

Conflict of interest

The authors declare no potential conflict of interests.


¸ois Franc

Jean-Patrick Rossigneux

Laurent Gelas

Marcos Lefèvre

Dias De Assunc

A Generic and Extensible Framework for Monitoring Energy Consumption of OpenStack Clouds

Although cloud computing has been transformational to the IT industry, it is built on large data centres that often consume massive amounts of electrical power. Efforts have been made to reduce the energy clouds consume, with certain data centres now approaching a Power Usage Effectiveness (PUE) factor of 1.08. While this is an incredible mark, it also means that the IT infrastructure accounts for a large part of the power consumed by a data centre. Hence, means to monitor and analyse how energy is spent have never been so crucial. Such monitoring is required not only for understanding how power is consumed, but also for assessing the impact of energy management policies. In this article, we draw lessons from experience on monitoring large-scale systems and introduce an energy monitoring software framework called KiloWatt API (KWAPI), able to handle Open-Stack clouds. The framework -whose architecture is scalable, extensible, and completely integrated into OpenStack -supports several wattmeter devices, multiple measurement formats, and minimises communication overhead.

I. INTRODUCTION

Cloud computing [START_REF] Armbrust | Above the clouds: A Berkeley view of Cloud computing[END_REF] has become a key building block in providing IT resources and services to organisations of all sizes. Among the claimed benefits of clouds, the most appealing derive from economies of scale and often include a pay-asyou-go business model, resource consolidation, elasticity, good availability, and wide geographical coverage. Despite these advantages when compared to other provisioning models, in order to serve customers with the resources and elasticity they need, clouds generally rely on large data centres that consume massive amounts of electrical power [START_REF] Baliga | Green cloud computing: Balancing energy in processing, storage, and transport[END_REF] [START_REF] Greenberg | The cost of a cloud: Research problems in data center networks[END_REF].

Although some data centres now approach a Power Usage Effectiveness (PUE) factor of 1.08 1 , such a mark means that the IT infrastructure is now responsible for a large part of the consumed power. Means to monitor and analyse how energy is spent are crucial to further improvement, but our previous work in this area has demonstrated that monitoring the power consumed by large systems is not always an easy task [START_REF] Orgerie | Save watts in your Grid: Green strategies for energy-aware framework in large scale distributed systems[END_REF]- [START_REF] Costa | Multifacet approach to reduce energy consumption in clouds and grids: The GREEN-NET framework[END_REF]. There are multiple power probes available in the market, generally with their own APIs, physical connections, precision, and communication protocols [START_REF] Diouri | Solving some Mysteries in Power Monitoring of Servers: Take Care of your Wattmeters![END_REF]. Moreover, cost related constraints can lead data centre operators to acquire and deploy equipments at multiple stages, or to monitor the power consumption of only part of an infrastructure.

From a cost perspective, monitoring the power consumption of only a small part of deployed equipments is sound, but it prevents one from capturing important nuances of the infrastructure. Previous work has shown that as a computer cluster ages, certain components wear out, while others are replaced, leading to heterogeneous power consumption among nodes that were seemingly homogeneous [START_REF] Dean | Some potential areas for future research[END_REF]. The difference between nodes that consume the least power and nodes that consume the most can reach 20% [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF], which reinforces the idea that monitoring the consumption of all equipments is required for exploring further improvement in energy efficiency and evaluate the impact of system-wide policies. Monitoring a great number of nodes, however, requires the design of an efficient infrastructure for collecting and processing the power consumption data. This paper describes the design and architecture of a generic and flexible framework, termed as KiloWatt API (KWAPI), that interfaces with OpenStack to provide it with power consumption information collected from multiple heterogeneous probes. OpenStack is a project that aims to provide ubiquitous open source cloud computing platform and is currently used by many corporations, researchers and global data centres 2 . We describe how KWAPI is integrated into Ceilometer; OpenStack's component conceived to provide a framework to collect a large range of metrics for metering purposes 3 . With the increasing use of Ceilometer as the de facto metering tool for OpenStack, we believe that such an integration of a power monitoring framework into OpenStack can be of great value to the research community and practitioners.

The remaining part of this paper is organised as follows. Section II describes background and related work, whereas Section III presents the KWAPI architecture. Section IV discusses experimental results on measuring the throughput of KWAPI, and Section V concludes the paper.

II. BACKGROUND AND RELATED WORK

This section provides an overview of Ceilometer's architecture and describes related work on monitoring power consumption of large-scale computing infrastructure.

A. OpenStack Ceilometer

Ceilometer -whose logical architecture4 is depicted in Figure 1 -is OpenStack's framework for collecting performance metrics and information on resource consumption. As of writing, it allows for data collection under three methods:

• Bus listener agent, which picks events on Open-Stack's notification bus and turns them into Ceilometer samples (e.g. cumulative type, gauge or delta) that can then be stored into the database or provided to an external system via publishing pipeline.

• Push agents, more intrusive, consist in deploying agents on the monitored nodes to push data remotely to be taken by the collector.

• Polling agents that poll APIs or other tools to collect information about monitored resources. The last two methods depend on a combination of central agent, computer agents and collector. The compute agents run on nodes and retrieve information about resource usage related to a given virtual machine instance and a resource owner. The central agent, on the other hand, executes pollsters on the management server to retrieve data that is not linked to a particular instance. Pollsters are software components executed, for example, to poll resources by using an API or other methods. The Ceilometer database, which can be queried via Ceilometer API, allows an external system to view the history of a resource's metrics, and enables the system to set and receive alarms.

The hmac module of Python's library can be used for signing metering messages, and a shared secret value can be provided in the configuration settings. The collector and systems accessing the API use signatures included in the messages for verification.

B. Energy Monitoring and Efficiency in Clouds

Over the past years, several techniques have been provided to minimise the energy consumed by computing infrastructure. At the hardware level, for instance, processors are able to operate at multiple frequency and voltage levels, and the operating systems or resource managers can choose the level that matches the current workload [START_REF] Laszewski | Poweraware scheduling of virtual machines in DVFS-enabled clusters[END_REF]. At the resource management level, several approaches are proposed, including resource consolidation [START_REF] Beloglazov | Openstack neat: A framework for dynamic and energy-efficient consolidation of virtual machines in openstack clouds[END_REF] and rescheduling requests [START_REF] Orgerie | Save watts in your Grid: Green strategies for energy-aware framework in large scale distributed systems[END_REF], generally with the goal of switching off unused resources or setting them to low power consumption modes. Attempts have also been made to assess the power consumed by individual applications [START_REF] Noureddine | Towards a better understanding of the energy consumption of software systems[END_REF].

A means to monitor the energy consumption is key to assess potential gains of techniques to improve software and cloud resource management systems. Cloud monitoring is not a new topic [START_REF] Aceto | Cloud monitoring: A survey[END_REF] as tools to monitor computing infrastructure [START_REF] Brinkmann | Scalable monitoring system for clouds[END_REF], [START_REF] Varrette | HPC performance and energy-efficiency of the openstack cloud middlewares[END_REF] as well as ways to address some of the usual issues of management systems have been introduced [START_REF] Ward | Varanus: In situ monitoring for large scale cloud systems[END_REF], [START_REF] Tan | Resilient self-compressive monitoring for large-scale hosting infrastructures[END_REF]. Moreover, several systems for measuring the power consumed by compute clusters have been described in the literature [START_REF] De Assunc ¸ão | The green Grid'5000: Instrumenting and using a Grid with energy sensors[END_REF]. As traditional system and network monitoring techniques lack the capability to interface with wattmeters, most approaches for measuring energy consumption have been tailored to the needs of projects for which they were conceived.

In our work, we draw lessons from previous approaches to monitor and analyse energy consumption of large-scale distributed systems [START_REF] Orgerie | Save watts in your Grid: Green strategies for energy-aware framework in large scale distributed systems[END_REF]- [START_REF] Costa | Multifacet approach to reduce energy consumption in clouds and grids: The GREEN-NET framework[END_REF], [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF], [START_REF] Carpentier | CompatibleOne: Designing an energy efficient open source cloud broker[END_REF]. We opt for creating a framework and integrating it with a successful cloud platform (i.e. OpenStack), which we believe is of value to the research community and practitioners working on the topic. To the best of our knowledge, this is the first generic energy monitoring framework to be integrated with OpenStack.

III. THE KWAPI ARCHITECTURE

An overview of the KWAPI architecture is presented in Figure 2. The architecture follows a publish/subscribe model based on a set of layers comprising:

• Drivers, considered data producers responsible for measuring the power consumption of monitored resources and providing the collected data to consumers via a communication bus; and

• Data Consumers -or Consumers for short -that subscribe to receive and process the measurement information.

The communication between layers is handled by a bus, as explained in detail later. Data consumers can subscribe to receive information collected by drivers from multiple sites. Both drivers and consumers are easily extensible to support, respectively, several types of wattmeters and provide additional data processing services. A REST API is designed as a data consumer to provide a programming interface for developers and system administrators. In this work it is used to interface with OpenStack by providing the information (i.e. by polling monitored devices) required by a KWAPI Pollster that feeds Ceilometer.

The following sections provide more details on the main architecture components and their relationship with OpenStack Ceilometer.

A. Driver Layer

Drivers are threads initialised by a Driver Manager with a set of parameters loaded from a file compliant with the Open-Stack configuration format. These parameters are used to query the meters (e.g. IP address and port) and determine the sensor ID to be used in the collected metrics. The measurements that a driver obtains are represented as JavaScript Object Notation (JSON) dictionaries that maintain a small footprint and that

{

"probe_id": "164", "w": 92, "a": 0.4, "v": 230.0, "message_signature": "bb45c88f3ba64..." } Fig. 3. Example of JSON payload. can be easily parsed. The size of dictionaries varies depending on the number of fields set by drivers (i.e. whether message signing is enabled).

Figure 3 shows a simple example of a JSON payload containing one measurement. Optional fields such as voltage and current can be included. ACK messages have a fixed size of 66 bytes when using TCP connection; drivers and data consumers communicate via IPC sockets when running on the same machine.

Drivers can manage incidents themselves, but the manager also checks periodically if all threads are active, restarting them if necessary. It is important to avoid losing measurements because the reported information is in W instead of kWh. The loss of a measurement may be significant.

Wattmeters available in the market vary in terms of physical interconnection, communication protocols, packaging and precision of measurements they take. They are mostly packaged in multiple outlet power strips called Power Distribution Units (PDUs) or enclosure PDUs (ePDUs), and more recently in the Intelligent Platform Management Interface (IPMI) cards embedded in the computers themselves. Support for several types of wattmeter has been implemented, which drivers can use to interface with a wide range of equipments. In our work, we used IPMI initially at Nova to shutdown and turn on compute nodes, but nowadays we also use it to query a computer chassis remotely.

Although Ethernet is generally used to transport IPMI or SNMP packets over IP, USB and RS-232 serial links are also common. Wattmeters that use Ethernet are generally connected to an administration network (isolated from the data centre main data network). Moreover, wattmeters may differ in the manner they operate; some equipments send measurements to a management node on a regular basis (push mode), whereas others respond to queries (pull mode). Other characteristics that differ across wattmeters include:

• refresh rate (i.e. maximum number of measurements per second);

• measurement precision; and

• methodology applied to each measurement (e.g. mean of several measurements, instantaneous values, and exponential moving averages).

Table I shows the characteristics of equipments we deployed and used with Kwapi in our cloud infrastructure.

B. Data Consumers

A data consumer retrieves and processes measurements taken by drivers and provided via bus. Consumers expose the information to other services including Ceilometer and visualisation tools. By using a system of prefixes, consumers can subscribe to all producers or a subset of them. When receiving a message, a consumer verifies the signature, extracts the content and processes the data. By default KWAPI provides two data consumers, namely the REST API (used to interface with Ceilometer) and a visualisation consumer.

1) REST API:

The API consumer computes the number of kWh of each driver probe, adds a timestamp, and stores the last value in watts. If a driver has not provided measurements for a long time, the corresponding data is removed. The REST API allows an external system to retrieve the name of probes, measurements in W or kWh, and timestamps. The API is secured by OpenStack Keystone tokens 5 , whereby the consumer needs to ensure the validity of a token before sending a response to the system.

2) Visualisation:

The visualisation consumer builds Round-Robin Database (RRD) files from received measurements, and generates graphs that show the energy consumption over a given period, with additional information such as average electricity consumption, minimum and maximum watt values, last value, total energy and cost in Euros. RRD files are of fixed size and store several collections of metrics with different granularities. A web interface displays the generated graphics and a cache mechanism triggers the creation of graphs during queries only if they are out of date. These visualisation resources offer quick feedback to administrators and users during execution of tasks and applications. Figure 4 shows an example of generated graph. [START_REF] Hintjens | ZeroMQ: Messaging for Many Applications[END_REF], a fast broker-less messaging framework written in C++, where transmitters play the role of buffers. ZeroMQ supports a wide range of bus modes, including cross-thread communication, IPC, and TCP. Switching from one mode to another is straightforward. ZeroMQ also provides several design patterns such as publish/subscribe and request/response. As mentioned earlier, in our publish/subscribe architecture drivers are publishers, and data consumers are subscribers. If no data consumer is subscribed to receive data from a given driver, the latter will not send any information through the network.

C. Internal Communication Bus

KWAPI uses ZeroMQ

Moreover, one or more optional forwarders can be installed between drivers and data consumers to minimise network usage. Forwarders are designed to act as especial data consumers who subscribe to receive information from a driver and multicast it to all normal data consumers subscribed to receive the information. Forwarders enable the design of complex topologies and optimisation of network usage when handling data from multiple sites. They can also be used to bypass network isolation problems and perform load balancing.

D. Interface with Ceilometer

We opted for integrating KWAPI with an existing open source cloud platform to ease deployment and use. Leveraging the capabilities offered by OpenStack can help in the adoption of a monitoring system and reduce its learning curve.

Ceilometer's central agent and a dedicated pollster (i.e. KWAPI Pollster) are used to publish and store energy metrics into Ceilometer's database. They query the REST API data consumer and publish cumulative (kWh) and gauge (W) counters that are not associated with a particular tenant, since a server can host multiple clients simultaneously.

Depending on the number of monitored devices and the frequency at which measurements are taken, wattmeters can generate a large amount of data, thus demanding storage capacity for further processing and analysis. Management systems often store and perform pre-processing locally on monitored nodes, but such an approach can impact on CPU utilisation and influence the power consumption. In addition, resource managers may switch off idle nodes or set them to stand by mode to save energy, which make them unavailable for processing. Centralised storage, on the other hand, allows for faster data access and processing, but can generate more traffic given that measurements need to be continuously transferred over the network to a central point.

Ceilometer uses its own central database, which is leveraged here to store the energy consumption metrics. In this way, systems that interface with OpenStack's Ceilometer, including Nova, can easily retrieve the data. It is important to notice that, even though Ceilometer provides the notion of a central repository for metrics, it also uses a database abstraction that enables the use of distributed systems such as Apache Hadoop HDFS 6 , Apache Cassandra7 , and MongoDB 8 .

The granularity at which measurements are taken and metrics are computed is another important factor because user needs vary depending on what they wish to evaluate. Taking one or more measurements per second is not common under certain scenarios, which can be a challenge in an infrastructure comprising hundreds or thousands of nodes, demanding efficient and scalable mechanisms for transferring information on power consumption. Hence, in the next section we evaluate the throughput of KWAPI under a few scenarios.

IV. PERFORMANCE EVALUATION

This section provides results of a performance evaluation carried out in our testbed. The goal is not to compare publish/subscribe systems since such work has already been performed elsewhere [START_REF] Eugster | The many faces of publish/subscribe[END_REF], [START_REF] Fabret | Filtering algorithms and implementation for very fast publish/subscribe systems[END_REF]. The evaluation demonstrates that the framework serves well the needs of a large range of users of the Grid'5000 platform [START_REF] Bolze | Grid'5000: a large scale and highly reconfigurable experimental Grid testbed[END_REF] -the infrastructure we use and where the KWAPI framework is currently deployed in production mode as the means for collecting and providing energy consumption information to users.

First we want to evaluate the CPU and network usage of a typical driver to observe the framework's throughput, since provisioning a large number of resources for monitoring purposes is not desirable. For this experiment we deployed the KWAPI drivers and API on a machine with a Core 2 Duo P8770 2.53Ghz processor and 4GB of RAM. We considered:

• a scenario where we emulated 1,000 IPMI cards, each card monitored by a driver thread placing a measurement per second on the communication bus.

• a case with 100 ten-outlet PDUs, each monitored by a driver thread placing ten values per second on the communication bus.

Under both scenarios, 1,000 measurements per second were placed on the bus, even though monitoring was done using different types of probes. We have evaluated these scenarios considering both with and without message signature. Table II summarises the considered scenarios.

Figure 5 shows the results of CPU usage of drivers under the evaluated scenarios. The socket type and number of driver threads do not have a distinguishable impact on the CPU usage. On the test machine, the KWAPI drivers with message signature disabled (i.e. IPMI cards unsigned and PDUs unsigned) consumed on average 20% of the total CPU power.

We also evaluated the CPU consumption of the REST API data consumer under the scenarios described in Table II. In addition to these scenarios, two conditions were assessed, namely (i) the REST API working as a consumer requesting data from drivers at a one-second time interval (REST API only); and (ii) the API requesting data at one-second interval and also answering a call every second to provide the collected data to an external system (REST API + 1 req/s). Figure 6 summarises the obtained results. The CPU consumption is in general low. Even when message signing is enabled and the API serves a query, its consumption is below 20%. The small variation between the scenarios without message signing is caused by the manner ZeroMQ accumulates data on nodes prior to transmission. Although the CPU usage often depends on the drivers, data consumers, and their complexity, and whether message signature is enabled, the experiments show that a large number of probes can be managed by a single machine. In our environment, a management machine per site is more than enough to accommodate the users' monitoring needs. The drivers and API can reuse a machine that already serves other monitoring purposes.

T im e (s e c .) Although a measurement interval of one second meets the requirements of users in our platform, we wanted to evaluate the impact of using a communication bus in the transfer of observations between drivers and the REST API consumer. In a second experiment we used two machines. On the first machine we instantiated 1,000 driver threads placing random observations on the communication bus. On the second machine we measured the number of measurements that the API is able to receive over a minute. We varied the time between measurements from 0.2 to 1.0 seconds. Figure 7 summarises the obtained results. Though the number of observations generated in this experiment is much higher than what we currently need to handle in our platform, we observe that the framework is able to transfer measurements from drivers to API under a 0.4 second interval without adding much jitter. Under smaller measurement intervals, however, observations start to accumulate and are transferred at large chunks. We believe that under small measurement intervals, and consequently a very large number of observations per second, an architecture based on stream processing systems that guarantees data processing might be more appropriate. Hence, although the framework suits the purposes of large range of users, if measurements are to be taken at very small time intervals, a stream processing architecture would probably yield better performance by enabling the placement of elements to pre-process data closer to where it is generated.

V. CONCLUSION

In this paper, we described a novel framework (KWAPI) for monitoring the power consumed by resources of an Openstack cloud. Based on lessons learned by monitoring the power consumption of large distributed infrastructure, we proposed an energy monitoring architecture based on a publish/subscribe model. The framework works in tandem with OpenStack's Ceilometer. Experimental results demonstrate that the overhead posed by the monitoring framework is small, allowing us to serve the users' monitoring needs of our large scale infrastructure.

As future work, we intend to explore means to increase the monitoring granularity and the number of measured devices by applying a hierarchy of plug-ins, and a stream processing system with guarantees on data processing [START_REF]Storm: Distributed and fault-tolerant realtime computation[END_REF], [START_REF]S4: Distributed stream computing platform[END_REF] for processing streams of measurement tuples.

Fig. 1 .

 1 Fig. 1. Overview of Ceilometer's logical architecture.

Fig. 2 .

 2 Fig. 2. Overview of KWAPI's architecture.

Fig. 4 .

 4 Fig. 4. Example of graph generated by the visualisation plug-in (4 monitored servers).

Fig. 5 .

 5 Fig. 5. Driver CPU usage under the evaluated scenarios.

Fig. 6 .

 6 Fig. 6. API consumer CPU usage under the evaluated scenarios.

Fig. 7 .

 7 Fig. 7. Number of observations received over a 60 second interval under different multiple intervals.

TABLE I .

 I

		WATTMETER INFRASTRUCTURE
	Device Name	Interface	Refresh Time (s)	Precision (W)
	Dell iDrac6	IPMI / Ethernet	5	7
	Eaton	Serial, SNMP via Ethernet	5	1
	OmegaWatt	IrDA Serial	1	0.125
	Schleifenbauer	SNMP via Ethernet	3	0.1
	Watts Up?	Proprietary via USB	1	0.1
	ZEZ LMG450	Serial	0.05	0.01

TABLE II .

 II SCENARIOS CONSIDERED IN THE EXPERIMENTS.

	Scenario name		Agent thread scheme	Message signature
	IPMI message signed	1 thread per card		Enabled
	IPMI message unsigned	1 thread per card		Disabled
	PDU message signed	1 thread per PDU		Enabled
	PDU message unsigned	1 thread per PDU		Disabled
	PDUs unsigned
	IPMI cards unsigned
	PDUs signed
	IPMI cards signed
	0	5	10	15 CPU Usage (%)	20	25	30

http://gigaom.com/2012/03/26/whose-data-centers-are-more-efficientfacebooks-or-googles/

http://www.openstack.org/user-stories/

https://wiki.openstack.org/wiki/Ceilometer

http://docs.openstack.org/developer/ceilometer/architecture.html

http://keystone.openstack.org

http://hadoop.apache.org/

http://cassandra.apache.org/

http://www.mongodb.org/

ACKNOWLEDGMENTS

This research is supported by the French Fonds national pour la Société Numérique (FSN) XLCloud project. Some experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the Inria ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://grid5000.fr). Authors wish to thank Julien Danjou for his help during the integration of KWAPI with OpenStack and Ceilometer.