Haifa Ben Saber
email: bensaberhaifa@yahoo.com

Mourad Elloumi
email: mourad.elloumi@fsegt.rnu.tn

Mohamed Nadif
email: mohamed.nadif@parisdescartes.fr

Block Mixture Model for the Biclustering of Microarray Data

Keywords: biclustering, block expectation maximization algorithm, block classification expectation maximization algorithm, block mixture model, fuzzy strategy, microarray data

An attractive way to make biclustering of genes and conditions is to adopt a Block Mixture Model (BMM). Approaches based on a BMM operate thanks to a Block Expectation Maximization (BEM) algorithm and/or a Block Classification Expectation Maximization (BCEM) one. The drawback of these approaches is their difficulty to choose a good strategy of initialization of the BEM and BCEM algorithms. This paper introduces existing biclustering approaches adopting a BMM and suggests a new fuzzy biclustering one. Our approach enables to choose a good strategy of initialization of the BEM and BCEM algorithms.

I. INTRODUCTION

Analysis of gene expression data obtained from microarray experiments can be made through biclustering. Indeed, gene expression data are usually represented by a data matrix M (see Table.1), where the i th row represents the i th gene, the j th column represents the j th condition and the cell x ij represents the expression level of the i th gene under the j th condition. In general, subsets of genes are coexpressed only under certain conditions but behave almost independently under others. Discovering such coexpressions can be helpful to discover genetic knowledge such as genes annotation or genes interaction. Hence, it is very interesting to make a simultaneous clustering of rows (genes) and of columns (conditions) of the data matrix to identify groups of rows (genes) coherent with groups of columns (conditions), i.e., to identify clusters of genes that are coexpressed under clusters of conditions, or clusters of conditions that make clusters of genes coexpress. This type of clustering is called biclustering [START_REF] Cheng | Biclustering of expression data[END_REF]. A cluster made thanks to a biclustering is called bicluster. Hence, a bicluster of genes (resp. conditions) is defined with respect to only a subset of conditions (resp. genes). Thus, a bicluster is a subset of genes showing similar behavior under a subset of conditions. Let us note that a gene/condition can belong to more than one bicluster or to no bicluster.

In other words, a bicluster can be defined as follows: Let I={1, 2, …, n} be a set of indices of n genes, J={1, 2, …, m} be a set of indices of m conditions and X(I,J) be a data matrix associated with I and J. A bicluster associated with the data matrix X(I,J) is a couple (I',J') such that I'⊆ I and J'⊆ J.

Actually, biclustering is a special case of clustering. Indeed, in biclustering, genes are clustered according to their expression levels under a number of conditions, not necessarily all the conditions. While in clustering, genes are clustered according to their expression levels under all the conditions. Similarly, in biclustering, conditions are clustered according to the expression levels of a number of genes, not necessarily all the genes.

The biclustering problem can be formulated as follows:

Given a data matrix X, construct a bicluster B opt associated with X such that:

f(B opt ) = max B∈BC(X) f(B) ( 1
)
where f is an objective function measuring the quality, i.e., coherence degree, of a group of biclusters and BC(X) is the set of all the possible groups of biclusters associated with X.

Condition 1 … Condition j … Condition m Gene 1 x 11 … x 1j … x 1m … … … … … … Gene i x i1 … x ij … x im … … … … … … Gene n x n1 … x nj … x nm
Clearly, biclustering is a highly combinatorial problem with a search space size O(2 |I|+|J| ). In its general case, biclustering is NPhard [START_REF] Cheng | Biclustering of expression data[END_REF].

In this paper, we introduce existing biclustering approaches adopting a BMM and suggest a new fuzzy biclustering one. Our approach enables to choose a good strategy of initialization of BEM and BCEM algorithms.

The rest of this paper is organized as follows: In sections II and III, we briefly review the different types of biclusters and groups of biclusters. In section IV, we present the Mixture Model (BMM). Section V is devoted to describe our new fuzzy biclustering approach. Finally, in section VI, we present our conclusion and perspectives II. TYPES OF BICLUSTERS A bicluster can be in one of the following cases:

1. Bicluster with constant values: It is a bicluster where all the values are equal to a constant c:

x ij = c (2)
2. Bicluster with constant values on rows or columns:

• Bicluster with constant values on rows: It is a bicluster where all the values can be obtained by using one of the following equations:

x ij = c + a i ( 3
)
x ij = c ×ai ( 4
)
Where c is a constant and a i is the adjustment for the row i, 1 ≤ i ≤ n.

• Bicluster with constant values on columns: It is a bicluster where all the values can be obtained by using one of these equations:

x ij =c + b j (5)

x ij = c × b j ( 6
)
Where c is a constant and b w is the adjustment for the column j, 1≤ j ≤ m.

Bicluster with coherent values:

It is a bicluster that can be obtained by using one of the following equations:

x ij = c + a i + b j (7) x ij = c × a i ×b j (8)
4. Bicluster with linear coherent values: It is a bicluster where all the values can be obtained by using the following equation:

x ij = c × a i + b j (9)
5. Or, bicluster with coherent evolution: It is a bicluster where all the rows (resp. columns) induce a linear order across a subset of columns (resp. rows).

III. GROUPS OF BICLUSTERS

A group of biclusters can be in one of the following cases [START_REF] Jollois | Speed-up for the expectation-maximization algorithm for clustering categorical data[END_REF]:

1. Single bicluster ((Figure 1 A natural way to visualize a group of biclusters consists in assigning a different color to each bicluster and in reordering the rows and the columns of the data matrix so that we obtain a data matrix with colored blocks, where each block represents a bicluster.

IV. BLOCK MIXTURE MODEL

An attractive way to make biclustering of genes and conditions is to adopt a Block Mixture Model (BMM) [START_REF] Dharan | Biclustering of gene expression data using reactive[END_REF][START_REF] Madeira | Biclustering algorithms for biological data analysis: A survey[END_REF]. This presents the advantage of selecting an adequate model for each bicluster, which leads to biologically pertinent biclustering of gene expressions data.

A. Definition of the model

Mixture Models (MM) [START_REF] Dharan | Biclustering of gene expression data using reactive[END_REF][START_REF] Madeira | Biclustering algorithms for biological data analysis: A survey[END_REF] are widely adopted in a variety of areas and particularly in unsupervised learning. With a clustering approach based on a MM, it is assumed that the data to be clustered are generated by a mixture of underlying probability distributions in which each component represents a different cluster. Given observations X = (x 1 ,…, x n ), the probability density function (pdf) is defined as follows:

) ; ( ) ; ( K 1 k k i k k i x p x f α ϕ θ ∑ = = , ( 10
)
where ϕ k (x i ;α k ) is the density of an observation x i from the k th component (cluster) and the p k 's, 1 ≤ k ≤ K, are the corresponding cluster parameters. These densities belong to the same parametric family (spherical, diagonal or general family). A parameter p k , 1≤k≤K, is the probability that a gene belongs to the k th component (cluster), where K, assumed to be known, is the number of the components (clusters) in the mixture. The parameter of this model is the couple θ = (p,α), where p = (p 1 ,…,p K ) represents the mixing proportions and α = (α 1 , …, α K ) represents the parameters of each component (cluster). The mixture density of the observation x is expressed as follows: )

; ( ) ; ( K 1 k k i k k i x p x f α ϕ θ ∑ = Π = . ( 11
)
Let z=( z ik ) be the (n × K) binary classification matrix indicating the labels of rows into K clusters (z ik =1 if the row i belongs to the kth cluster and 0 otherwise). Govaert and Nadif [START_REF] Govaert | Clustering with block mixture models[END_REF] have shown that equation ( 10) can be written in the following way:

( ; ) ∑ ( ; ) ( ; , ) (12)
where Z is the set of all possible assignments z of I and, p( ; ) ∏ , and ( ; , ) ∏ ( ; ) , .

In the context of the biclustering problem, equation ( 13) takes the following expression:

( ; ) ∑ ( ; ) ( ; ) ( ; , , ) ( , ) Z W (14)
where w=(w jl ) corresponds to the (m × L) binary classification matrix indicating the labels of columns into L clusters (w jl =1 if the column j belongs to the lth cluster and 0 otherwise) while W denotes the set of all possible assignments w of J. A MM associated with equation ( 14) is called Block Mixture Model (BMM) [START_REF] Dharan | Biclustering of gene expression data using reactive[END_REF][START_REF] Madeira | Biclustering algorithms for biological data analysis: A survey[END_REF].

B. Soft and hard approaches

Approaches based on a BMM operate thanks to a Block Expectation Maximization (BEM) algorithm and/or a Block Classification Expectation Maximization (BCEM) one [START_REF] Scott | Clustering methods based on likelihood ratio criteria[END_REF].

Let's note that BCEM algorithm differs from BEM one by the fact that we have inserted a classification step between estimation and maximization ones.

1) Block EM algorithm

A Block EM algorithm (BEM) has the same convergence properties as Generalized EM one [START_REF] Govaert | An EM algorithm for the block mixture model[END_REF]. And, like an EM one, a BEM algorithm converges slowly in some situations. The second important drawback of this kind of algorithms is that their solutions can highly depend on their starting position and, consequently, produce sub-optimal maximum likelihood estimates. To cope with this high dependency, Govaert and Nadif [START_REF] Govaert | An EM algorithm for the block mixture model[END_REF] have proposed to use the em-EM strategy to initialize parameters estimations of a BEM algorithm.

This strategy consists in several short runs of BEM from random positions followed by a long run of BEM from the solution maximizing the likelihood. A BEM algorithm can be defined as follows:

Step 1: Set r: = 0, choose s from [2..n] and t from [2..m];

Choose z (0) , w (0) and θ (0) as initial values of z, w and θ

Step 2: Compute z (r+1) , w (r+1) and θ (r+1) thanks to the alternated application of EM on intermediate data matrices; for details see for instance [START_REF] Govaert | An EM algorithm for the block mixture model[END_REF].

Step 3: Iterate Step 2 until convergence.

2) Block CEM algorithm

As we said earlier, a BCEM algorithm differs from a BEM one by the fact that we have inserted a classification step between estimation and maximization ones. That is why, we need to adopt classification log-likelihood function L c (z,w,θ) instead of loglikelihood one L(z,θ). A BCEM algorithm [START_REF] Govaert | An EM algorithm for the block mixture model[END_REF] differs from a BEM one by the fact that in Step 2 we call a CEM algorithm instead of an EM one.

C. Fuzzy approaches

Hard approaches adopting a BMM present the following drawbacks:

(i) All initializations strategies, including em-EM, of parameters estimation adopted by hard approaches are still inefficient and therefore do not contribute to generate pertinent groups of biclusters, i.e., groups of biclusters with optimal values of ML/CML function.

(ii) These approaches converge slowly and therefore are time consuming.

(iii) Concerning hard approaches that operate separately on the sets of genes I and the sets of conditions J (denoted by 2EM and 2CEM), the correlations that exist between subsets of genes and subsets of conditions are ignored.

That is why researchers opted for fuzzy approaches adopting a BMM to deal with the biclustering problem. In this context, Govaert and Nadif [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] proposed Block Fuzzy C-Methods (BFCM) based on Neal and Hinton interpretation of EM [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. The obtained approach enables to maximize a fuzzy clustering criterion. Miin-Shen Yang and Chih-Ying Lin [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF], on their side, proposed a Block Fuzzy K-Modes (BFKM) clustering approach based on Huang and Ng Fuzzy K-Modes (FKM) approach [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF].

1) Neal and Hinton interpretation of EM

As pointed out by Hathaway [START_REF] Hathaway | Another interpretation of the EM algorithm for mixture distributions[END_REF], in a classical MM context, an EM algorithm can be viewed as an alternated maximization of the following fuzzy clustering criterion:

( , θ) ( , θ) ∑ log , ( 15
)
where s is a n × g data matrix, with n is the number of genes and g is the maximum number of gene clusters to generate. The i th row of s represents the i th gene, the k th column represents the k th cluster and the cell s ik expresses the probability that the i th gene belongs to the k th cluster. The function

( , ) log( ( , ))
, is the fuzzy complete data log-likelihood. Indeed, maximizing L c with respect to s, and consequently F C with respect to s, yields the expectation step and maximizing L c with respect to θ, and consequently F C with respect to θ, yields the maximization one.

To extend this interpretation to all models in which an EM algorithm may be employed, Neal and Hinton [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] proposed to use the following function:

( , θ) ( ( , θ)) + ( ) (16)
where z represents the missing data, P is a distribution over the space of z and H is an entropy function. One can also relate F C to the Kullback-Leibler divergence between P(z) and ( ) ( ; , ), as follows:

( , θ) (θ) ( , P ) (17)
The alternated maximization of the F C function is then simple to set up:

(i) Maximization of F C (P,θ) with respect to P for fixed θ : Equation (17) leads to the minimization of ( , )) and then to P=P θ .

