D Chapelle

M Fragu

V Mallet

P Moireau

Fundamental principles of data assimilation underlying the Verdandi library: applications to biophysical model personalization within euHeart

We present the fundamental principles of data assimilation underlying the Verdandi library, and how they are articulated with the modular architecture of the library. This translates -in particular -into the definition of standardized interfaces through which the data assimilation library interoperates with the model simulation software and the so-called observation manager. We also survey various examples of data assimilation applied to the personalization of biophysical models, in particular for cardiac modeling applications within the euHeart European project. This illustrates the power of data assimilation concepts in such novel applications, with tremendous potential in clinical diagnosis assistance.

Introduction

In order to obtain some information -as detailed as possible -on a natural system, such as in geophysics, or regarding the important subcategory of living systems which constitute the objects of study in biology and medicine, the most straightforward strategy consists in obtaining measurements on the system at hand. Note that we deliberately use the term measurement, related but not reduced to experiments, to signify that it is not in general possible to design specific experiments allowing to determine all kinds of physical properties of the system -as is done for instance for industrial systems. By contrast, natural systems induce drastic limitations in measurements, in that they generally must be "taken as they are", namely, observed in their current operating conditions, whether this is due to the practical impossibility of apprehending them comprehensively (e.g. in geophysics), or to the undesirable character of any strong perturbation of the system (invasiveness in living systems).

Specializing now our discussion to the human body, abundant measurements are frequently at hand, in the form of clinical images and signals of various origins. Despite the diversity and rich information contents of these data, they are also inevitably limited in many respects. Beside considerations on sampling and noise, this holds in particular as regards their extent. For example, only 2D measurements, or boundary information, may be available for a 3D system, or only part of the whole domain. Limitations also pertain to measurement types, as some quantities are never measured, such as internal stresses in a living tissue, or various physical constitutive parameters (stiffness, contractility, etc.).

Nevertheless, we may want to consider mathematical models to describe and predict the behavior of these systems. Clearly, these models may be seen as providing some complementary information on the system. However, their predictivity requires the careful adjustment of many parameters -in particular regarding the detailed geometry (anatomy), physical properties, boundary conditions and initial conditions needed in the model -most of which being out of reach of the available measurements.

The purpose of data assimilation is then to combine the information available from these two sources -measurements on the one hand, and models and the other hand -by seeking an adequate compromise between, on the one hand, the discrepancy computed between simulations of the model and the corresponding measurements and, on the other hand, the a priori confidence in the model, since errors are also present in the measurements. The desired output of this procedure is an estimation of the unknown quantities of interest, namely, (1) state variables (the "trajectories" of the system), and in particular their initial values at a given reference time, and (2) physical parameters which must be prescribed in the model equations.

In terms of clinical applications, the expected benefits are to assist and improve both diagnosis and prognosis. Diagnosis can be enhanced by providing more complete information on the patient (spatially-distributed quantities, and various otherwise unreachable indicators), and with improved accuracy. Benefits also extend to prognosis, since once data assimilation has been performed the model can be more confidently used to predict natural or artificial evolutions of the system, for example to simulate the effect of various possible therapeutic strategies. Note that, in a prognosis perspective, depending on the application at hand, modeling time scales can vary from -typically -microseconds at the sub-cellular level to seconds (e.g. a heart beat), to hours or weeks when longer evolutions are considered, for instance tissue remodeling.

In this paper, we present fundamental concepts of data assimilation, for both variational and sequential families of methods. We also survey some application examples pertaining to the personalization of biophysical models, in particular for cardiac modeling applications within the euHeart European project. This leads us to discussing the architecture of the Verdandi library, and how it relates to the fundamentals of data assimilation and addresses applicative needs.

Definitions and notation

We now introduce some basic notation necessary to discuss the fundamental principles of data assimilation.

Physical model definition -First of all, we consider physical models in the form of dynamical systems governed by equations of the type ẋ = A(x, θ, t).

(

In this equation, x denotes the so-called state variable, namely, the physical quantity which the model aims at describing in its time-wise evolution -hence, the time derivative in the left-hand side -and also frequently spatial variations for distributed quantities. In this generic notation, the whole model is essentially summarized in the so-called dynamical operator A, which applies on the state variable itself, and may depend on time t as well as on a set of physical parameters denoted by θ. This operator may arise from various types of physical formulations, e.g. in solid and fluid mechanics or electrophysiology. Mathematically speaking, it may take the form of partial differential equations (PDEs) or ordinary differential equations (ODEs, namely, only differentiated with respect to the time variable), or algebraic systems, in particular.

Clearly, in such model formulations we need to prescribe -hence to estimate when unknown via the data assimilation procedure -the initial condition x(0) and the parameter vector θ. We dissociate the estimation of the state through x(0) and the estimation of the parameter θ (also called identification), but formally the two types of estimation can be considered together by defining an augmented state x = ( x θ ) verifying ẋ = A(x), and x(0

) = x(0) θ , (2)
with a slight abuse of notation since the new augmented operator A has two components A(•) 0 . Then every estimation can be considered in the light of state estimation.

The variables x and θ can represent fields in a space-distributed form, but they are then space-discretized. In that case we will denote by

˙ X = A( X, θ, t), (3)
the space discretization of (1). Typically, in the models considered the state variable may contain a large number of scalar coefficients -typically 10 3 to 10 7 degrees of freedom in a continuum mechanics model -whereas the size of the parameter vector is generally much more limited, and in practice we seldom have to estimate more than a few hundreds of parameter values. Once the initial condition and the parameter vector associated with [START_REF] Auroux | A nudging-based data assimilation method: the Back and Forth Nudging (BFN) algorithm[END_REF] are estimated, the model can be simulated in time, using appropriate numerical techniques. We also have to consider time discretizations of the model for solving the dynamics in practice. Every time scheme can be summarized in the form

X n+1 = A n+1|n ( X n , θ), (4)
where we frequently have X n = X(t n ), even though in some specific discretization schemes we may gather more variables. Finally, let us point out that every effort to model a physical system suffers from modeling errors. In the data assimilation formalism the modeling error can be considered as a timedependent quantity ω(t) that also needs to be estimated. Then we have ẋ = A(x, θ, ω, t).

(

) 5
with, most of the time, a linear expression of the form

ẋ = A(x, θ, t) + Bω(t). (6)
This quantity can be considered in a deterministic framework as a perturbation variable in a specific space, or as a probabilistic uncertainty.

Measurements description -Another important notation concerns the measurements, typically represented by an equation of the type

z = H(x, t) + χ(t), (7)
where z denotes the actual data, H is the so-called observation operator, and χ accounts for the error inherent to the measurement process, often called the noise, and for representativeness errors. Note that the quantity z will frequently correspond to pre-processed -not raw -data, for example medical images processed with segmentation or optical flow techniques in order to extract some position, displacement or velocity information. We further emphasize that this equation also represents a model -in this case of the measurements -where modeling ingredients are embedded both in the expression of H and the characterization of χ, which may be of probabilistic or deterministic nature.

In practice, we have to take into account that the measurements are also discretized in space and time with a specific level of space discretization and time sampling. It is still possible to rewrite [START_REF] Delingette | Personalization of Cardiac Motion and Contractility From Images Using Variational Data Assimilation[END_REF] to take into account these discretizations and to adapt the observation operator to model discretization choices. The new observation definition reads

Z n = H n ( X n ) + χ n,∆T , (8)
where in χ n,∆T we gather the contributions coming from the measurement noise and the discretization errors (possibly deterministic and biased). For example, for data available at a constant time-sampling ∆T -in general different from, and frequently larger than, the computational time step of the model simulation -we can consider that the observations are only available at simulation time steps coinciding with (or closest to) data sampling times, i.e.

Z n = Z(k∆T ) if |t n -k∆T | ≤ obs ∆T undefined otherwise
Alternatively, we can regenerate an observation at every simulation time step by considering some time interpolation strategy, for instance a linear scheme, viz.

Z n = (k+1)∆T -tn ∆T Z(k∆T ) + tn-k∆T ∆T Z((k + 1)∆T ) , for t n ∈ k∆T, (k + 1)∆T .

Fundamental principles and methods

We now proceed to present some fundamental concepts of data assimilation. For the sake of simplicity, we mostly restrict our presentation to a deterministic point of view, whereas many concepts can also be interpreted in a probabilistic light, namely, within the so-called Bayesian framework, see e.g. [START_REF] Bensoussan | Filtrage Optimal des Systèmes Linéaires[END_REF] and references therein.

Least square minimization: the BLUE algorithm

Let us start by considering only one time and one observation, and let us also assume that the observation operator is linear. The most straightforward approach consists in minimizing the following criterion min

x J (x) = 1 2 z -Hx 2 M + 1 2 x -x 2 N , (9)
where M and N denote the norms associated with the observation space and state space, respectively. In other words, M and N characterize our understanding of the underlying quantities z and x -for example their regularities if they correspond to distributed fields -and they also weigh our confidence in each term in the global balance expressed by the criterion. When considering probabilistic fields, the inverse of these norm operators can be associated with the measurement noise covariance W for the observations, and the state covariance P for the state when also choosing x = E(x) as the expected (or mean) value of x.

In practice this minimization is solved after adequate spatial discretization -when x represents a field -leading to a criterion

J( X) = 1 2 Z -H X 2 M + 1 2 X -X 2 N = 1 2 ( Z -H X) W -1 ( Z -H X) + ( X -X ) P -1 ( X -X ), (10)
where M and N denote the matrix counterparts of the above norms, and W and P the respective inverse matrices, namely, covariance matrices. It is easy to show that the result of this minimization is

ˆ X = X + (P -1 + H W -1 H) -1 H W -1 ( Z -H X ). (11)
By defining

P + = (P -1 + H W -1 H) -1 , (12)
we see that [START_REF] Imperiale | Constitutive parameter estimation methodology using tagged-MRI data[END_REF] becomes

ˆ X = X + P + H W -1 ( Z -H X ). ( 13
)
Using the matrix inversion lemma (see e.g. [START_REF] Moireau | Filtering-based Data Assimilation for Second-Order Hyperbolic PDEs[END_REF]), we obtain that ( 11) can be rewritten as

ˆ X = X + P H (W + HP H ) -1 ( Z -H X ). ( 14
)
These formulations provide a linear least square estimator ˆ X of the state X given one observation Z and an a priori X . Everything can be summarized in the filtering form

       ˆ X = X + K( Z -H X ) with K = P + H W -1 = (P -1 + H W -1 H) -1 H W -1 = P H (W + HP H ) -1 (15)
This estimator is classically referred to as the BLUE (best linear unbiased estimator) method.

In the case of a non-linear observation operator H(•), it is still possible to compute the minimum via e.g. a Newton-like algorithm. We obtain

       ˆ X (k+1) = X (k) + K (k) ( Z -H( X (k) )) with K (k) = P + (k) ∂H (k) ∂ X W -1 = P -1 (k) + ∂H (k) ∂ X W -1 ∂H (k) ∂ X -1 ∂H (k) ∂ X W -1
but we can also be satisfied with an approximate version using only the first iteration of the Newton algorithm. In this case, it is sufficient to use [START_REF] Moireau | Filtering-based Data Assimilation for Second-Order Hyperbolic PDEs[END_REF] with H substituted with ∂H ∂ X in the definition of the gain K.

BLUE with multiple observations

Now, let us explain how this can be extended when considering multiple observations ( Z n ). Using minimization principles we would find that the estimator is given by

ˆ X n = n k=0 H k W -1 k H k -1 n k=0 H k W -1 k Z k . (16)
Then, defining the sequence of symmetric matrices (P n )

P -1 n = n k=0 H k W -1 k H k , (17)
we can show as in [START_REF] Moireau | Filtering-based Data Assimilation for Second-Order Hyperbolic PDEs[END_REF] that ˆ X n can be computed recursively by

ˆ X n = ˆ X n-1 + K n ( Z n -H n ˆ X n-1 ), (18)
with the filter given by

K n = P n-1 H n (W n + H n P n-1 H n ) -1 = (H n W -1 n H n + P -1 n-1 ) -1 H n W -1 n = P n H n W -1 n . (19)
This recursive form can only be defined in a linear context but an approximate version can then be formulated for nonlinear operators by again substituting H n with ∂Hn ∂ X in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF].