(ii) Maximization of F C (P,θ) with respect to θ for fixed P:

Equation ( 16) shows that θ maximizes the expectation

E p (L c (z,θ)).
These two steps are precisely those of an EM algorithm. Moreover, after the first step, we have F C (P,θ) = F C (P θ ,θ) = L M (θ), which shows that each iteration increases the loglikelihood L M (θ).

One can easily see that the Neal and Hinton interpretation extends Hathaway's relation (14). Indeed, for a classical MM, the missing data z represents the assignments of the genes of the set Z, and the assignments of each x i obtained in the expectation step are independent. Then, the conditional distribution P(z) is given by the vector that is the concatenation of the rows of the n×g data matrix s. Moreover, thanks to the Neal and Hinton interpretation [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF][START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF] of equations ( 15) and ( 16), E p (L c (z,θ)) can be transformed to L c (s,θ), and H(P) can be transformed to ∑ log , . For latent block model, Govaert and Nadif [START_REF] Jollois | Speed-up for the expectation-maximization algorithm for clustering categorical data[END_REF] adopted Neal and Hinton fuzzy criterion (15) defined by:

( , θ) ( ( , θ)) + ( ) (18) where ( , ; ) , + , + , , , ( ; )

2) Block Fuzzy C-Methods

In [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] the authors proposed Block Fuzzy C-Method (BFCM) as a new biclustering method based on a BMM and fuzzy cpartitions [START_REF] Ruspini | A new approach to clustering[END_REF]. The new optimized objective function is defined as follows: where c ik and d jl represent the fuzzy partitions of genes and conditions, respectively.

( ) ( ) ( ) , , , , , c c F c d L c d H c d θ θ = + (
To maximize F c (c,d;θ), in [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] the authors shown that this aim can be performed by maximizing alternatively two conditional fuzzy criteria F c (c,θ|d) and F c (d,θ|c). Their maximization can be viewed as an EM algorithm with a ML function associated with a classical MM on intermediate matrices. Thus, we can summarize BFCM algorithm as follows:

Step 1: Choose c (0) , d (0) and θ (0) as initial values of c, d and θ

Set r: = 0, choose s from [2..n] and t from [2..m];

Step 2: Compute c (r+1) , d (r+1) and θ (r+1) thanks to the alternated application of EM on intermediate data matrices.

Step 3: Iterate Step 2 until convergence.

We present now the Block Fuzzy K-Modes (BFKM) clustering algorithm [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF].

3) Block Fuzzy K-Modes algorithm

For clustering categorical data, i.e., types of data which may be divided into groups, Huang and Ng [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF] proposed the Fuzzy K-Modes (FKM) algorithm using simple matching dissimilarity measures. Huang and Ng [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF] added the concept of biclustering with FKM, and then proposed the Block FKM (BFKM) method.

Let Y = {y 1 ,..,y n } be a set of data. Let each data be defined by a set of attributes A 1 ,…,A m such that y i = {y i1 ,..,y in }, i = 1,…,n. Each attribute A j , j = 1, … , J describes a domain of values denoted by DOM(A j )={a j (1) , …, a j (Lj) } where L j denotes the cardinality of DOM(A j ). Suppose that ( , . . , ) is the centroid of the kth cluster where each component , … , for k=1,…,K , j=1,…,L with 1 and 0 for l'≠l, 1≤j≤m, 1≤l', l≤ L j . Huang and Ng [START_REF] Dharan | Biclustering of gene expression data using reactive[END_REF] used the following matching dissimilarity measure:

) ) .

We note is the cluster center, it is defined by: 1, if ∑ ∑ 0, otherwise

Where f' and l' are the numbers of the attribute j A (21) The BFKM clustering algorithm aims at minimize the following objective function:

  (a)), 2. Exclusive row and column group of biclusters (Figure 1 (b)), 3. Non-overlapping group of biclusters with checkerboard structure (Figure 1 (c)), 4. Exclusive rows group of biclusters (Figure 1 (d)), 5. Exclusive columns group of biclusters (Figure 1 (e)), 6. Non-overlapping group of biclusters with tree structure (Figure 1 (f)), 7. Non-overlapping non-exclusive group of biclusters (Figure 1 (g)), 8. Overlapping group of biclusters with hierarchical structure (Figure 1 (h)), 9. Or, arbitrarily positioned overlapping group of biclusters (Figure 1 (i)).

Figure 1

 1 Figure 1 Possible structures of a group of biclusters in a data matrix.

δ 1 For i = 1 ,

 11 The fuzzy k-modes clustering algorithms is to minimize the following objective function: …,n. The update equations for fuzzy k-modes are as follows: b>1,

TABLE 1 Gene

 1

Expression Data Matrix.

 ( ,,,,,Subject to ∑ 1 , for k=1,…, K and ∑ 1, for j=1,…, L. and are respectively the clusters centers, X is a data matrix with n genes and m conditions, and Y is the transpose of X. The update equations for BFKM are as follows:

)

) and (∑ (

) )

(23) We note and are the clusters centers, there are defined by:

Where f' and h' are the numbers of the attribute j A Thus, we can summarize BFKM algorithm as follows:

Step 1: Set r: = 0, choose k from [2..K], l from [2..L] and ε from [0..0.5];

Choose µ (0) and σ (0) as initial values of µ and σ

Step 2: Compute µ (r+1) , σ (r+1) , (r+1) , (r+1) from µ (r) , σ (r) , (r) and (r) thanks to equations (23), ( 24) and (25).

Step 3: If |µ (r+1) -µ (r) | + |σ (r+1) -σ (r) | < ε then stop else set r := r + 1, go to Step 2.

V. A NEW FUZZY ALGORITHM ADOPTING BMM

Both BFCM [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] and BFKM [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF] algorithms present the following drawbacks:

(i) They do not really apply the principal of fuzzy logic. Indeed, they do not adopt Fuzzy Membership Functions (FMF).

(ii) To act against this high dependency on its initial position, in [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] the authors propose to use the "em-EM" strategy.

Unlike BFCM and BFKM algorithms, our biclustering algorithm:

(i) Really adopts FMF (ii) Adopts fuzzy logic techniques in the strategy of initialization of parameters estimation. Indeed, fuzzy logic enables flexible initialization of parameters estimation.

Besides, we propose Fuzzy Logic Block Lazy EM algorithm (FLBLEM) to be adopted by our biclustering approach instead of BEM one. FLBLEM is based on a Lazy EM (LEM) [START_REF] Jollois | Speed-up for the expectation-maximization algorithm for clustering categorical data[END_REF] which is faster than EM. We can summarize our biclustering approach as follows:

During the first step, we fix a maximum number of biclusters, via the FMF, to initialize the parameters of our BMM.

During the second step, we cluster genes and/or conditions simultaneously, or consecutively, by applying a Block Lazy EM (BLEM) algorithm, if we are in an estimation phase.

During the last step, we tune the obtained group of biclusters either by merging biclusters or by creating others. If the algorithm does not converge, we repeat the second step.

V. CONCLUSION AND PERSPECTIVES

In this paper, we have reviewed existing biclustering approaches adopting a Block Mixture Model (BMM) and suggested a new fuzzy biclustering one.

Our approach enables to choose a good strategy of initialization of Block Expectation Maximization (BEM) and Block Classification Expectation Maximization (BCEM) algorithms.

This an ongoing work, as perspectives, we have to: (i) Handle missing values in microarray data, (ii) Speed-up our approach by defining a good strategy of choice of the initial number of biclusters.

(iii) Make an experimental study of our approach and compare the obtained results with those obtained by Block Fuzzy C-Methods (BFCM) [START_REF] Govaert | Fuzzy clustering to estimate the parameters of block mixture models[END_REF] and Block Fuzzy K-Modes (BFKM) [START_REF] Huang | A fuzzy k-modes algorithm for clustering categorical data[END_REF] ones.


Keywords: microarrays, classification, association

the classification and extraction of knowledge resulting from microarrays are a growing field which gave rise to several research tasks, however these works have not been able yet to be freed with the problems dependant between relations of all genes implied in breast cancer's disease. Moreover, the suggested methods generated some limiting dependency on the format of the treated data and extraction knowledge method. In this article, we present a method of classification and extraction of knowledge to handle data resulting from the microarrays. This method makes it initially possible to classify genes or the patients who have similar profiles of expression, and then to interpret the groups obtained. It is about a coupling between an algorithm of classification and an algorithm of extraction's rules associations. This work was validated by an expert. The obtained results seem to be interesting and promising owing to the fact that extracted knowledge can be used in the clinical medium to discriminate a disease or the state of a patient.

I. INTRODUCTION

he bioinformatics is a domain in full rise ; it results from he need to apply algorithms of biological data to analyze hem [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]. Thus, the principal bioinformatics task is to manage, exploit and analyze the production of voluminous genomics data [START_REF] Statnikov | A comprehensive evaluation of multicategory classi_cation methods for microarray gene expression cancer diagnosis[END_REF]. In our work, we are interested in the data analysis of genes or patients. The data used to carry out our analysis result from the microarrays: a technology simultaneously to analyze the activity level and the interactions of thousands of genes [START_REF] Trang | Data mining of gene expression microarray via weighted pre_x trees[END_REF]. The exploitation of the voluminous data by means of the computer programs is essential to extract interesting knowledge.

II. PROBLEMS

The traditional data-mining is used to extract data from the significant clustering of the attributes (columns) or objects (lines) or bi-clustering. This task gave place to many works [START_REF] Dr | Data analysis tools for DNA mi-croarrays[END_REF]3,[START_REF] Slonim | From patterns to pathways : gene expression data analysis comes of age[END_REF][START_REF] Alizadeh | Distinct types of 6di_use large B-cell lymphoma identi_ed by gene expression pro_ling[END_REF][START_REF] Yeung | Validating clustering for gene expression data[END_REF][START_REF] Ramaswamy | cancer diagnosis using tumor gene expression signatures[END_REF][START_REF] Steenman | Distinct, "molecular portraits of human failing hearts identi_ed by dedicated cDNA microarrays[END_REF][START_REF] Dudoit | A prediction-basedresampling method for estimating the number of clusters in a dataset[END_REF]. Nevertheless, this work is not adapted H.Ben Saber is with the Informatics Department, Hight School of Sciences and Technologie, University of Tunis, 5 Av Taha Hussein, Tunis 1089, Tunisie F.Harrathi is with the Informatics Department,ISSAIT, University 7 th November M.M.Gammoudi with the Informatics Department ISSAIT, University 7 th November to data of certain contexts to know disproportionately amongst attributes in front of the objects number where we find few columns only comparing lines, which is atypical in the field of data-mining [START_REF] Datta | Comparisons and validationof statistical clustering techniques for microarray gene expression data[END_REF].

III. RELATED WORK

Multiple methods of classification were proposed to analyze the data resulting from microarrays presenting some biological limits. Indeed, it proves that the discovery of the relations between all genes is partial [START_REF] Trang | Data mining of gene expression microarray via weighted pre_x trees[END_REF], moreover the associations between genes is not solved. From data processing point of view, for: -The statistical methods, although they are rather effective, certain works [START_REF] Alizadeh | Distinct types of 6di_use large B-cell lymphoma identi_ed by gene expression pro_ling[END_REF]2] show that the tests on a reduced number of samples do not give satisfactory results for the analysis of differentially expressed genes.

-The supervised classification methods, the method of K neared neighbors is not directly applicable to more than two classes [START_REF] Ramaswamy | cancer diagnosis using tumor gene expression signatures[END_REF]. Moreover, the centroid classification method is sensitive to the aberrant values and disturbed from the microarrays [START_REF] Yeung | Validating clustering for gene expression data[END_REF][START_REF] Ramaswamy | cancer diagnosis using tumor gene expression signatures[END_REF]. Then for the linear discriminating analysis method, it is not very effective and sensitive corresponding to data to treat [START_REF] Dudoit | A prediction-basedresampling method for estimating the number of clusters in a dataset[END_REF].

-The non supervised classification methods, the hierarchical classification method involves the main issue considering the interpretation of a graph as a delicate task [START_REF] Eisen | Cluster analysis and display of genome-wide expression patterns[END_REF]. Moreover, the choice of a distance from similarity is not permitted [START_REF] Steenman | Distinct, "molecular portraits of human failing hearts identi_ed by dedicated cDNA microarrays[END_REF]. To avoid the limits of classification methods, we propose a classification and extraction of knowledge's method using Formal Concepts Analysis.

IV. MATHEMATICAL FOUNDATION USED FOR THE CLASSIFICATION AND THE EXTRACTION OF KNOWLEDGE

A. Binary relation

A binary relation is a subset of the Cartesian product