General variational minimization: The 4D-Var

The previous section allowed us to understand in a simple case the equivalence between criterion minimization and recursive estimation formulae, and also the link between the probabilistic point of view -involving mean and covariances computations -and the deterministic point of view. We can now consider the more general configuration with a model given as a dynamical system and a time-sequence of observations. We demonstrate the detailed strategy for the time-continuous case, but similar principles also apply to time-discrete systems that can reflect the time-discretization of an underlying time-continuous model.

In variational procedures, we consider a criterion to be minimized in order to achieve the above-mentioned compromise between simulation-measurements discrepancy and model confidence, see e.g. [START_REF] Bensoussan | Filtrage Optimal des Systèmes Linéaires[END_REF][START_REF] Chavent | Nonlinear Least Squares for Inverse Problems[END_REF] and references therein. A typical criterion would read

J T (ξ x , ξ θ ) = T 0 z -H(x) 2 M dt + ξ x 2 (P ) -1 + ξ θ 2 (P * ) -1 , (20)
where . M , . (P ) -1 and . (P * ) -1 denote suitable norms for each quantity concerned, and associated with the operators appearing as subscripts. In this criterion, x(t) is constrained to satisfy the model equation ( 1) starting from the initial condition x(0) = x 0 + ξ x and with parameter values given by θ = θ 0 + ξ θ . In order to perform this minimization, a classical strategy consists in computing the gradient of the criterion, which requires the simulation of the so-called adjoint model. The adjoint model equation is an evolution system closely related to -and inferred from, indeed -the direct model equation ( 1), viz.

           ṗx + ∂A ∂x p x = -∂H ∂x M (z -H(x)) p x (T ) = 0 ṗθ + ∂A ∂θ p x = 0 p θ (T ) = 0 (21)
These equations must be simulated backwards in time from the final time T to the initial time in order to obtain the gradient value criterion expressed as

d ξx J • δξ x = ξ x (P ) -1 δξ x -p x (0) δξ x d ξ θ J • δθ = ξ θ (P * ) -1 δξ θ -p θ (0) δξ θ ,
Hence, each gradient computation requires the forward simulation of the direct model and the backward simulation of the adjoint, and this must be repeated until convergence of the minimization algorithm. This type of variational procedure is also referred to as "4D-Var" in the data assimilation community, while "3D-Var" is used to refer to minimization estimation performed for static models, or for dynamic models at a given time (namely, without time integral).

Concerning the time discretization, it is classical to formulate an optimal discrete time minimization criterion and find its corresponding adjoint rather than discretizing directly [START_REF] Pham | Stochastic methods for sequential data assimilation in strongly nonlinear systems[END_REF]. Hence, we consider a criterion of the form

J N (ξ x , ξ θ ) = N k=0 Z k -H( X k ) 2 M k + ξ x 2 (P ) -1 + ξ θ 2 (P * ) -1 , (22)
with for example M k = ∆tM.

3.4 Sequential optimal procedure: The Kalman algorithms family

By contrast, sequential procedures -also often referred to as filtering -proceed by simulating equations closely resembling the direct model equations, with an additional correction term taking into account the discrepancy between the simulation and the actual measurements, namely z -H(x), quantity called the innovation. For example when only the initial condition is unknown the filtering equation would be of the type

ẋ = A(x, t) + K(z -H(x)), (23)
where the operator K -frequently linear -is called the filter. The filtering equation simulation is then started from the candidate initial condition x , and the aim of the correction is to bring the simulated trajectory close to the target system. This type of strategy was made extremely popular by the Kalman theory, which formulated an optimal setting for deriving the filter operator, initially when the operators A and H are both linear. In this case, the Kalman equations read

           ẋ = Ax + P H M (z -H x) Ṗ -P A -AP + P H M HP = 0 P (0) = P x(0) = x (24)
where W = M -1 and P denote the so-called covariance operators of the measurement noise and initial condition uncertainty, respectively. We can see that the filter expression is based on the computation of the time-dependent covariance P which satisfies a Riccatti equation. In fact, in the linear case the variational and Kalman procedures can be shown to be exactly equivalent, and the minimizing direct and adjoint states (x ∞ and p ∞ , respectively) are related to the Kalman filter equations by the identity [START_REF] Bensoussan | Filtrage Optimal des Systèmes Linéaires[END_REF] x ∞ = x + P p ∞ .

When nonlinearities are to be considered, various extensions are available, and in particular the Extended Kalman Filtering (EKF) approach, in which the linearized forms of the operators are used in the filter equation. However, in such a case the approach is no longer equivalent to the variational setting. Nevertheless, some alternative filter equations can be derived from the variational formulation, but the filter computation then requires solving a Hamilton-Jacobi-Bellman equation in a space which has the dimension of the state variable [START_REF] Moireau | Filtering-based Data Assimilation for Second-Order Hyperbolic PDEs[END_REF], which is in general not practical.

Note that, in practice, the Kalman (or EKF) approach itself is also quite limited as regards the size of the system which can be handled, since the covariance operator P has the size of the state variable, and is "dense", unlike the dynamical operators. In order to circumvent this limitation some alternative approaches can be proposed, see Section 3.5, or discretizations involving much fewer degrees of freedom must be considered. This is particularly the case in the context of reduced basis or POD discretizations [START_REF] Chapelle | A Galerkin strategy with Proper Orthogonal Decomposition for parameter-dependent problems -Analysis, assessments and applications to parameter estimation[END_REF].

Concerning time discretization, it is classical to formulate the optimal filter directly from the optimal time-discrete minimization criteria rather than discretizing directly [START_REF] Smith | euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF]. Hence after some quite tedious computations very similar to those presented in Section 3.2, we can formulate a prediction-correction scheme 1. Prediction:

   ˆ X - n+1 = A n+1|n ( X + n ) P - n+1 = ∂A n+1|n ∂ X P + n ∂A n+1|n ∂ X (25a)
2. Correction:

       P + n+1 = ∂H n+1 ∂ X W -1 n+1 ∂H n+1 ∂ X + (P - n+1 ) -1 -1 K n+1 = P + n+1 ∂H n+1 ∂ X W -1 n+1 ˆ X + n+1 = ˆ X - n+1 + K n+1 ( Z n+1 -H n+1 ( ˆ X - n+1 )) (25b)
For the sake of simplicity, we have introduced the Kalman filter in a deterministic context, but the probabilistic counterpart exists. In a linear framework the expressions are exactly the same, albeit with the additional interpretation

             a priori mean: ˆ X - n+1 = E( X n+1 | Z 0 , . . . , Z n ),
a priori covariance:

P - n+1 = E ( X n+1 -ˆ X - n+1 )( X n+1 - ˆ X - n+1 ) , a posteriori mean: ˆ X + n+1 = E X n+1 | Z 0 , . . . , Z n+1 ,
a posteriori covariance:

P + n+1 = E(( X n+1 -ˆ X + n+1 )( X n+1 - ˆ X + n+1 ) ), (26)
This results extend to the non-linear context with the EKF (25) but the identities (26) are then only approximate. To improve the quality of this approximation, the Unscented Kalman Filter has then been introduced [START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filter and estimators[END_REF], based on the idea of substituting means and covariances by empirical quantities computed from sample points:

             ˆ X - n+1 = d i=1 α i X [i]- n+1 , P - n+1 = d i=1 α i ( X [i]- n+1 - ˆ X - n+1 )( X [i]- n+1 - ˆ X - n+1 ) , ˆ X + n+1 = d i=1 α i X [i]+ n+1 , P + n+1 = d i=1 α i ( X [i]+ n+1 - ˆ X + n+1 )( X [i]+ n+1 - ˆ X + n+1 ) , (27)
with d i=1 α i = 1. In practice, the correction particles are sampled around the mean ˆ X + n with a covariance P + n and the prediction samples then verify

X [i]- n+1 = A( X [i]+ n ). ( 28a
)
Then, by computing

Z [i]- n+1 = H( X [i]- n+1 ), Z - n+1 = d i=1 α i Z [i] n+1 (28b)
The gain is defined by

       K n+1 = (P X, Z n+1 ) • (P Z, Z n+1 ) -1 P X, Z n+1 = d i=1 α i ( X [i]- n+1 - ˆ X - n+1 )( Z [i]- n+1 -Z - n+1 ) P Z, Z n+1 = d i=1 α i ( Z [i]- n+1 -Z - n+1 )( Z [i]- n+1 -Z - n+1 ) + W n+1 (28c)
so that we keep having

ˆ X + n+1 = ˆ X - n+1 + K n+1 ( Z n+1 -Z - n+1 ) P + n+1 = P - n+1 -P X, Z n+1 (P Z, Z n+1 ) -1 (P X, Z n+1 ) (28d)
and proceed recursively with new correction particles

X [i]+
n+1 . The Ensemble Kalman Filter, introduced in [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF], follows the same principles of approximating the covariances by sampled particles with, most of the time, an increased number of particles with respect to the UKF filter, and various ways of sampling the particles around the mean value. Finally Monte Carlo strategies exploit a very large number of particles to give a better approximation of the non-linear optimal filter but the practical details of these methods are beyond the scope of this review focused on large dimensional systems coming from the discretization of PDEs.

Reduced-order sequential strategies

Reduced-Order Extended Kalman Filtering (ROEKF)

In order to deal with the limitations of sequential strategies due to the system size, a classical strategy consists in assuming a specific reduced-order form for the covariance operators. For example, making the ansatz ∀t,

P (t) = L(t)U (t) -1 L(t) (29)
with U an invertible matrix of small size r and L an extension operator, we can show that within linear assumptions the solution of the Riccatti equation in [START_REF] Smith | euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF] reduces to L = AL and U = L H M HL.

which is actually computable in practice.

In a non-linear framework, we can then approximate the covariance dynamics by extending (30) as

L = ∂A ∂x L and U = L ∂H ∂x M ∂H ∂x L. (31)
These strategies have relevant applications in the case of parameter identification per se. It is common, indeed, to assume more space regularity for the parameters than for the state initial condition. Hence, after discretization the parameters can be represented by a small number of degrees of freedom. Assuming that we can limit U (0) to the parametric space, the extension operator can then be decomposed into two components

L = Lx L θ with L θ = 1 and Lx = ∂A ∂x L x + ∂A ∂θ . (32)
We recognize in this expression the dynamics of the sensitivity operator ∂x ∂θ , which provides a nice interpretation for this strategy of uncertainty covariance reduction. Furthermore, in the linear framework we can prove that this sequential estimator corresponds to the optimal filter associated with the criterion

J (ξ θ ) = T 0 z -H(x) 2 M dt + ξ θ 2 (P * ) -1 , (33)
where x follows the trajectory associated with ξ θ and fixed initial condition, which is commonly used in variational identification procedures. All this justifies naming this strategy Reduced-Order Extended Kalman Filter (ROEKF), but it is also known as Singular Evolutive Extended Kalman Filter following the work of [START_REF] Pham | A singular evolutive extended Kalman filter for data assimilation in oceanography[END_REF]. This concept can be applied directly on time and space discretized versions of the equations, which leads to a discrete time Reduced-Order Extended Kalman Filter.

Prediction:

ˆ X - n+1 = A n+1|n ( X + n ) L n+1 = ∂A n+1|n ∂ X L n (34a) 2. Correction:        U n+1 = U n + L n+1 ∂H n+1 ∂ X W -1 n+1 ∂H n+1 ∂ X L n+1 K n+1 = P + n+1 ∂H n+1 ∂ X W -1 n+1 ˆ X + n+1 = ˆ X - n+1 + K n+1 ( Z n+1 -H n+1 ˆ X - n+1 ) (34b)

Reduced-Order Unscented Kalman Filtering (ROUKF)

Alternatively, this strategy can be coupled with the UKF approach by showing that particles can be generated only in the space of small dimension and the computation made in the UKF filter (28) can be compatible with the time discretized counterpart of (29). This was proven in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF][START_REF] Moireau | Erratum of article "reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] which also provided a very general version of the Reduced Order Unscented Kalman Filter (ROUKF). In particular this algorithm is close to the Singular Evolutive Interpolated Kalman Filter [START_REF] Pham | Stochastic methods for sequential data assimilation in strongly nonlinear systems[END_REF][START_REF] Hoteit | A simplified reduced order Kalman filtering and application to altimetric data assimilation in Tropical Pacific[END_REF] for a choice of particles d = r + 1 and reads as follows.

Algorithm -Given an adequate sampling rule, we store the corresponding weights in the diagonal matrix D α and precompute so-called unitary sigma-points (i.e. with zero mean and unit covariance) denoted by ( I [i] ) 1≤i≤r+1 ; we then perform at each time step 1. Sampling:

   C n = (U n ) -1 ˆ X [i]+ n = ˆ X + n + L n • C n • I [i] , 1 ≤ i ≤ r + 1 (35a) 2. Prediction:    ˆ X [i]- n+1 = A n+1|n ( ˆ X [i]+ n ), 1 ≤ i ≤ r + 1 ˆ X - n+1 = r+1 i=1 α i ˆ X [i]- n+1 (35b)
3. Correction:

                           L n+1 = [ ˆ X [ * ]- n+1 ]D α [ I [ * ] ] Z [i]- n+1 = H n+1 ( ˆ X [i]- n+1 ) Z - n+1 = r+1 i=1 α i Z [i]- n+1 Γ n+1 = [ Z [ * ]- n+1 ]D α [ I [ * ] ] U n+1 = 1 + Γ n+1 W -1 n+1 Γ n+1 ˆ X + n+1 = ˆ X - n+1 -L n+1 U n+1 Γ n+1 W -1 n+1 ( Z n+1 -Z - n+1 ) (35c)
where we denote by [ ] ] the matrix concatenating the ( I [i] ) vectors side by side, and similarly for other vectors.