F E  . An element of a relation R is noted ) , ( y x
, x is said an argument of R , and y is an image of x by R . We note y R x the fact that an element x of E is related to an element y of F [8].

An approach to the classification and Knowledge Extraction via microarrays by using Formal Concepts Analysis

' ' ' A A R B A B A       and ' B B  . For a context ) , , ( R F E the paire ) , ( B A such as F B E A   , is a concept if ) ( ) ( B f A f  and A B g  ) (
; with f and g are the Galois Connexion. A is called the extension; B is called the intension [8].

D. Galois lattice

Let R be a binary relation, T is the whole of extracted concepts from this relation.

) , (  T is a Galois lattice. In other words, the Galois lattice T for a binary relation R is the whole of all extracted concepts from this relation of partial order  [8].

E. Order relation

Let ) , ( 1 1 1 B A RE  and ) , ( 2 2 2 B A RE  be two concepts of R , ) , ( ) , ( 2 2 1 1 B A B A  2 1 A A   and 2 1 B B  .  is a partial relation[8].

F. Graph of concepts

Let R be a binary relation definite on the sets E and F . let G a set of optimal concepts obtained from the conceptual decomposition of R .

) , (  G forms a graph of concepts. The obtained concepts are a subset of covered Galois Lattice. If G would be the set of concepts, ) , (  G form a galois lattice [START_REF] Gammoudi | Décomposition conceptuelle des relations binaires et ses applications[END_REF].

G. Definitions related to associations rules

An association rule allows correlations between the elements of a formal context which are related. The rule

j i X X 
occurs in the set of transactions with a support S and a confidence C such as: The support is the percentage of transaction which contains

j i X X  , called P , where : S X X P X X Support j i j i     ) ( ) (
and confidence is the percentage of transaction containing i X which also contains j X , called P , where :

C X X P X X Confidence i j j i    ) / ( ) (
. The rules having the

form j X X  are known as relevant if Support MinSupp X X j i  ) (  and Confiance MinConf X X j i  ) ( 
such as MinConf and MinSupp are minimal thresholds. Association rules are said to be valid if the measurement of confidence is higher or equal to a minimal threshold, which is called MinConf. Two types of association rules are related to the measurement of confidence. Thus, a rule is called exact if

1 ) (   j i X X Confidence
, and a rule is known as

approximate if 1 ) (  j i X X Confiance  : [11, 9] V. CONTRIBUTION
In this section we present some suggestions to explain the process of coupling algorithms "DCRB -INC" and "AprioriT ID" method.

A. Structure of the method

Our method is composed of four stages (Fig 1), namely the matrix transposition, the discretization of the selected data, the extraction of concepts from the binary relation and the generation of association rules. Through this method, we initially compute many lines and columns, we use the transposed matrix in the objective to use a reduced number of columns compared to the number of lines [2]. We generate in the second time associations between genes or patients by using an algorithm of conceptual graph generation covering the associated binary relation. After that, we apply an algorithm generating the association's rules to these concepts.

a-Data organization

We present an example of microarrays extraction context.

G is the set of L genes (itemsets) presenting columns of the matrix called concept intension. C is the set of N conditions (objects) presenting line set called the concept extension. R is a binary relation between C and G . This step consists in extracting from knowledge of transposing matrix instead of making the extraction starting from original data to control proportionately between lines and columns.

b-Data selection

The selection includes the discretization which consists in transforming the numerical data in a binary matrix. The form of discretization is established after having selected the data input, by using intervals of discretization published by the expert. As the end result depends on the choice of the discretization threshold, we left this task to the expert specialized in the domain. To classify genes, we suppose that each column represents a gene and that each line gives the results of gene level experiment measures .The genes which are over-expressed; in other words those which have important biological activities. For our study, the genes which are over expressed in a biological situation are coded by the Boolean value 1, and the others, are coded by the Boolean value 0. The obtained matrix form a formal context of the genes G set and the conditions C set.

B. Classification and generation of the association rules prototype

This section presents the software components developed through a process of knowledge extraction applied to a case study which is Breast cancer. Through this prototype, we show the contributions of this practice on data containing a large number of attributes, like the data of genome. With this intention, we treat a data file rather large compound of 77 blades of 9216 genes of breast cancer [START_REF] Trang | Data mining of gene expression microarray via weighted pre_x trees[END_REF]. To meet this need, we implemented a prototype, called EC -ACF@DNArray.

a-Management of data

We use a form to be filled with the names of the input files (data file, file of the gene names, and files of the blade names), the values of the thresholds to separate the levels from kind of genes and the choice of the expression profile from genes or of the blades. Indeed, to manage the diagonal data of our case study, these data are standardized and cleaned for the study of the case of expression gene in two different cultures. We point out that the first culture (of reference) is marked with the green fluorochrome (CY 3 ) and that of the second culture (of test) is with the red fluorograme (CY 5 ): The level of gene expression is determined by the According to the obtained results, we confirm that we could minimize the concept number presented by the nodes compared to the other methods. At this level, we could reach our primary purpose which consists in generating the gene classes or the patient classes to have a certain gene characterization or patients. This characterization facilitates to the specialised in this field the making of certain decisions, such as the setting of the patients under the same treatment, or the application of certain analyzes for co-expressed gene units. To give more significations to the generated classes, we propose to release the whole association rules. It is the extraction of relevant knowledge starting from the microarrays.

d-The results visualization by generating the association rules

In our approach, we work on scanned data, quantified, and deferred on a matrix level gene differentially expressed. With this intention, we carry out gene regroupings and/or samples according to their expression profile. In figure ( We can also express the results obtained in the form of histogram (Fig 6). Moreover, the association rules generated by the gene profiles detect the whole of Co-on-expressed genes. For our case study, we notice that genes BRCA1, BRCA2, BRCATA, BRCA3, BWSCR1A, TP53 BRIP1, RB1CC1, RAD51, CHEK2, BARD1 are implied in the disease of breast cancer. Moreover, genes PPM1D, PIK3CA, AKT1, PALB2, CASP8, TGFB1, NQO1, HMMR, PTEN, STK11, ATM, NCOA3, ZNF217 are also implied in the disease but with a low support of confidence.

C. Results obtained

To test our method, we developed a prototype where we examined data resulting from the micro arrays. The input of our prototype is the rough raw data files and the choice of the thresholds applied to the matrix, namely the profile of expression of genes or profile of expression of blades. The prototype discretizes the data, selects the data in which the specialist is interested then extracted the concepts by applying algorithm DCRB-INC and generates the association rules. Indeed, a concept includes genes which are over-expressed in the same experiments. If the appearing experiments in the extension share same properties, so again, we can make assumptions on genes possible. Moreover, one gene (or an experiment) can appear in several concepts (in opposition to what occurs the" clustering in the case of"). If the biologist is interested in a particular gene, we can thus study according to the biological situations which are the genes related to this one. Currently, the results which we obtained seem to be promising. In fact, our classification enabled us to raise what follows : the obtained classes by our method and organized in the lattice form the interpretation of our results make easier by the biologist thanks to the lattice structure. These results appear important the fact that such an information can be used in the clinical medium.

VI. Conclusion

In this paper, we presented our method of knowledge generation considered to be relevant by the expert starting from the data resulting from the microarrays. The suggested method generates conceptual graph either for the blades (individuals, or experiments), or for genes, as it can generate the whole of the association rules. To achieve this purpose, we selected data to which we applied a thresholding to release the binary context, and association rules relative to this concept lattice. We evaluated our work by making tests on real data files provided by Pasteur institute of Lille 2. Moreover, the interpretation of the results seems to be more reliable compared with former works.

FUTURE WORKS

This work permits to consider several prospects. Indeed, it would be interesting to finalize this work on other biological data by exploiting other diseases and while referring large number of accessible knowledge on the Web by experts, annotated by a several properties and this by considering scaling. Moreover, we propose to cover the subject of classification under-constraints to take account of the constraint types with the level of individuals such as the constraint" must-link" (two individuals must be in the same class) and the constraint"cannot-link" (two individuals should not be placed in the same class). Finally, to improve the quality of analysis, while considering the semantic aspect between the objects, what could play a big role in the classification and classes modeling to improve quality of the partition in clusters objects.

Fig. 1 .

 1 Fig.1. General principle of EC-ACF@DNA-Array

  c-Graph of concept extraction We generate according to algorithm DCRB-INC, the graph of concepts which represents in our case sub lattice of concepts d-Generation of the association rules In the context of the microarrays extractions, we use the rules of association to identify groups of co-expressed genes or classifications of blades. Let conditions (objects). For our study case, a reason is an itemset presenting a subset of G. A reason G I  is included in an object C J  , if I and J are in relation by binary relation R . Let an itemset G I  , the support of itemset I in the formal context C

N

  experimental different measurement forming a matrix from profile of gene expression. We determine the level of gene expression change between two samples CY 3 and CY 5 . For example, a spot representing a gene in sample CY 5 must be controlled twice as much as same gene in sample CY 3, what amounts saying that the log 2 of the report is equal to 1 or to -1. Therefore, if 1 M ij  , gene i is called over regular under the condition j associated with culture CY 5 is called unregular under the condition j associated with culture CY 5 . If 1 M ij   , gene i is known as lower regular under the condition j associated with culture CY 5 . b-Discretization This stage consists in generating of the context binary contexts structure, by carrying out the choice of the thresholds according to the biological instructions. We present in the following figure a formal context (Fig 2) to treat patients.

Fig. 2 .

 2 Fig.2. Generated binary context It is the same to give of a matrix which refers to the genes and which represents the transposed matrix. c-Classification method of the microarrays data The EC-ACF@DNArray proposal is effective research tools genes and blood groups. During the development of our prototype DCRB-INC, we generated for a context (of 125 Objects and 25 attributes) 100 concepts (Fig 3) and 246 associations' rules. In our case study, a bi-cluster is a node in

Fig. 3 .

 3 Fig.3. Resulting binary context

  Fig 5), we present the association's rules generated by the blades profiles where each gene plays the part of an object and each sample represents an attribute of gene.

Fig. 4 .Fig. 5 .

 45 Fig.4. Hasse Diagram generated to classify blades .

Fig. 6 .

 6 Fig.6. Illustration of the rules


H ∞ Robustification Control of Existing Piezoelectric-Stack Actuated Nanomanipulators

Hamid Ladjal, Jean-Luc Hanus, Antoine Ferreira

Abstract-In current AFM-based nanomanipulation systems, the commercial position closed-loop controller for piezoelectric nanopositioning stages are implemented with success in a wide range of industrial applications. Even if these controllers operate with satisfactory nominal tracking performance, considerable attention has been focused on appropriate control strategies to compensate hysteresis, nonlinearities, drift and creep for high bandwidths and large scanning regimes. As these closedloop controllers are very cost-effective, a special interest in robust plug-in compensators seems to be a solution. We proposed in this paper a robust plug-in compensator using the H-infinity loop-shaping techniques which can be plugged into the existing controller without affecting the already satisfactory nominal tracking performance of the existing closed-loop system. Dynamic modeling, identification and robust control of a 3 d.o.f. piezoelectric nanorobotic positioner are presented in this paper in order to improve the nanorobot performance under plant parameter variations and in the presence of external disturbances. Simulation and experimental results are given to validate the proposed plug-in robust compensator in the case of a nanorobotic manipulation task.

I. INTRODUCTION

I N biotechnology field, the tele-nanomanipulation control of microobjects (biological cells, viruses, MEMS) with human handling is difficult due to dynamics and uncertainties which are strongly impacted by the user's gesture. In applications where high performance and accuracy are not critical constitutive nonlinearities of piezoelectric nanopositioning stages and hysteresis can be compensated by standard Proportional-Integral (PI) or Proportional-Integral-Derivative (PID) controllers, but this can potentially lead to bandwidth limitation and inefficiencies. Even if these controllers have proved their performance in nanoscale positioning systems, challenging problems of nanoscale control remain due to nonlinear dynamics, actuator's modeling uncertainties, instabilities and lack of robustness against external perturbations and sensor noise [START_REF] Ge | Tracking control of a piezoceramic actuator[END_REF]. As these industrial closed-loop controllers are cost-effective and dedicated, a special interest in robust plug-in compensators seems to be a solution. The key idea of this paper is to robustify existing controllers by "Plug-In" attachable robust compensators for piezoelectric nanopositioning systems. Principle of the new methods is to settle on real plants as their nominal models with local compensators. Since the additional local compensators are designed independent of previously designed controllers, they are applicable for any existing control systems including nonlinear and/or non-closed form control scheme.

Hamid Ladjal is with PRISME Institute, ENSI de Bourges, F-18000 Bourges, France.

Antoine Ferreira is with PRISME Institute, ENSI de Bourges, F-18000 Bourges, France(e.mail: antoine.ferreira@ensi-bourges.fr).

A brief analysis of robust control techniques shows that a considerable number of feedback design schemes based on linear robust control techniques have been already proposed for nanopositioning systems. Development of inversion-based feedforward control with robust feedback control have proved their efficacy in output tracking in Atomic Force Microscope [START_REF] Zou | Control issues in high-speed AFM for biological applications: Collagen imaging example[END_REF], [START_REF] Zou | Preview-based optimal inversion for output tracking: Application to scanning tunneling microscopy[END_REF]. The performance of the inverse feedforward control, however, is strongly limited by modeling errors/uncertainties and disturbances [START_REF] Devasia | Should model-based inverse input be used as feedforward under plant uncertainty[END_REF]. Closed-loop linear H ∞ control technique seems to be an efficient alternative technique. These schemes have provided improvements in bandwidth and robustness to non linearities and hysteresis. In [START_REF] Tsai | Robust Tracking control of a piezoactuator using a new approximative hysteresis model[END_REF], the authors proposed a Smith predictor-based H ∞ controller for a piezoactuator with an emphasis in reducing the hysteresis. In [START_REF] Salapaka | High bandwidth nano-positioner: a robust control approach[END_REF], a H ∞ controller design for one-dimensional nanopositioning system performing high closed-loop bandwidths and robustness against nonlinearities has been synthesized with success. Finally, a robust Glover-McFarlane H ∞ scheme [START_REF] Sebastian | H ∞ loop shaping design for nanopositioning[END_REF] to simultaneously achieve performance and robustness where neither specific tracking requirement nor a characterization of uncertainty are available a priori.

Considering these robust control schemes, we proposed a framework for increasing the robustness of existing industrial control schemes with a quantifiable compromise on performances. We proposed in this paper a robust plug-in compensator using the H ∞ loop-shaping techniques which can be plugged into the existing controller without affecting the already satisfactory nominal tracking performance of the existing closed-loop system. A similar approach has been introduced by Salapaka in [START_REF] Sebastian | Design Methodologies for Robust Nano-Positioning[END_REF].In the present study, the H ∞ loop-shaping technique is applied for robustification of complex industrial controllers achieving robustness with marginal reduction in performance. Dynamic modeling, identification and robust control of a 3 degree of freedom piezoelectric nanorobotic positioner are presented in this paper in order to improve the nanorobot performance under plant parameter variations and in the presence of external disturbances.

The paper is organized in the following way. In Section 2, a description of the nanorobotic device is given. This is followed by the frequency-domain-based system identification of the existing controller and nanorobotic positioner. The control design and the experimental results are then presented in Section 4. Experimental results are given in Section 5 to validate the proposed plug-in robust compensator in the case of a nanorobotic manipulation task.

II. PIEZOELECTRIC NANOMANIPULATOR DESCRIPTION

The nanomanipulator structure is composed of three linear translation stages (x,y,z) driven by DC motors for coarse motion (range: 8 mm, accuracy: 15 nm) combined with a 3 d.o.f ultra-high-resolution piezomanipulator (x,y,z) for fine positioning (range: 100 μm, accuracy: 1 nm). This hybrid nanopositioning system combines the advantages of ultra-low inertia, high-speed and long travel range (Fig. 1(a). The micro-endeffector is constituted by a piezoresistive AFM cantilever integrating a full-bridge strain gauge sensors ((Fig. 1(b)). Considering these limitations, we proposed in the following a robust plug-in compensator using the H ∞ loop-shaping technique which can be plugged into the existing PI-controller without affecting the already satisfactory nominal tracking performance.

III. CONTROLLER DESIGN AND IDENTIFICATION

A. Controller Design

In this paper, we simply assume that the controller has been already designed by the manufacturer and is in operation in a real micromanipulation system. However, it is often the case that this controller may not work well when the plant is perturbed and/or external disturbances presents. In this situation, an additional controller is needed to improve the robustness of the overall system against plant uncertainties. It is desirable that this additional controller can be plugged into the existing controller and without affecting the already satisfactory tracking performance. Considering these limitations, we proposed in the following a robust plug-in compensator using the H ∞ loop-shaping technique which can be plugged into the existing PI-controller without affecting the already satisfactory nominal tracking. This the reason why we called that controller a plug-in controller.

B. Identification

The digital system to be identified is constituted of both components: the piezoelectric nanopositioner and the PIcontroller embedded into the acquisition card. The acquisition board is composed of digital inputs (DAC) and outputs (ADC), the identification should be considered as a discrete system. Among the various modes of identification for numerical systems ( [START_REF] Longchamp | Commande numrique de systmes dynamiques[END_REF], [START_REF] Rivoire | Cours d'automatique Tome 3[END_REF]), we choose the least squares simple method. This choice is justified by the simplicity of the method during implementation, accuracy of the identified parameters and off-line parameters adjustment. The least square method is based on the determination of a vector of parameters so as to minimize a error vector ε:

J(ϑ ) = ε(k) 2 = (y(k) -φ (k)ϑ ) T (y(k) -φ (k)ϑ ) (1)
The vector parameters minimizing the criterion J(ϑ ), denoted ϑ (k), is called the estimated vector parameters defined by ( [START_REF] Rivoire | Cours d'automatique Tome 3[END_REF]):

ϑ (k) = (φ T (k)φ (k)) -1 φ T (k)y(k) (2)
where φ (k) is the observation matrix φ (k) of p-order. In to order to identify the wide dynamics bandwidth of the servonanopositioner, we chose a Pseudo-Random Pattern Generator signal (PRPG). It has been already proved reliable for a good identification in many dynamic systems. The scheme of Fig. 2 depicts the principle of the numerical identification. U(Z) represents the numerical input signal (PRPG input signal), Y (Z) the numerical output signal of the system. The system is provided with analog-to-digital and digital-to-analog converters which make it possible to pass from discrete towards the continuous and the continuous towards the discrete space. In Fig. 2, the terms G(S) and H(z) represent the transfer functions of the continuous system and the discrete system, respectively. The parameters provided by the identification represent the a i and the b j of the discrete transfer function which is given by the following expression:

' $ $ ' 3, $ 6WDJH 6HQVRU Q [ + ] 8 ] * V < ]
0 z 2 + b 1 z + b 2 z 3 + a 1 z 2 + a 2 z + a 3 ( 3
)
The identification was carried out for several operating points, ranging from 0 to 60μm with elementary steps of 10μm order and driving frequencies lower than the resonance frequencies. The discrete models have been identified with averaged parameter values for x, y and z axes, respectively

H x (z) = 0.00535z 2 + 0.0271z + 0.0507 z 3 -0.00557z 2 + 0.0196z -0.93 (4)
H y (z) = 0.003z 2 + 0.00741z + 0.0117 z 3 -0.00579z 2 + 0.025z -0.997

(5)

H z (z) = 0.000556z 2 + 0.00535z + 0.00761 z 3 -0.00566z 2 + 0.0222z -1 (6)
Then, we used the Matlab function d2c to convert the discrete functions to continuous functions. The results presented in Fig. 3 show the experimental frequency responses. These results show good tracking performances with negligible error. However, it should be noticed that the identified models vary strongly due to uncertain dynamics and are subject to modeling uncertainties.

IV. PLUG-IN ROBUST CONTROLLER DESIGN

The existing controller is a proportional integral (PI) and for the robustification of the existing controller we use the H ∞ loop shaping technique and how it can be integrated in a robust plug-in controller.

A. Coprime Factorization Approach

An approach was developed by ( [START_REF] Glover | Robust stabilization of normalized coprime factors: An explicit H ∞ solution[END_REF]) and ( [START_REF] Glover | Robust Stabilization of normalized coprime factor plant descriptions with H ∞ bounded uncertainly[END_REF]) starting from the concept of the coprime factorization of transfer matrix. This approach presents interesting properties and its implementation calls upon notions traditional of automatic control.

B. Robust Controller Design using Normalized Coprime Factor

We define the nominal model of the system to be controlled starting from the coprime factors on the left: G(s) = M(s) -1 N(s). Then a perturbed model is written (see Fig . 4). where G is a left coprime factorization (LCF) of G, and Δ m , Δ n are unknown and stable transfer functions representing the uncertainty. We can then define a family of models with the following expression :

G = ( M + Δ m ) -1 ( N + Δ n ) (7)
ξ ε = { G = ( M + Δ m ) -1 ( N + Δ n ) : (Δ m Δ n ) ∞ < ε} (8)
where ε max represent the margin of maximum stability. The robust problem of stability is thus to find the greatest value

) (s K ) 0 ( ) 0 ( 1 W K ⋅ référence 1 W ) (s G 2 W 1 W ) (s G 2 W (b) (a)
G a (s) + -Fig. 5. Loop Shaping design procedure.

of ε = ε max , such as all the models belonging to ξ ε can be stabilized by the same corrector K. The problem of robust stability H ∞ amounts finding γ min and K(s) stabilizing G(s) such as:

I K (I -GK) -1 (I W 2 GW 1 ) ∞ = γ -1 min = ε -1 max (9)
However, [START_REF] Mc Farlane | A Loop Shaping design procedure using H ∞ synthesis[END_REF] showed that the minimal value γ of is given by:

γ min = ε -1 max = 1 + λ sup (XY ) (10)
where λ sup indicates the greatest eigenvalue of XY . Moreover, for any value ε < ε max a corrector stabilizing all the models belonging to ξ ε is given by:

K(s) = B T X(sI -A + BB T X -γ 2 ZYC T C) -1 γ 2 ZYC T Z = (I +Y X -γ 2 I) -1 (11) γ = ε -1
where A, B and C are state matrices of the system defined by the function G and X, Y are the positive definite matrices and solution of the Ricatti equation.

C. The Loop Shaping Design Procedure

Contrary to the approach of Glover-Doyle, no weight function can be introduced into the problem. The adjustment of the performances is obtained by affecting an open modeling (Loop Shaping) process before calculating the compensator. The design procedure is as follows :

• We add to the matrix G(s) of the system to be controlled a pre-compensator W 1 and/or a post-compensator W 2 . The nominal plant G(s) and shaping functions W 1 and W 2 are combined in order to improve the performances of the system such as G a (s) = W 2 (s)G(s)W 1 (s) (see Fig. 5.a). • From coprime factorizations of G a (s), we apply the previous results to calculate ε max , and then synthesize a stabilizing controller K ensuring a value of ε slightly lower than ε max .

• The final feedback controller is obtained by combining the H ∞ controller K with the shaping functions W 1 and

W 2 such that G a (s) = W 2 (s)G(s)W 1 (s) (See Fig.5.b).
V. CHARACTERIZATION OF PIEZOELECTRIC-STACK ACTUATED NANOMANIPULATOR

A. Implementation of the controllers

The controller K is obtained by combining the pre-filter and the post-filter . The pre-filter and post-filter are used to shape the open-loop plant to achieve a desired frequency responses according to some well defined design specifications such as bandwidth and steady-state error ( [START_REF] Lundstrom | Weight selection for H ∞ and mu-control methods: insights and examples from process control[END_REF]). To obtain a high gain at low frequency, a PI controller is synthesized for the x-axis like pre-filter W 1 and in order to obtain a small gain at high frequency, a low-pass filter is synthesized like pre-filter W 2 . In order to obtain a high performance and a good robustness, we add the following weight functions:

W 1 = 50 × 10s+35.2 10s ,W 2 = 1 s+11.3 .
Using this weight functions we obtain six-order the H ∞ controller. In order to implement this controller we use the function of Matlab balmr to reduce the order and obtain a third-order approximation of the plug-in robust controller given by : K r (s) = -72.69s 2 -1824s -6706 s 3 + 66.88s 2 + 625.7s -5.07 × 10 -11 .

(

) 12
and the discrete controller :

K r d (s) = 0.1726z 2 -0.3347z + 0.1622 z 3 -2.842z 2 + 2.688z -0.846 (13)
The Fig. 6 and Fig. 7

B. Nanomanipulator Characterization

In this section, the nanopositioning is characterized in terms of range, sensitivity and resolution in the open-and closed-loop configurations. The calibration data showed some hysteresis in open-loop. Hysteresis is primarily due to the nonlinear relationship between applied voltage and displacement which are important for large deflections. To present the effectiveness of the H ∞ closed-loop design, the hysteresis curves obtained in open-loop are compared with the closedloop design using PI controller and the H ∞ controller. We can see clearly that for a nanopositioning displacement of 40μm, a maximum output hysteresis of 10μm (25%) was observed. The same experiment with the closed-loop controller showed that the effects were practically eliminated (Fig. 9). An important observation comes from the fact that the operating applied voltage is settled to 100V while the system was identified for a low value of input signal (PRPG signal). It shows clearly the linearity relationship between the input-output signals. Similar linearity results have been measured when considering creep effects. It should be noticed that important differences were observed for different displacement with PI and H ∞ controller. High precision is obtained with robustified PI controller.

C. Robust Experimentation in Microsphere pushing

To illustrate the behavior of the piezoelectric-stack actuated nanomanipulator system and to confirm the findings of the previous section. The working space of the nanomanipulator is visualized through two orthogonal optical microscopes (topview and side-view) and a high magnification microscope for a zooming view (see Fig. 10). The micromanipulation task is to realize micro pushing tasks of microparticles of 5 μm of diameter on a glass substrate. During the micro pushing tasks, we made a comparison between the responses with PI and H ∞ controllers when pushing a microparticle. We give different responses of the system with plug-in robust controller in the presence of the disturbance and position sensor noise. The set of experiments of Fig. 11 shows the position with the plugin robust controller. The experimental results demonstrate the excellent position tracking with the proposed plug-in robust controller and the output track the reference and we can see clearly the better tracking PI controller (Fig. 11.b). The experimental response with a step perturbation demonstrate the excellent tracking and the position follows well the reference and not affected by the application of a disturbance (Fig. 12.a). To emphasize the noise problem of the closed loop system with the plug-in robust, Fig. 12.b demonstrate the influence of this parameter on the performance of the robust controller. as perturbations to the nominal model, and robustness is guaranteed by ensuring that the stability specifications (nanometer resolution) are satisfied in the worst-case uncertainty. As conclusion, a very good point of the design is the remarkable features that achieves robustness with marginal reduction of performance. Our future work will consist of using the piezoelectric-stack actuated nanomanipulators in robust nanoteleoperation of mouse embryonic stems cells (mESCs) for regenerative cell construction.

Fig. 1 .

 1 Fig. 1. Structure of the hybrid 6 dof AFM-based nanomanipulator and its sensorized force cantilever. The 3 d.o.f ultra-high-resolution piezomanipulator (P-611.3S NanoCube from Physics Instruments) is a versatile, multi-axis piezo-nanopositioning system . Its 100 x 100 x 100 μm positioning and scanning range comes in an extremely compact package. Equipped with a zero-stiction, zero-friction guiding system, this piezomanipulator provides motion with ultra-high resolution and settling times of only a few milliseconds. Single-axis nanopositioning stage with anti-accurate-motion flexure design. The best flexure designs provide guiding precision in the low nanometer range. In open-loop operation, the platform's position is roughly proportional to the drive voltage. In the closed-loop version, the Proportional-Integral controller allows absolute position control. However, existing closed-loop controller is designed to achieve specific tracking and performance requirements (such as zero steady-state tracking) in a narrow frequency closed-loop bandwidth. During nanomanipulation tasks, the xy-z nanopositioner servo-controller lacks in robustness against: a) Modeling uncertainties: 1) Mechanical nonlinearities of mechanical-guiding systems for large travel ranges; 2) Modeling uncertainties due to operating point, temperature effects and time execution; 3) Hysteresis and creep effects du to piezoelectric ceramics. b) External perturbations: 1) Noise measurement perturbations; 2) Force/torque perturbations during nanomanipulation. Considering these limitations, we proposed in the following a robust plug-in compensator using the H ∞ loop-shaping technique which can be plugged into the existing PI-controller without affecting the already satisfactory nominal tracking performance.

Fig. 3 .

 3 Fig. 3. Experimental frequency responses of (a) x-axis, (b) y-axis and (c) z-axis.

Fig. 2 .

 2 Fig. 2. Identification system of the transfer function.

  b

Fig. 4 .

 4 Fig. 4. Coprime factor robust stabilization problem.

  show the frequency responses of the weight functions W 1 , W 2 , of the open-loop system W 1 .P.W 2 and of the controller K . The results show that the open-loop remains close to the step response obtained after the choice of the shaping functions and K r ensures correct margins of stability.

Fig. 6 .

 6 Fig. 6. Bode diagram of W 1 and W 2 .

Fig. 7 .

 7 Fig. 7. Bode diagram of open-loop system W 1 .P.W 2 and K r .

Fig. 9 .

 9 Fig. 9. PI and H ∞ controller to remove the hysteresis and the zoom of the response.

Fig. 10 .

 10 Fig. 10. Different optical microscope views (top view, side view and zoom view) of the nanomanipulator pushing microparticles of 5 μm of diameter.

Fig. 11 .

 11 Fig. 11. (a) Step responses for PI and H ∞ controllers and (b) corresponding static errors.

Fig. 12 .Fig. 8 .

 128 Fig. 8. PRPG responses with (a) PI controller (b) robustified PI controller.

Fig. 12 .

 12 Fig. 12. Step responses with robustified PI controller in presence of (a) step perturbation (b) white noise.


Hamid Ladjal

Jean-Luc Hanus

Anand Pillarisetti

Carol Keefer

Antoine Ferreira

Jaydev P Desai

Realistic Visual and Haptic Feedback Simulator for Real-time Cell Indentation

Comprehensive and training simulators that provide realistic visual and haptic feedback during cell indentation tasks are currently investigated. Complex cell geometry inherent to biological cells and intricate mechanical properties drive the need for precise mechanical and numerical modeling to assure accurate cell deformation and force calculations. Advances in alternative finite element formulation, such as mass-tensor approach, have reached the state where they are applicable to model soft cell deformation in real time. The geometrical characteristics and the mechanical properties for different cells are determined with AFM indentation. A real-time, hapticsenabled simulator for cell centered indentation has been developed which utilizes the atomic force microscopy data (mechanical and geometrical properties of embryonic stem cells (mESC)) to accurately replicate the indentation task and predict the cell deformation during indentation in real-time. This tool can be used as a mechanical marker to characterize the biological state of the cell. The operator is able to feel in real-time the change in the stiffness during cell deformation between fixed and live cells. A comparison study with finite element simulations using a commercial software and the experimental data demonstrate the effectiveness of the proposed physically-based model.

I. INTRODUCTION

Mechanical manipulation and characterization of biological cells is currently one of the promising research areas in the field of medical robotics applied to cellular level interactions. The mechanical properties such as elasticity, membrane tension, cell shape, and adhesion strength may play an important role in the cell fate and differentiation [START_REF] Settleman | Tension precedes commitment-even for a stem cell[END_REF]. Embryonic stem cells (mESCs) can give rise to every cell type in the body unlike adult stem cells and thus have huge therapeutic potential. Stem cell based systems offer a very promising and innovative alternative for obtaining large number of cells for early efficacy and higher toxicity screening. The mechanical properties of biological cells have been studied with different techniques [START_REF] Desai | Engineering Approaches to Biomanipulation[END_REF], the most popular are optical tweezers [START_REF] Arai | Micromanipulation based on micro physics, strategy based on attractive force reduction and stress measurement[END_REF], magnetic beads [START_REF] Bausch | Measurement of local viscoelasticity and forces in living cells by magnetic tweezers[END_REF], and micropipette aspiration [START_REF] Ohasti | Intracellular mechanics and mechanotransduction associated with chondrocyte deformation during pipette aspiration[END_REF]. AFM has been widely used in the study of micro-and nanostructures including living cells [START_REF] Alonso | Feeling the forces: atmoic force microscopy in cell biology[END_REF]. Modern hamid.ladjal, jean-luc.hanus, antoine.ferreira, @ensi-bourges.fr A. Pillarisetti is with Department of Mechanical Engineering, University of Pennsylvania, USA anandpi@seas.upenn.edu J.P. Desai is with the Department of Mechanical Engineering, University of Maryland, College Park, MD, USA jaydev@umd.edu C. Keefer is with the Department of Animal and Avian Sciences, University of Maryland, College Park, MD, USA ckeefer@umd.edu AFM techniques allow solving a number of problems of cell biomechanics by simultaneous evaluation of the local mechanical properties and the topography of the living cells, at a high spatial resolution and force sensitivity [START_REF] Touhami | Nanoscale mapping of the elasticity of microbial cells by atomic force microscopy[END_REF]. In these experiments AFM cantilever serves as a microindenter to probe the cell, and further analysis of force-indentation data yields the local Young's modulus. In addition, AFM indentation technique can be used to characterize the viscoelastic behavior of the cell cytoskeleton [START_REF] Darling | Viscoelastic properties of zonal articular chondrocytes measured by atomic force microscopy[END_REF], including viscosity [START_REF] Mahaffy | Quantitative analysis of the viscoelastic properties of thin regions of fibroblasts using atomic force microscopy[END_REF], loss and storage moduli [START_REF] Alcaraz | Microrheology of human lung epithelial cells measured by atomic force microscopy[END_REF], and stress relaxation times. AFM-based indentation techniques with force feedback have been applied here to characterize single-cell mouse embryonic stem cells (mESCs) [START_REF] Pillarisetti | Mechanical Response of embryonic stem cells using haptics-enabled atomic force microscopy[END_REF]. We developed finite element models that include the topological information of adherent stem cells (shape and dimensions) and AFM tips (conical, spherical), and biological structure (cells models with cytoplasm layers, cytoskeletons and nucleus) [START_REF] Lulevich | Cell Mechanics using Atomic Force Microscopy-based Single-Cell Compression[END_REF]. The paper consists of five sections. In Section II, we present the bio-mechanical finite element model in real-time indentation. In Section III, we present the methodology of the virtual environment system for cell indentation. In Section IV, we present the experimental setup using the force microscopy. In Section V, we make a comparative study between the physics-based model and experimental data. In Section VI, we present our real-time cell indentation simulator using a spherical tip. Finally, we give some concluding remarks and the directions for future work.

II. PHYSICALLY BASED MODELING OF DEFORMABLE OBJECTS FOR CELL NANOINDENTATION SIMULATION

In this section we present different steps for the linear elastic mass-tensor model introduced by [START_REF] Bro-Nielsen | Real-Time Volumetric Deformable Models for Surgery Simulation Using Finite Elements and Condensation[END_REF] and extended by the authors in [START_REF] Schwartz | Modeling Liver Tissue Properties Using a Non-linear Visco-Elastic Model for Surgery Simulation[END_REF]. In our work, the cell is modeled by the volume object discretized into a conformal tetrahedral mesh as defined by finite element theory. Using classical notations, inside each tetrahedron T k , the displacement field, defined by a linear interpolation N k of the nodal displacement vector u k of the four vertices of tetrahedron, is written as follows:

U(x) k = N k (x) u k (1)
For a linear elastic material, the relation between the Cauchy stress tensor and the linearized strain tensor is written with Lame's coefficient in condensed vector notation as : where I is the identity matrix.

σ σ σ = λ ε ε ε 1 + ε ε ε 2 + ε ε ε 3 I +2µ ε ε ε (2)
The principle of virtual work applied to a single tetrahadron T k leads to the elementary stiffness matrix K k such that the elementary nodal force vector acting on a tetrahedron is:

f k = K k u k (3)
This stiffness matrix is composed of a plurality of elementary submatrices each connecting the elementary force acting on the node, i, to the displacement of the node, j:

K k ij = 1 36 V k λ m i m j T + (4)
µ m j m i T + µ m i T m j I
where m are unit outward-pointing normals to triangular faces and V k is the volume of the tetrahedron T k . Taking into account the contribution of all adjacent tetrahedra, the global internal force acting on a node l can be expressed as follows:

F l int = k∈V l 4 j=1 K k ij u j (5)
where V l is the neighborhood of vertex l (i.e. the tetrahedra containing node l).

The tensors K k ij , depending on the rest geometry and Lame's coefficients, are constant. They can be pre-computed in an off-line phase. It's the essential advantage of the mass-tensor approach which makes it useful for real-time application.

A. Cell modeling using finite element model

In our work, the mESC is meshed with 3D first-order tetrahedral elements as shown in Fig. 1. The mesh that has been used for the following simulations is constituted by a non regular volume mesh composed of 835 vertices and 2456 tetrahedra where all vertices are free to move except the ones in contact with the petridish which are fixed.

B. Dynamic Model

The equation of motion of a vertex, l, of the cell mesh can be written:

M l ül + γ l ul + k∈V l 4 j=1 K k ij u j = F l ext (6)
where M l and γ l are respectively the mass and damping coefficients of each vertex. To solve the dynamic system, we tested different integration schemes (implicit and explicit) taking into account the tradeoff between real-time simulation and haptic stability requirements. We choose the explicit centered finite-difference scheme defined by: Fig. 2 shows the architecture of the real-time virtual reality based cell indentation simulator. A significant difficulty of using the finite element technique for real-time simulation is that it is computationally costly. These mesh-based schemes also require an expensive numerical integration operation for the computation of the system stiffness matrices. For this reason we adopted an off-line computation of entries of stiffness matrix for each triangular element and a real-time visual and haptic interaction for indentation simulation.

u(t + h)=u(t)+ u(t) h + 1 2 ü(t) h 2 + O(h 3 ) u(t + h)= u(t)+ 1 2 ü(t) h + 1 2 ü(t + h) h + O(h 3 ) ( 7

A. Off-line computation

All the operations costly in time are realized during this step. The mechanical and geometrical properties of the cell structure are determined through AFM indentation experiments and validated using analytical Hertz contact modeling. The off-line calculation of all tensors K ij allows precalculation stiffness values for a wide range of cell deformations.

B. Real-time Indentation

During the indentation task, the contact between the tip and the cell must occur at a special set of points (nodal points).

= k (R -d) n i ← i +1 j ← j +1

1) Collision Detection:

The collision detection is the first step to carry out realistic interactions. However, locating the contact primitive (e.g. facets) between two objects may be computationally expensive, especially if the objects are composed of a large number of polygons. In general, to determine a contact between two virtual objects we compute the distance between them. If this distance is negative then the objects are in contact. For the purposes of detecting collisions, we used a point-based representation of the cantilever tip as a simple spherical object and employed a simple spherepoints intersection algorithm with local search technique (See Fig. 3). It calculates the distance between the center of the spherical tip and the nodes of the cell surface (see Fig. 4). Fig. 4. Scheme illustrating the multiple collision detection algorithm using a simple sphere-points intersection problem: the distance between the center of the spherical tip and the nodes of the cell surface is calculated in real-time.

2) Penalty based models: In an interaction model, a collision between two objects may generate local and global deformations of the colliding objects. In our simulations we used the penalty based models. This model considers the reactive force that first reduces the relative speed of the two objects and finally repulses them from each other, as proportional to the deformation. Thus the reactive force is a function of the local interpenetration depth δ. The normal repulsive force F rep , can be written:

F rep = λδ n ( 8
)
where λ is the penalty of the interpenetration and n the outward-pointing normal.

3) Virtual coupling for stability of the haptic rendering: The internal operating loop of the haptic interface requires an update frequency neighboring of 1KHz. However, the finite element model update frequency between the 25 and 30Hz. This frequency difference does not ensure the coherence between both systems leading to instabilities of the user haptic rendering. We adopted the solution to use a virtual coupling model defined in [START_REF] Adams | Stability and performance of haptic displays : theory and experiments[END_REF].

IV. EXPERIMENTAL SETUP

A. Description

The experimental tests of the cells indentations were performed using an atomic force microscope (Model: MFP-3D-BIO, Asylum Research, Santa Barbara, CA) Fig. 5. The top view module enables viewing of cells and easy alignment of the laser beam on the AFM cantilever. The XY stage (manual) allows the user to position the cell beneath the cantilever tip of AFM. The entire AFM set up is enclosed in an acoustic isolation chamber to prevent acoustic noise from interfering with the AFM measurements. The xand y-axes range of the scan head is 90µm. The z-axis scan range is 40µm. When Fig. 5. Atomic Force Microscopy: The experimental set up for indentation studies and spherical geometrical Tip used for(mESC): The cantilevers attached with a spherical probe (5µm in diameter), the silicon nitride cantilever (Novascan Technologies, Inc.) the tip of the cantilever is in contact with the cell, the initial cantilever deflection, (d 0 ), and initial cantilever movement in (z) direction (z 0 ) are stored. As the cantilever moves in the z direction and deforms the cell, the final cantilever deflection, (d 1 ), and the cantilever movement, (z 1 ) are obtained. In the following experiments, we employed two types of cantilevers attached with a spherical probe (5 µm in diameter):

• Silicon nitride cantilever with a spring constant of 0.06 N/m (Novascan Technologies, Inc., Ames,IA) for live cells.

• Silicon cantilever with a spring constant of 1.75 N/m (Novascan Technologies, Inc) for fixed cells.

V. M OUSE EMBRYONIC STEM CELLS CHARACTERIZATION

A. Geometry of the cell and cantilever interaction with the surface

The cells were 10 and 15 µm in diameter using the cantilever tip to measure the interaction force between the cell and the tip. The geometric modeling of the cell requires the determination of the height, adhesion surface and cell's contour (see Fig. 6). The mESC height was calculated from the force-displacement curve obtained on: hard surface and on mESC surface. First, initial height h 0 is measured by bringing the AFM tip in contact with the hard surface near mESC and then the second height h 1 by bringing the AFM tip on mESC surface. The mESC height is the difference between h = h 0h 1 . Second, to determine mESC shape, AFM imaging experiments are carried out in order to determine mESC's contour and diameter using the image processing techniques [START_REF] Ladjal | Atomic Force Microscopy-Based Single-Cell Indentation: Experimentation and Finite Element Simulation[END_REF]. From Fig. 6 we selected mESC shapes with

B. Comparison between experimental data and analytical models

To estimate the mechanical properties of biological cells using the AFM, various analytical models can be used to identify Young's modulus of mESC in live as well as fixed cells. We performed single indentation studies with a spherical tip on mESC. The main goal is to determine whether it appropriately predicts the force-indentation relationship of mESC [START_REF] Ladjal | Atomic Force Microscopy-Based Single-Cell Indentation: Experimentation and Finite Element Simulation[END_REF]. The figure Fig. 7 shows the force versus cell indentation for live and fixed undiff mESC using the Hertz model. The relationship between the indentation δ, and the loading force F is given by:

F = 4 3 E * R 1 2 δ 3 2 (9)
where E * , R are the combined modulus and the relative curvature of the tip and the cell respectively

1 E * = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 (10) 1 R = 1 R 1 + 1 R 2 (11)
where (E 1 , ν 1 , R 1 ) and (E 2 , ν 2 , R 2 ) represent the elastic modulus, Poisson's ratio, and the radius of the cell and spherical indenter respectively. The elastic modulus of the silicon nitride cantilever and silicon cantilever were 222.22 GPa and 168.17 GPa used for the live and fixed cells respectively. The elastic modulus of the cells is in range of KPa. Hence, our assumption that the tip used for probing is infinitely stiff compared to the cell and Hertz contact model is valid. Thus, Eq. 9 can be rewritten as:

F = 4E 1 3(1 -ν 2 1 ) R 1 2 1 δ 3 2 (12)

VI. FINITE ELEMENT MODEL (FEM) SIMULATION OF CELL INDENTATION

The Fig. 8 presents a photo of visual and haptic user interface setup. To test the accuracy and reliability of the proposed user interface system with haptics enabled simulation, a set of experiments is designed using the setup given in Section 2. The geometrical dimensions of the cell are determined through AFM image processing and force curve deflection (diameter and height) and the mechanical properties through AFM indentation studies thanks to the use of an analytical Hertz contact model.

In order to demonstrate the validity of our finite element indentation model in real-time, we compared the simulation indentation using the physically-based FEM model to the incremental data provided by AFM indentation tool and to the response simulated with a commercial FE software. The figures (Fig. 9, Fig. 10) present the nonlinear relationship between the penetration distance (indentation) and the reaction force for live and fixed undiff cells. The finite element simulations using the commercial software show good agreement with both Hertz contact model and experimental indentation data. The rendered reflection force is represented by a cone VII. CONCLUSIONS We have developed a computer-based training system to simulate real-time cell indentation procedures in virtual environments for training biologists residents. The simulation of this procedure involves the real-time rendering of computer generated graphical images of AFM tip, physically-based modeling of soft biological tissue, and display of touch and force sensations to the training through the simulation of haptic interactions. We first investigated the challenging issues in the real-time modeling of the bio-mechanical properties of the cell indentation through finite element models. Compared to experimental AFM indentation results performed on mESC. Currently, we are working on integrating nanoscale effects such as friction, viscosity, punction. All modalities will be merged in an ergonomic tool and intelligent biological simulator for stem cell characterization.

H

  . Ladjal, J-L Hanus and A. Ferreira are with In-

Fig. 1 .

 1 Fig. 1. Representation of the 3D mESC mesh using P1 tetrahedral finite elements.

Fig. 2 .

 2 Fig. 2. Computational architecture for simulating force-reflecting deformable cell indentation in virtual environment. The figure shows the two simulation phases used for the real-time indentation of the cell: (i) off-line precalculations of stiffness matrices and (ii) simulation of visual and haptic interaction.

  START New tip position M = M + △M loop over triangles loop over vertices END i<n surf j<2 i ← 0 j ← 0 Compute distance between current node (P ) and tip (M ) d = MP False False False d<R Compute penalty force Fp

Fig. 3 .

 3 Fig. 3. Multi collision detection flowchart between a spherical tip and a stem cell: M represents the position of the center of the spherical tip and R its radius, {△M } represents the step of the imposed displacement, d is the computed distance between the center of the tip and a node of the external cell surface, Frep represents the repulsive force and {n} the outward-pointing normal

Fig. 6 .

 6 Fig. 6. How to determine the diameter of the stem cell through image analysis (ImagePro software from Asylum). circular adhesion surfaces. As illustration, the left corner inset shows a typical circular mESC where D 1 and D 2 stands for the horizontal cells diameters (D 1 ≃ D 2 close to 12µm).

Fig. 7 .

 7 Fig. 7. Force versus indentation for live and fixed mESC.

Fig. 8 .

 8 Fig. 8. Real-time indentation simulator using the Omni haptic interface.

Fig. 9 .

 9 Fig. 9. Force versus indentation for live mESC: Comparison between the real-time FEM, the Hertz contact, the quasi-static FEM and experimental data .

Fig. 10 .

 10 Fig. 10. Force versus indentation for fixed mESC: Comparison between the real-time FEM, the Hertz contact, the quasi-static FEM and experimental data .

Fig. 11 .Fig. 12 .Fig. 13 .

 111213 Fig. 11. Depth deformation of the cell; the left picture presents the depth indentation (2µm) of the cell, the right picture show the zoom of the cell deformation .

TABLE I

 I

	MECHANICAL AND GEOMETRICAL PROPERTIES USED IN THE
	HAPTICS-ENABLED SIMULATOR
		Live mESC Fixed mESC
	Young modulus (KPa)	0.169	26.7
	Poisson coefficient	0.49	0.49
	Diameter(µm)	≃ 12	≃ 12
	Height (µm)	8	8

ACKNOWLEDGMENT

We would like to acknowledge the support of National Science Foundation grant: CMMI 0826158, and STMD-Maryland Technology Development Corporation grant: 08071517, for part of this work. We would like acknowledge the support of the Centre de Recherches biologiques (CERB) de Baugy (France) for funding support.