I [ *

Luenberger observers

The so-called observer theory, initiated by Luenberger [START_REF] Luenberger | Determining the State of a Linear with Observers of Low Dynamic Order[END_REF], is based on the simple realization that, defining the estimation error x = x -x, we obtain in the linear case, when subtracting the direct model and filtered equations, the dynamics ẋ = (A -KH)x -Kχ.

This type of dynamical equation is well-known in control theory: it is similar to the closedloop controlled equation of a system of natural dynamics governed by A, and submitted to a feedback control defined by the operator K applied on the quantity observed through H. Hence, obtaining an accurate estimation of the state variable is exactly equivalent to driving the estimation error x to zero -namely, stabilizing this error -using the feedback control K. This approach opened new avenues for formulating novel filtering approaches, because for actual dynamical systems control and stabilization motivations have frequently already led to the formulation of effective feedback controls, used in a large variety of industrial systems. Hence, these approaches can be quite directly adapted to obtain adequate filters which -unlike the Kalman filter -are tractable in practice for large systems. Moreover, these filters are often deeply rooted in the physics of the system considered, hence the computational building blocks needed are likely to be at hand in the system simulation software. However, in the Luenberger approach we lose the Kalman optimality, which of course only holds in quite restricted (linear) cases.

For examples of such approaches applicable to beating heart models, with detailed descriptions, theoretical analyses and extensive assessments, we refer in particular to [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF][START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF]. We also point out that the Luenberger observer approach is also sometimes referred to as "nudging" in the data assimilation community [START_REF] Auroux | A nudging-based data assimilation method: the Back and Forth Nudging (BFN) algorithm[END_REF].

Joint state-parameter estimation with Luenberger observers

As apparent in the above discussion, Luenberger observers were originally designed for state estimation. When parameters are to be jointly estimated, it is quite classical in the filtering context to complement the state equation ( 1) with the artificial parameter dynamics θ = 0, (37) as already discussed when introducing [START_REF] Bensoussan | Filtrage Optimal des Systèmes Linéaires[END_REF]. Then the whole estimation objective is to estimate the initial condition of the augmented state (x, θ). Of course, when a Kalman approach is out of reach for the state variable alone, it holds a fortiori for the augmented state. On the other hand, devising a Luenberger observer for the augmented state is difficult, because part of the dynamics is non-physical, hence feedback controls are not readily available for the joint system. Nevertheless, an effective approach for joint state-parameter estimation was proposed in [START_REF] Moireau | Joint state and parameter estimation for distributed mechanical systems[END_REF], based on a Luenberger observer applied on the state equation alone. In essence, this firststage state estimation reduces the uncertainty to the parameter space, which allows to consider a ROEKF approach for handling the remaining parameter uncertainty. This algorithm can be summarized as

                               ẋ = A(x, θ) + K x (z -Hx) + L x θ θ = U -1 L x H M (z -Hx) Lx = ∂A ∂x -K x H L x + ∂A ∂θ U = L x H M HL x x(0) = x 0 θ(0) = θ 0 L x (0) = 0 U (0) = U 0 ( 38
)
where K x denotes the state filter (Luenberger observer operator), L x represents the sensitivity of the state variable with respect to the parameters, and U the inverse of the parameter estimation error covariance. This methodology derives from the above-discussed reduced-order filtering approaches, see also e.g. [START_REF] Pham | A singular evolutive extended Kalman filter for data assimilation in oceanography[END_REF], since only the part of the dynamics concerning parameters is handled using optimal filtering. This joint estimation approach was later extended in [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] towards strategies inspired from the ROUKF crucially allowing to avoid the computation of differentiated operators required in (38).

Application examples 4.1 Within the euHeart consortium

Two recent results have illustrated the great potential of data assimilation in the context of the biophysical personalization of cardiovascular systems as purported in the euHeart project [START_REF] Smith | euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF].

The first example in [START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reducedorder unscented Kalman filter[END_REF] consists in estimating a parameter associated with the passive part of a cardiac constitutive law. Here, the tissue passive behavior is assumed to be described by the Guccione law, one of the main material laws used in cardiac mechanics. The overall mechanical model is assumed to be quasi-static and submitted to time-dependent external forces. In this context, the strategy followed by [START_REF] Xi | Myocardial transversely isotropic material parameter estimation from in-silico measurements based on a reducedorder unscented Kalman filter[END_REF] is similar to a BLUE filter with multiple observations, since there is no dynamics associated with the state. In this context, the uncertainty was reduced to the parameter space related to the Guccione law, and an ROUKF algorithm was chosen for its combined robustness and ease of implementation. The results on synthetically generated data demonstrate the potential of the approach for the estimation of such complex constitutive laws.

The second example in [START_REF] Bertoglio | Sequential parameter estimation for fluidstructure problems. application to hemodynamics[END_REF] also used a ROUKF filter, albeit in its full version, since the model considered is a complete dynamical fluid-structure interaction problem in large strains. Here again, the motivation is to estimate some parameters in the solid constitutive law designed for the arterial walls. The authors show in particular how a complex time-discretization scheme can be processed in the data assimilation formalism presented in the previous section. In this case also, results on synthetically generated data demonstrate the effectiveness of the approach.

Other examples in the cardiovascular context

The previous two estimation problems have also been considered within variational procedures, in [START_REF] Xi | An automatic data assimilation framework for patient-specific myocardial mechanical parameter estimation[END_REF] for the estimation of a cardiac constitutive law and in [START_REF] D'elia | A variational data assimilation procedure for the incompressible Navier-Stokes equations in hemodynamics[END_REF] for the estimation of wall parameters in blood flows with fluid-structure interaction models. In the second case synthetic data experiments are performed but theoretical stability results are provided whereas in [START_REF] Xi | An automatic data assimilation framework for patient-specific myocardial mechanical parameter estimation[END_REF], the Guccione law parameters are obtained from real data and compared to the literature. Concerning the estimation of the active part of a mechanical model, variational methods have also been employed in [START_REF] Delingette | Personalization of Cardiac Motion and Contractility From Images Using Variational Data Assimilation[END_REF] to identify some contractility parameters. In all these results, the data are assumed to be displacements extracted by adequate processing techniques from the image cine-MRI or tagged-MRI.

Data assimilation is also becoming more popular in electrophysiology, where variational strategies have been employed e.g. in [START_REF] Moreau-Villeger | Building maps of local apparent conductivity of the epicardium with a 2-D electrophysiological model of the heart[END_REF][START_REF] Relan | Coupled Personalization of Cardiac Electrophysiology Models for Prediction of Ischaemic Ventricular Tachycardia[END_REF] or filtering methods in [START_REF] Wang | A reduced-rank square root filtering framework for noninvasive functional imaging of volumetric cardiac electrical activity[END_REF]. The identification of parameters related to the cell electric properties are the main objectives, but state estimation of the propagating wave can also be valuable, in particular when the electrical activity is an input in a mechanical model to initiate the contraction. Therefore, Eikonal approximations of the wave propagation have been combined with probabilistic filtering in [START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to Eikonal-Diffusion models in cardiac electrophysiology[END_REF] to estimate the wave trajectory and some propagation parameters.

Specific examples with image data

Image data provide a very rich source of information, albeit are frequently difficult to directly associate with a measurement operator H, e.g. providing some displacements corresponding to a biomechanical model in a Lagrangian formulation. In general, it is much easier to extract from the cardiac or cardiovascular images only some surfaces of interest describing the contours of the tissue -or the tagged planes in tagged MRI. This kinematical information is fundamentally Eulerian, as opposed to Lagrangian displacements in the model. Therefore, in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF] and [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF] for domain contours or in [START_REF] Imperiale | Constitutive parameter estimation methodology using tagged-MRI data[END_REF] for tagged data, the authors proposed to extend the classical data assimilation observation operator to a so-called discrepancy operator, in order to directly use the information provided by the segmentation processes. The key idea is to note that the data and observation operator are always considered together in data assimilation methods in the form z -H(x). As initiated in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspectives in medical imaging[END_REF], this quantity is thus replaced by a non-linear form by considering an operator D, z -H(x) → D(x, t),

In the case of segmentations and model comparison we can then consider a discrepancy operator based on distances computations

D(x, t) = dist(x(ξ, t), t),
with, e.g., interpolation rules between successive images dist :

L 2 (Σ) × [0, T ] → L 2 (Σ) (x(ξ), t) → dist(x(ξ), t) = t k+1 -t ∆T dist S k (x(ξ)) + t-t k ∆T dist S k+1 (x(ξ)) × t k+1 -t ∆T n S k (x(ξ)) + t-t k ∆T n S k+1 (x(ξ)) , if t ∈ [t k , t k+1 ] (39)
This methodology was applied in particular in [START_REF] Chabiniok | Estimation of tissue contractility from cardiac cine-mri using a biomechanical heart model[END_REF] to estimate regional contractility parameters in an infarcted heart based on actual MR images, thereby demonstrating the potential of the complete modeling and estimation chain for diagnosis assistance purposes.

Verdandi: a generic data assimilation library

There exists data assimilation software with different scientific contents and technical designs. For reference, we list the main software that is adapted to high-dimensional applications that we are aware of. The OpenDA library (http://www.openda.org) focuses on Kalman filters and methods for parameter estimation. It is implemented in Java and relies on XML configuration files. It originates from Delft University of Technology and is distributed under the GNU LGPL. The Parallel Data Assimilation Framework (PDAF, http://pdaf.awi.de) is also focused on Kalman filters, with special attention given to parallelism. It is written in Fortran 90 and distributed by the Computing Center of the Alfred Wegener Institute under the GNU GPL. The Data Assimilation Research Testbed (DART, http://www.image.ucar.edu/DAReS/DART), from the National Center for Atmospheric Research (NCAR), provides another Fortran implementation of the Kalman ensemble filter. Several major organizations involved in data assimilation have their internal solutions. For instance, the European Centre for Medium-Range Weather Forecasts (ECMWF) develops the Object-Oriented Prediction System (OOPS), written in C++ and Fortran. Different tools are also available to help implementing data assimilation algorithms. An example is the open-source coupler PALM (http://www.cerfacs.fr/ globc/PALM_WEB) distributed by Cerfacs. It can ease the coupling between a data assimilation algorithm, a numerical model and an observation module.

Realizing that data assimilation principles and practical experiences present a generic character was the starting point of Verdandi. In order to guarantee genericity and performance at the same time, we took different directions than other data assimilation software. The library also provides unique tools to carry out data assimilation experiments. The design of the library is detailed below.

We implemented a wide range of estimation algorithms within a common modularity paradigm, with the major objective to provide biophysical personalization tools for the various abovediscussed models, but more generally to be of interest for the modeling community concerned with merging models and data. Despite their diversity, all these models have in common the large dimension of the corresponding numerical systems to be solved, and the Verdandi library is targeted at providing generic data assimilation procedures well-adapted to large-dimensional systems. This library gathers until now optimal interpolation, Kalman filter and non-linear extensions, reduced-order versions, ensemble Kalman filter, non-linear filters in probabilistic or deterministic forms but also 4D-Var adjoint-based variational strategies.

The Verdandi library is distributed1 under the GNU LGPL license. It is based on generic C++ programming -a choice guided by performance -but is directly linkable with Python programs in order to ease scripting, debugging, visualization of results, and so on. Full compliance with the C++ standard warrants straightforward portability, and it has been thoroughly tested on Linux, MacOS, and Windows platforms.

Library overall design

Following the above methodological discussion, we identify three main components in the overall design of Verdandi :

-the numerical model, which provides the state vector X n , the dynamics with A n+1|n (and possibly its tangent linear and adjoint versions), and statistics for the initial state error P and the model error (when applicable); -the observation manager, which provides the observations Z n , their error covariance matrix W n and the observation operator H n ;

-the data assimilation algorithm, which drives the model and the observation manager along the assimilation procedure.

Each component is implemented in a separate C++ class. A class is an entity that contains variables, called attributes, and member functions that usually operate on the attributes. For instance, the model class has the state vector X n as attribute, and a member function called GetState to access this state vector.

An assimilation algorithm encapsulates the model and the observation manager. It is in charge of all assimilation computations, such as the correction step (25b) in the Kalman filter. It delegates the specific computations to the model and the observation operator, such as the time integration over one time step ( 4). An assimilation algorithm is implemented with essentially no assumption on the inner implementation of the model and the observation manager. Verdandi provides several data assimilation algorithms and expects that the user plugs a model and an observation manager with specified interfaces. Figure 1 illustrates the links between the three components. In Verdandi, the data assimilation algorithm is implemented in a class that encapsulates a model and an observation manager. The model defines the state vector, the dynamics, the associated error statistics, . . . The observation manager provides the observations, the observation operator, the associated error statistics, . . . These variables are accessed through member functions like GetState or GetObservation.

Data assimilation algorithm

Standard interfaces

The model and the observation manager are supposed to provide a standard interface so that Verdandi algorithms may be applied with them. In the model interface, one may find, among by the user. Out of the box, Verdandi supports a wide range of data structures, from dense to sparse structures, from contiguous to distributed memory blocks. For example, within the model, the state vector may be composed of a collection of non-contiguous memory blocks -one block per physical variable. In case even more structures were needed, advanced users could define their own data structures, provided they would implement the corresponding basic linear algebra.

The core library is implemented in C++, but the core model and observation manager can be implemented in another language. Users however need to implement the adequate C++ interface to their software. This task is rather straightforward when the software is in Fortran, C or C++. Verdandi also provides a C++ interface to a generic Python model, so that writing the interface to a Python model should be an easy task. Note that the high-dimensional data (especially the state vector) can always be provided to C++ by reference (or with pointers) so that no memory duplication is required, even with an underlying Python model.

In order to guarantee the best performance along with the genericity, the type of the data structures must be known at compile time. Hence an assimilation algorithm is always a C++ class template, with the types of the model and observation manager as arguments. The data types are defined by the user, inside the model and the observation manager, with typedef declarations. All key variables can have their own data structure -for example, the state vector and the observation vector may be stored in two different formats.

In addition to its C++ core (mainly, the assimilation algorithms available in C++ classes), Verdandi supports the automatic generation of Python interfaces. The interfaces are generated by SWIG, for any assimilation algorithm, and also for the numerical model and the observation manager. The C++ interfaces are first compiled, and SWIG generates wrappers in Python, so that all member functions of the model, the observation manager and the assimilation algorithm are exposed in Python.

Other features

All Verdandi objects are configured with Lua scripts. It means that the configuration files are interpreted and can execute non-trivial operations (system calls, numerical computations, . . . ).

The library includes perturbation schemes which are primarily used in Monte Carlo simulations and Ensemble Kalman filter. The perturbations can be applied to input parameters of the model. Input fields can be perturbed with spatially correlated perturbations.

Different tools are provided to help users implement their interfaces. One tool collects and aggregates observations within given time spans, so that observations can be assimilated even if the model times do not coincide with the observation times. Another tool allows the classes (model, observation manager, assimilation algorithm and possibly others) to send/receive messages to/from one or several other classes. In particular, it provides a direct communication channel between the model and the observation manager, which is otherwise impossible since both objects are entirely driven by the assimilation algorithm. Once a message is received by an object, this can trigger an action that is not part of the assimilation algorithm.

Conclusions

We have presented the fundamental principles of data assimilation underlying the Verdandi library, and how they are articulated with the modular architecture of the library. This translates -in particular -into the definition of standardized interfaces through which the data assimilation library interoperates with the model simulation software and the observation manager.

We also discussed various examples of data assimilation applied to the personalization of biophysical models, in particular for cardiac modeling applications within the euHeart European project. Whereas such applications are somewhat specific in some respects -for example in the type of data considered -this also shows that they can benefit from data assimilation methodologies valid within a wider scope. This justifies pursuing the development of the library in a long-term perspective, for the benefit of the VPH community -for example Verdandi is now also used and further developed in the VPH-Share European project -as well as in other application fields and in the industry, which is why the LGPL open-source license was selected for the software distribution.

Figure 1 :

 1 Figure1: In Verdandi, the data assimilation algorithm is implemented in a class that encapsulates a model and an observation manager. The model defines the state vector, the dynamics, the associated error statistics, . . . The observation manager provides the observations, the observation operator, the associated error statistics, . . . These variables are accessed through member functions like GetState or GetObservation.

http://verdandi.gforge.inria.fr

Acknowledgments: This work has been partially supported by the European Commission (FP7-ICT-2007-224495: euHeart and FP7-ICT-2009-269978: VPH-Share).

other member functions: GetState() that returns a reference on the state vector; Forward() that carries out the time integration for one time step; ApplyTangentLinearOperator(x) that applies the tangent linear operator to the vector x; GetStateErrorVarianceRow(i, v) that returns in v the row i of the state error covariance matrix. A simulation without assimilation may be carried out with the following C++ lines: Similarly, the observations and associated variables are provided by a few member functions of the observation manager. Notice that the observation operator depends on the model since it maps from the model state space into observation space. As a consequence, several member functions of the observation manager take the model as argument. Among the member functions of the interface, one may find: SetTime(model, t) that prepares the observation manager so that subsequent calls return values at time t; GetObservation(y) that returns the observation vector; ApplyOperator(x, y) that applies the observation operator; GetErrorVariance() that returns a reference to observation error covariance matrix.

The optimal interpolation, which replaces the state vector with BLUE whenever observations become available, is roughly implemented with these lines:

Implementation strategies

The algorithms are implemented with few assumptions on the data structures. The state vector, the error statistics, the observations, and so on, are provided in data structures defined


Mehwish Alam
email: mehwish.alam@inria.fr

Adrien Coulet
email: adrien.coulet@inria.fr

Amedeo Napoli
email: amedeo.napoli@inria.fr

Malika Smaïl-Tabbone
email: malika.smail@inria.fr

Formal Concept Analysis Applied to Transcriptomic Data

Keywords: Formal Concept Analysis, Knowledge Discovery, Data Integration, Transcriptomic Data

Identifying functions or pathways shared by genes responsible for cancer is still a challenging task. This paper describes the preparation work for applying Formal Concept Analysis (FCA) to biological data. After gene transcription experiments, we integrate various annotations of selected genes in a database along with relevant domain knowledge. The database subsequently allows to build formal contexts in a flexible way. We present here a preliminary experiment using these data on a core context with the addition of domain knowledge by context apposition. The resulting concept lattices are pruned and we discuss some interesting concepts. Our study shows how data integration and FCA can help the domain expert in the exploration of complex data.

Introduction

Over past few years, large volumes of transcriptomic data were produced but their analysis remains a challenging task because of the complexity of the biological background. In the field of transcriptomics, biologists analyze routinely the transcription or expression of genes in various situations (e.g., in tumor samples versus non-tumor samples).

Some earlier studies aimed at retrieving sets of genes sharing the same transcriptionl behaviour with the help of Formal Concept Analysis (see, e.g., [START_REF] Kaytoue-Uberall | Two FCA-Based Methods for Mining Gene Expression Data[END_REF][START_REF] Rioult | Using Transposition for Pattern Discovery from Microarray Data[END_REF][START_REF] Rioult | Mining Concepts from Large SAGE Gene Expression Matrices[END_REF]). Further studies analyze gene expression data by using gene annotations to determine whether a set of differentially expressed genes is enriched with biological attributes [START_REF] Gabriel | Characterizing gene sets with FuncAssociate[END_REF][START_REF] Doniger | MAPPFinder: using Gene Ontology and GenMAPP to Create a Global Gene-expression Profile from Microarray Data[END_REF][START_REF] Zhong | GoSurfer: a graphical interactive tool for comparative analysis of large gene sets in Gene Ontology space[END_REF]. Many useful resources are available online and several efforts have been made for integrating heterogeneous data [START_REF] Michael | The 2012 Nucleic Acids Research Database Issue and the online Molecular Biology Database Collection[END_REF][START_REF] Khatri | Ontological Analysis of Gene Expression Data: Current Tools, Limitations, and Open Problems[END_REF]. A recent example is of the Broad Institute where biological data were gathered from multiple resources to get thousands of predefined gene sets stored in the Molecular Signature DataBase, MSigDB [START_REF] Liberzon | Molecular Signatures Database (MSigDB) 3.0[END_REF]. A predefined gene set is a set of genes known to have a specific property such as their position on the genome, their involvement in a biological process (or a molecular pathway) etc. Subsequently, given an experimental gene list as input the GSEA (Gene Set Enrichment Analysis) program is used to asses whether each predefined gene set (in the MSigDB database) is significantly present in the input list by computing an enrichment score [START_REF] Subramanian | Gene Set Enrichment Analysis: A Knowledge-based Approach for Interpreting Genome-wide Expression Profiles[END_REF].

In this paper, we are interested in applying knowledge discovery techniques for analyzing a differentially expressed gene set and identifying functions or pathways shared by these genes assumed to be responsible for cancer. Knowledge discovery aims at extracting relevant and useful knowledge patterns from a large amount of data. It is an interactive and iterative process involving a human (analyst or domain expert) and data sources. We show how various gene annotations and domain knowledge are integrated in a database which is then queried for building in a flexible way formal contexts. We present here a preliminary experiments using these data. It was performed on a core context with the addition of domain knowledge (by context apposition). The considered domain knowledge are the hierarchical relationships between molecular pathways. Pruning the obtained lattices allows us to retrieve interesting concepts which we discuss. The results obtained from both experiments are also compared.

The plan of the paper is as follows: Section 2 introduces Formal Concept Analysis, Section 3 explains the data resources which are integrated, Section 4 focuses on the application of FCA, Section 5 discusses the results and Section 6 concludes the paper and presents future Work.

Formal Concept Analysis

We introduce here the basics of Formal Concept Analysis that are needed to understand what follows. Let G and M be the set of objects and set of attributes respectively and I be the relation between the objects and the attributes I ⊆ G × M , where g ∈ G, m ∈ M , gIm is true iff the object g has the attribute m. The triple K = (G, M, I) is called a formal context. If A ⊆ G, B ⊆ M are arbitrary subsets, then a Galois connection denoted by is given by:

A := {m ∈ M | gIm ∀ g ∈ A} (1)
B := {g ∈ G | g I m ∀ m ∈ B} (2)
FCA framework is fully described in [START_REF] Ganter | Formal Concept Analysis: Mathematical Foundations[END_REF]. FCA helps in defining concepts which are composed of a maximal set of objects sharing a maximal set of attributes. However, given an input context, the resulting concept lattice can be very large leading to computational and interpretation problems. In order to have reduced and meaningful concepts, one can select concepts whose support is greater than a certain threshold, i.e., the iceberg lattice. For a concept (A,B), the support is the cardinality of the extent A. An alternative is to use the notion of stability that was proposed in [START_REF] Kuznetsov | On stability of a Formal Concept[END_REF][START_REF] Roth | Towards Concise Representation for Taxonomies of Epistemic Communities[END_REF]. The stability index measures how much the concept intent depends on particular objects of the extent.

Complex Biological Data Integration

In this section, we introduce and describe the biological data on which we are working.

Molecular Signature Database (MSigDB)

Molecular Signature Database (MSigDB) is an up-to-date database which contains data from several resources such as KEGG, BIOCARTA, REACTOME, and Amigo [START_REF] Liberzon | Molecular Signatures Database (MSigDB) 3.0[END_REF]. It is a collection of 6769 predefined gene sets. A predefined gene set is a set of genes having a specific property such as their position on the genome (e.g., the genes at position chr5q12, i.e., band 12 on arm q of chromosome 5), their involvment in a biological process or a molecular pathway (e.g., the genes which are involved in the KEGG APOPTOSIS pathway)... A pathway is a series of actions among molecules in a cell that leads to a certain change in a cell. KEGG is a database storing hundreds of known pathways4 . Besides, the MSigDB gene sets are grouped into five categories (Table 1). For instance, all the gene sets which are defined on the basis of gene position belong to the category C1. The category C5 groups the gene sets defined on Gene Ontology (GO) terms annotating the genes (with respect to their molecular function or their housing cellular component).

For our study, we used MSigDB Version 3.0. One entry, shown below in XML format, describes the gene set corresponding to the GO term 'RNA Polymerase II Transcription Factor Activity Enhancer Binding' (all the attribute names are underlined). The Members attribute contains the list of gene symbols belonging to the gene set. MSigDB was chosen as the main source for describing genes because it gathers up-to-date informations about many aspects of human genes.

<GENESET Standard Name ="RNA Polymerase II Transcription Factor Activity Enhancer Binding" Systematic Name = "M900" Historical Names ="" Organism ="Homo sapiens" Geneset Listing URL ="" Chip = "Human Gene Symbol" Category Code ="c5" Sub Category Code ="MF" Contributor ="Gene Ontology" Contributor Org ="GO" Description Brief ="Genes annotated by the GO term GO:0003705. Functions to initiate or regulate RNA polymerase II transcription by binding an enhancer region of DNA." Description Full ="" Members =" MYOD1, TFAP4, EPAS1, RELA, MYF5, MYEF2, NFIX, PURA, HIF1A" Members Symbolized = "MYOD1, TFAP4, EPAS1, RELA, MYF5, MYEF2, NFIX, PURA, HIF1A" Members EZID =" 7023, 2034, 5970, 3091" Members Mapping = " MYOD1, 4654-TFAP4, TFAP4, 7023-EPAS1, EPAS1, 2034-RELA, RELA, 5970-MYF5, MYF5, 4617-MYEF2, MYEF2, 50804-NFIX, NFIX, 4784-PURA, PURA, 5813-HIF1A" Status ="public" > </GENESET>

Domain Knowledge

Besides the gene annotations included in MSigDB, many types of domain knowledge are interesting to use when analyzing genes. The first type of such do- AmiGO main knowledge are the hierarchical relationships between GO terms or between KEGG pathways. Indeed, the KEGG hierarchy for human groups the KEGG pathways into 40 categories and 6 upper level categories. Figure 1 illustrates the KEGG hierarchy detailing on one upper-level category and one category. In our study we have genes described by pathways involving them which may in turn be present in some category of pathways. For example, if a gene is involved in a pathway apoptosis it will also be in the category 'Cell Growth and Death'. In order to facilitate the knowledge discovery, it is important to identify the relevant data sources, organize, and integrate the data at one single database. In our case, the relevant primary data sources are MSigDB, KEGG PATHWAYS database, and AmiGO database.

From Data to Knowledge

Once the data are integrated in our database the next step is to build formal contexts for applying FCA. Our experiment focuses on applying FCA to a core context describing genes by MSigDB-based attributes and shows its extension based on the addition of domain knowledge.

Test Data Sets

The experiments described here are based on three published sets of genes corresponding to Cancer Modules defined in [START_REF] Segal | A Module Map Showing Conditional Activity of Expression Modules in Cancer[END_REF]. The authors compiled gene sets from various resources and a large collection of micro-array data related to cancers. These modules correspond to gene sets whose expression significantly change in a variety of cancer conditions (they are also defined as MSigDB gene sets in the C4 category). Our test data are composed of three lists of genes corresponding to the Cancer Modules 1 (Ovary Genes), 2 (Dorsal Root Ganglia Genes), and 5 (Lung Genes).

Using FCA for Analyzing Genes

We apply FCA for analyzing a context describing genes of each Cancer Module with MSigDB-based attributes. Table 2 shows five genes (involved in Cancer Module 1) as a set of objects described by attributes which are the memberships to gene sets from MSigDB. For example, CCT6A is in the set of genes (gene set) whose standard name is Reactome Serotonin Receptors. Interestingly, by querying our integrated database the analyst is able to select the predefined gene sets to include in the formal context.

In order to extend the analysis of a list of genes, we need to take into account the domain knowledge. Hence, the same experiment was conducted with the addition of the KEGG hierarchy knowledge to the core contexts resulting in three extended contexts. All KEGG categories and upper-level categories were added as a set of attributes. If a gene is member of a KEGG pathway which in turn belongs to a category and an upper level category then a cross '×' is added in the corresponding cells in the extended context.

Table 2 shows five genes (from Cancer Module 1) with the addition of one KEGG category (kc) and one KEGG upper level category (kuc). In the given example CCT6A is involved in pathway KEGG PPAR Signaling Pathway which belongs to the category kc:Endocrine System and upper level category kuc:Organismal Systems. The lattices were generated and the statistics for each Cancer Module are given in Table 3. The concepts were filtered and ranked based on same criteria as in the first experiment.

× × PSPHL × × × × CCT6A × × QNGPT1 × × × × MYC × ×

Results

In this study, biologists are interested in links between the input genes in terms of pathways in which they participate, relationship between genes and microRNAs etc. We obtained concepts with shared transcription factors, pathways, positions of genes and some GO terms. After the selection of concepts with higher support, we observed that there were some concepts with pathways from KEGG and RE-ACTOME as their intent. These pathways are either related to cell proliferation or apoptosis (cell death). The addition of domain knowledge effectively gives an opportunity to obtain the pathway categories shared by larger sets of of genes. 4) .

To sum up, we were able to discover interesting biological properties of subsets of genes in the three test data sets. As for example, the Focal Adhesion pathway was found to be associated to 17 genes in modules 1 and 2; the KEGG category Immune System was found to be shared by 11 to 25 genes in the three cancer modules (Table 4). Given the test data sets, these results are hopeful and constitute interesting positive control. This confirms that FCA-based analysis offers a powerful procedure to deeply explore sets of genes.

Conclusion and Future Work

Our study shows how Formal Concept Analysis can be applied to complex biological data. Data integration and FCA give the flexibility of using various types of attributes (pathways, GO terms, positions, microRNAs and Transcription Factor Targets) for analyzing a list of genes. Our approach gives an insight into how domain knowledge can be introduced in the analysis with the help of FCA. As for future work, we plan to apply our approach to experimental gene lists and take into account gene-gene relationships (physical Protein Protein Interactions), term-term relationships (Gene Ontology relationships, namely is-a, part-of, and regulates) and relationships between gene positions. Moreover, in order to efficiently deal with the relationships present within the data we can use Relational Concept Analysis.

Fig. 1 .

 1 Fig. 1. Hierarchical Relationship in KEGG

Table 2 .

 2 A Toy Example of Formal Context with Domain l u l a r C o m p o n e n t A s s e m b l y c h r 5 q 1 2 k c : E n d o c r i n e S y s t e m k u c : O r g a n i s m a l S y s t e m s BTB03 ×

Table 1 .

 1 Categories of MSigDB Gene Sets

	Category	Description	Data Provenance
	C1: Positional Gene	Location of the gene on the	Broad Institute
	Sets	chromosome.
	C2: Curated Gene	Pathways	KEGG,	REAC-
	Sets		TOME, BIOCARTA
	C3: Motif Gene Sets microRNAs, Transcription	Broad Institute
		Factor Targets.
	C4: Computational	Cancer Modules	Broad Institute
	Gene Sets
	C5: Gene Ontology	Biological Process, Cellu-
	(GO) Gene Sets	lar Components, Molecular
		Functions

Table 3 .

 3 Concept Lattice Statistics for the Cancer Modules with Domain Knowledge

	Data Sets No. of Genes No. of Attributes No. of Concepts Levels
	Module 1	361	3496	9,588	12
	Module 2	378	3496	6,508	11
	Module 5	419	3496	5,004	12

Table 4

 4 shows the top-ranked concepts found in each module. For example, in module 5, we have confirmation that Cytokine Cytokine Receptor Interaction pathway comes under the category Signaling Molecules and Interaction and upper level category Environmental Information Processing (Concept ID 4938). The absolute support and stability of the concept containing only the category Signaling Molecules and Interaction and upper level category Environmental Information Processing as its intent are higher (Concept ID 4995, Table

Table 4 .

 4 Top-ranked Concepts with Domain Knowledge

	Dataset Concept	Intents					Absolute	Stability
	ID						Support
	Module 1 9585	M2192:GGGAGGRR V$MAZ Q6		51	0.99
	9571	M2598:GO Membrane Part			27	0.99
	9566	kc:Immune System, kuc:Organismal Sys-	25	0.99
		tems
	9402	chr19q13					10	0.99
	9078	M10792:KEGG	MAPK	Signaling	12	0.87
		Pathway,	kc:Signal	Transduction,
		kuc:Environmental	Information	Pro-
		cessing
	Module 2 6502	M2192:GGGAGGRR V$MAZ Q6		44	0.99
	6496	kc:Immune System, kuc:Organismal Sys-	15	0.99
		tems
	6388	chr6p21					10	0.97
	6335	M10792:KEGG	MAPK	Signaling	11	0.89
		Pathway,	kc:Signal	Transduction,
		kuc:Environmental	Information	Pro-
		cessing
	Module 5 5002	kuc:Cellular Processes			48	0.99
	4995	kc:Signaling Molecules and Interaction,	26	0.99
		kuc:Environmental Information Process-
		ing
	4933	chr19q13					11	0.99
	4985	kc:Immune System, kuc:Organismal Sys-	11	0.99
		tems
	4938	M9809:KEGG Cytokine Cytokine Recep-	11	0.87
		tor Interaction, kc:Signaling Molecules and
		Interaction, kuc:Environmental Informa-
		tion Processing

http://www.genome.jp/kegg/pathway.html


Mehwish Alam
email: mehwish.alam@inria.fr

Adrien Coulet
email: adrien.coulet@inria.fr

Amedeo Napoli
email: amedeo.napoli@inria.fr

Malika Smaïl-Tabbone
email: malika.smail@inria.fr

Formal Concept Analysis Applied to Transcriptomic Data

Keywords: Formal Concept Analysis, Knowledge Discovery, Transcriptomic Data

Identifying functions shared by genes responsible for cancer is a challenging task. This paper describes the preparation work for applying Formal Concept Analysis (FCA) to complex biological data. We present here a preliminary experiment using these data on a core context with the addition of domain knowledge. The resulting concept lattices are explored and some interesting concepts are discussed. Our study shows how FCA can help the domain experts in the exploration of complex data.

Introduction

Over past few years, large volumes of transcriptomic data were produced but their analysis remains a challenging task because of the complexity of the biological background. Some earlier studies aimed at retrieving sets of genes sharing the same transcriptional behavior with the help of Formal Concept Analysis [START_REF] Kaytoue-Uberall | Two FCA-Based Methods for Mining Gene Expression Data[END_REF][START_REF] Rioult | Using Transposition for Pattern Discovery from Microarray Data[END_REF]. Further studies analyze gene expression data by using gene annotations to determine whether a set of differentially expressed genes is enriched with biological attributes [START_REF] Berriz | Characterizing gene sets with FuncAssociate[END_REF][START_REF] Doniger | MAPPFinder: using Gene Ontology and GenMAPP to Create a Global Geneexpression Profile from Microarray Data[END_REF]. Several efforts have been made for integrating heterogeneous data [START_REF] Galperin | The 2012 Nucleic Acids Research Database Issue and the online Molecular Biology Database Collection[END_REF]. For example, at the Broad Institute, biological data were recently gathered from multiple resources to get thousands of predefined genesets stored in the Molecular Signature DataBase (MSigDB) [START_REF] Liberzon | Molecular Signatures Database (MSigDB) 3.0[END_REF]. A predefined geneset is a set of genes known to have a specific property such as their position on the genome, their involvement in a molecular pathway etc.

This paper focuses on the preparation of biological data to data mining guided by domain knowledge. The objective is to apply knowledge discovery techniques for analyzing a list of differentially expressed genes and identifying functions or pathways shared by these genes assumed to be responsible for cancer. Section 2 explains the proposed approach for FCA-based analysis of biological data. Section 3 focuses on the conducted experiment. Section 4 discusses the results. Section 5 concludes the paper.

The Proposed Framework

We rely on the standard definition of FCA fully described in [START_REF] Ganter | Formal Concept Analysis: Mathematical Foundations[END_REF] and adapt it according to the current problem. Let G be the set of genes {g 1 , g 2 , g 3 , ..., g n }, and M be a set of attributes of MSigDB for describing genes. M will be considered as a partition of three points of view,

M = M 1 ∪ M 2 ∪ M 3 , with M i ∩ M j = ∅ whenever i = j.
The first set of attributes M 1 refers to four types of attributes, "Location", "Pathway", "Transcription Factors" and "GO Terms" (see Table 1). For our convenience we have named MSigDB categories as types of attributes and used only C 1 , C 2 , C 3 and C 5 . The category C 4 was not used as it keeps information on sets of genes related to a certain kind of cancer, which is not useful for the current problem. Thus we have a first context K 1 = (G, M 1 , I 1 ) where I 1 denotes the relation stating that gene g i has an attribute m j in M 1 . The second set of attributes M 2 is related to the so-called "categories" where a category makes reference to a set of attributes with the "Pathway" type. For example, "Cell Growth and Death" is an example of category (see Figure 1). The categories in M 2 determine a second context, K 2 = (G, M 2 , I 2 ) where M 2 is the set of categories and I 2 denotes the relation between a gene and a category. It can be noticed that the categories are only related to the "Pathway" type and that they can be considered as domain knowledge.

Types of

Moreover, the third set of attributes, namely M 3 , refers to the so-called "upper categories", which are defined as groupings of categories. Actually, we have for the type "Pathway" a hierarchy of categories with two levels, categories and upper categories (see Figure 1). The upper categories in M 3 define a third context K 3 = (G, M 3 , I 3 ) where M 3 is the set of upper categories and I 3 denotes the relation between gene g i and an upper level category m j . Upper categories are also related to the "Pathway" type and as categories, they can be considered as domain knowledge too. Now we consider the apposition of the three contexts K 1 , K 2 and K 3 , which yields the final context K = (G, K 1 ∪ K 2 ∪ K 3 , I 1 ∪ I 2 ∪ I 3 ). For example the context in Table 2 shows five genes described by attributes of M 1 , M 2 and M 3 .

Using FCA for Analyzing Genes

The framework described above was applied on three published sets of genes corresponding to Cancer Modules defined in [START_REF] Segal | A Module Map Showing Conditional Activity of Expression Modules in Cancer[END_REF]. Our test data are composed of three lists of genes corresponding to the so-called "Cancer Module 1" (Ovary Genes), "Cancer Module 2" (Dorsal Root Ganglia Genes), and "Cancer Module 5" (Lung Genes). For example, "PSPHL" is one gene with "Pathway" attribute as "PPAR Signaling" which belongs to category "kc:Endocrine System" and upper category "kuc:Organismal System". Considering the three lists of genes given by "Cancer Module 1", "Cancer Module 2" and "Cancer Module 5", we built three different contexts having the same form as the context in Table 2). Then we obtained three associated concept lattices with the help of the Coron Plate-form (http://coron.loria.fr). The concept lattice for Table 2 is given in Figure 2. The global characteristics of the three concept lattices are given in Table 3.

The exploration of a given concept lattice is carried out following the "Iceberg metaphor", i.e., the lattice is explored level by level according to the support of

× × PSPHL × × × × CCT6A × × QNGPT1 × × × × MYC × × Table 2.
A toy example of formal context including domain knowledge.

each concept, where the the support of a concept is the cardinality of the extent.

In addition, we also used stability for extracting interesting frequent and stable concepts [START_REF] Kuznetsov | On stability of a Formal Concept[END_REF].

Results

In this study, biologists are interested in links between the input genes in terms of pathways in which they participate, relationships between genes and their positions etc. We obtained concepts with shared transcription factors, pathways, locations of genes and GO terms. After the selection of concepts with a high support (≥ 10), we observed that there were some concepts with pathways either related to cell proliferation or apoptosis (expert interpretation). The addition of domain knowledge gives an opportunity to obtain the pathway categories shared by larger sets of genes (as categories and upper categories are there for maximizing the grouping of objects, see below). Table 4 shows the top-ranked concepts found in each module. For example, in Table 4, we have the concept C 4938 :(KEGG Cytokine Cytokine Receptor Interaction, kc:Signaling Molecules and Interaction, kuc:Environmental Information Processing) and the concept C 4995 :(kc:Signaling Molecules and Interaction, kuc:Environmental Information Processing). These two concepts are such as C 4938 ≤ C 4995 , meaning that C 4995 has greater support than C 4938 . Moreover, we observed that the introduction of categories and upper categories in the global context allows us to consider concepts that otherwise would not be frequent. Actually, the role of categories and upper level categories is to facilitate the observation of sets of related genes. This is a general way of obtaining larger sets of objects to interpret. When available, one can introduce a hierarchy of attributes -this is domain knowledgeand then insert the levels of each attribute in this hierarchy as a new attribute in the context. As a result, some classes of objects, that could not emerge before, will appear based on these hierarchical indications. Given the test data sets, the preliminary results obtained here constitute an interesting and positive control, and confirm that FCA-based analysis offers an efficient and practical procedure to explore complex and large sets of genes.

Conclusion

The preliminary study presented here shows how FCA can be applied to complex biological data and can give flexibility in using various types of attributes for analyzing a list of genes. In addition, domain knowledge can be introduced and guide the analysis. As for future work, we plan to take into account relationships between genes and between terms (Gene Ontology relationships) and use the framework of relational concept analysis.

Fig. 1 .

 1 Fig. 1. Categories and upper categories in KEGG.

  c a t i o n ) k c : E n d o c r i n e S y s t e m k u c : O r g a n i s m a l S y s t e m s BTB03 ×

Fig.

  Fig.2. The concept lattice corresponding to table 2.

Table 1 .

 1 Types of attributes from MSigDB

	Attributes	Description	Data Provenance
	C1: Positional Gene Sets	Location of the gene on the chro-	Broad Institute
		mosome.
	C2: Curated Gene Sets	Pathway	KEGG,	REAC-
			TOME, BIOCARTA
	C3: Motif Gene Sets	Transcription Factors	Broad Institute
	C4: Computational Gene	Cancer Modules	Broad Institute
	Sets
	C5: Gene Ontology (GO)	Biological Process, Cellular Com-	AmiGO
	Gene Sets	ponents, Molecular Functions

Table 3 .

 3 Concept lattice statistics for the cancer modules with domain knowledge.

	Data Sets No. of Genes No. of Attributes No. of Concepts Levels
	Module 1	361	3496	9,588	12
	Module 2	378	3496	6,508	11
	Module 5	419	3496	5,004	12

2. The concept lattice corresponding to table 2.

Table 4 .

 4 Top ranked concepts for each cancer module

	Dataset Concept	Intents	Absolute	Stability
	ID		Support
	Module 1 9585	GGGAGGRR V$MAZ Q6	51	0.99
	9571	GO Membrane Part	27	0.99
	9566	kc:Immune System, kuc:Organismal Systems	25	0.99
	9402	chr19q13	10	0.99
	9078	KEGG MAPK Signaling Pathway, kc:Signal Transduction,	12	0.87
		kuc:Environmental Information Processing
	Module 2 6502	GGGAGGRR V$MAZ Q6	44	0.99
	6501	AACTTT UNKNOWN	38	0.99
	6496	kc:Immune System, kuc:Organismal Systems	15	0.99
	6388	chr6p21	10	0.97
	6335	KEGG MAPK Signaling Pathway, kc:Signal Transduction,	11	0.89
		kuc:Environmental Information Processing
	Module 5 5002	kuc:Cellular Processes	48	0.99
	5000	GGGAGGRR V$MAZ Q6	44	0.99
	4995	kc:Signaling Molecules and Interaction, kuc:Environmental In-	26	0.99
		formation Processing
	4933	chr19q13	11	0.99
	4985	kc:Immune System, kuc:Organismal Systems	11	0.99
	4938	KEGG Cytokine Cytokine Receptor Interaction, kc:Signaling	11	0.87
		Molecules and Interaction, kuc:Environmental Information Pro-
		cessing