Francisco De A T De Carvalho

Renata M C R De Souza

Marie Chavent
email: marie.chavent@math.u-bordeaux.fr

Yves Lechevallier
email: yves.lechevallier@inria.fr

Adaptive Hausdorff Distances and Dynamic Clustering of Symbolic Interval Data

Keywords: Symbolic data analysis, dynamic clustering, interval data, Hausdorff distance, adaptive distances

This paper presents a partitional dynamic clustering method for interval data based on adaptive Hausdorff distances. Dynamic clustering algorithms are iterative two-step relocation algorithms involving the construction of the clusters at each iteration and the identification of a suitable representation or prototype (means, axes, probability laws, groups of elements, etc.) for each cluster by locally optimizing an adequacy criterion that measures the fitting between the clusters and their corresponding representatives. In this paper, each pattern is represented by a vector of intervals. Adaptive Hausdorff distances are the measures used to compare two interval vectors. Adaptive distances at each iteration change for each cluster according to its intra-class structure. The advantage of these adaptive distances is that the clustering algorithm is able to recognize clusters of different shapes and sizes. To evaluate this method, experiments with real and synthetic interval data sets were performed. The evaluation is based on an external cluster validity index (corrected Rand index) in a framework of a Monte Carlo experiment with 100 replications. These experiments showed the usefulness of the proposed method.

Introduction

Cluster analysis seeks to organize a set of items (usually represented as a vector of quantitative values in a multidimensional space) into clusters such that items within a given cluster have a high degree of similarity, whereas items belonging to different clusters have a high degree of dissimilarity [START_REF] Bock | Classification and Clustering: Problems for the Future[END_REF], [START_REF] Jain | Data Clustering: A Review[END_REF]). Cluster analysis techniques can be divided into hierarchical and partitional methods [START_REF] Spaeth | Cluster analysis algorithms[END_REF], [START_REF] Gordon | Classification. Chapman and Hall/CRC[END_REF], [START_REF] Everitt | Cluster Analysis[END_REF]).

Hierarchical methods yield complete hierarchy, i.e., a nested sequence of partitions of the input data. Hierarchical methods can be agglomerative or divisive. Agglomerative methods yield a sequence of nested partitions starting with trivial clustering in which each item is in a unique cluster and ending with trivial clustering in which all items are in the same cluster. A divisive method starts with all items in a single cluster and performs a splitting procedure until a stopping criterion is met (usually upon obtaining a partition of singleton clusters).

Partitional methods seek to obtain a single partition of the input data into a fixed number of clusters. They usually produce clusters by (locally) optimizing an adequacy criterion. To improve cluster quality, the algorithm is run multiple times with different starting points and the best configuration obtained from the total runs is used as the output clustering. This paper addresses the partitioning of interval data often present in real applications. Table 1 shows an example of an interval data table. This kind of data have been studied mainly in Symbolic Data Analysis (SDA), a new domain in the area of knowledge discovery and data management related to multivariate analysis, pattern recognition and artificial intelligence. The aim of SDA is to provide suitable methods (clustering, factorial techniques, decision trees, etc.) for managing aggregated data described by multi-valued variables, where the cells of the data table contain sets of categories, intervals, or weight (probability) distributions [START_REF] Diday | The symbolic approach in clustering and related methods of data analysis[END_REF], [START_REF] Bock | Analysis of Symbolic Data, Exploratory methods for extracting statistical information from complex data[END_REF], [START_REF] Billard | From the statistics of data to the statistics of knowledge: Symbolic Data Analysis[END_REF]).

SDA has provided partitioning methods for clustering symbolic data. [START_REF] Diday | Symbolic Cluster Analysis[END_REF] proposed a clustering approach based on a transfer algorithm. [START_REF] El-Sonbaty | Fuzzy Clustering for Symbolic Data[END_REF] proposed a fuzzy k-means algorithm for clustering different types of symbolic data. [START_REF] Verde | A Dynamical Clustering Algorithm for symbolic data[END_REF] introduced a dynamic cluster algorithm for symbolic data considering context-dependent proximity functions. [START_REF] Gordon | An Iteractive Relocation Algorithm for Classifying Symbolic Data[END_REF] presented an iterative relocation algorithm that minimizes the sum of the description potentials of the clusters. [START_REF] Bock | Clustering algorithms and kohonen maps for symbolic data[END_REF] gives clustering strategies based on a clustering criterion and presents a sequential clustering and updating strategy for constructing a Self-Organizing Map in order to visualize symbolic interval-type data. [START_REF] Chavent | Dynamical Clustering Algorithm of Interval Data: Optimization of an Adequacy Criterion Based on Hausdorff Distance[END_REF] proposed a dynamical clustering algorithm for interval data where the prototype is defined by the optimization of an adequacy criterion based on the Hausdorff distance. Moreover, in [START_REF] Souza | Clustering of interval data based on city-block distances[END_REF], an adaptive dynamic clustering algorithm is presented for interval data based on City-block distances.

The main contribution of this paper is the proposal of a new partitional dynamic clustering method for interval data based on the use of an adaptive Hausdorff distance at each iteration.

The partitioning dynamical cluster algorithms [START_REF] Diday | La méthode des Nuées dynamiques[END_REF]) are iterative two-step relocation algorithms involving the construction of clusters at each iteration and the identification of a suitable representation or prototype (means, axes, probability laws, groups of elements, etc.) for each cluster by locally optimizing an adequacy criterion between the clusters and their corresponding representations [START_REF] Diday | Clustering Analysis[END_REF]). An allocation step is performed to assign individuals to classes according to their proximity to the prototypes. This is followed by a representation step where the prototypes are updated according to the assignment of the individuals in the allocation step, until the convergence of the algorithm, when the adequacy criterion reaches a stationary value.

The idea of dynamical clustering with adaptive distances [START_REF] Govaert | Classification automatique et distances adaptatives[END_REF], [START_REF] Diday | Classification Automatique avec Distances Adaptatives[END_REF]) is to associate a distance to each cluster, which is defined according to its intra-class structure. The advantage of this approach is that the clustering algorithm recognizes different shapes and sizes of clusters. In this paper, the adaptive distance is a weighted sum of Hausdorff distances. Explicit formulas for the optimum class prototype, as well as for the the weights of the adaptive distances, are found. When used for dynamic clustering of interval data, these prototypes and weights ensure that the clustering criterion decreases at each iteration.

In this paper, we present a dynamic clustering method with adaptive Hausdorff distances for partitioning a set of interval data. This method is an extension of the dynamic clustering algorithm based on non-adaptive Hausdorff distances proposed in [START_REF] Chavent | Dynamical Clustering Algorithm of Interval Data: Optimization of an Adequacy Criterion Based on Hausdorff Distance[END_REF]. In Section 2, a description is given of the classical dynamic clustering method with adaptive distances. This is followed by the presentation of the dynamic clustering method based on adaptive Hausdorf distances for interval data (Section 3). To validate this new method, Section 4 presents experiments with real and synthetic interval data sets. In Section 5, the concluding remarks are given.

Introduction to partitional dynamic clustering with adaptive distances

Let Ω be a set n of objects indexed by i and described by p variables indexed by j. Each object i is represented by a vector of feature values x i = (x 1 i , . . . , x p i ). Throughout this paper, we consider the problem of clustering Ω into K disjoint clusters C 1 , ..., C K such that the resulting partition P = (C 1 , ..., C K ) is optimum with respect to a given clustering criteria.

In dynamic clustering [START_REF] Diday | Clustering Analysis[END_REF]), we represent each cluster C k ∈ P by a prototype y k , which is also a vector of feature values. We measure the quality of this cluster by the sum of the dissimilarities d(x i , y k ) between objects i ∈ C k and the prototype y k . This measure of quality i∈C k d(x i , y k ) is called the adequacy criterion of the cluster C k . The classification problem is to find a partition P and a set L of K prototypes that minimize the following clustering criterion:

∆(P, L) = K i=1 i∈C K d(x i , y k ) ( 1 )
over all partitions P = (C 1 , ..., C K ) of Ω and all choices of set L = (y 1 , ...y K ) of cluster prototypes.

In this context, the dynamic clustering algorithm iteratively performs both a representation step and an allocation step: a) Representation step (the partition P is fixed).

Finding L that minimizes ∆(P, •) is equivalent to finding for k = 1, ..., K, the prototype y k that minimizes the adequacy criterion i∈C k d(x i , y k ). For con-tinuous quantitative data in p and the city-block distance, y k is the median vector of the cluster C k .

b) Allocation step (the set of prototypes L is fixed).

Finding P that minimizes ∆(•, L) is equivalent to finding for k = 1, ..., K, the cluster

C k = {i ∈ Ω | d(x i , y k ) ≤ d(x i , y m ) , ∀m = 1, ..., K}
Once these two steps properly defined, the partitioning criterion (1) decreases at each iteration and the algorithm converges to a stationary value of this criterion under the two following conditions:

i) Unicity of the choice for the cluster affectation of each object of Ω; ii) Unicity of the choice of the prototype y k that minimizes the adequacy cri-

terion i∈C k d(x i , y k )
The main idea of dynamic clustering with adaptive distances is to associate a distance d k to each cluster C k (and its prototype y k ) such that the sum of the distances d k (x i , y k ) between objects i ∈ C k and the prototype y k is as small as possible. The distances used in the dynamic algorithm are therefore not determined once and for all. Moreover, they are different from one cluster to another. The clustering criterion is:

∆(P, L) = K k=1 i∈C k d k (x i , y k ) ( 2 )
where P = (C 1 , ..., C K ) and now L = (G, d), where G = (y 1 , ..., y K ) and d = (d 1 , . . . , d K ).

In our context, the distance d k is a weighted sum of distances d j , where d j compares a pair of objects according to variable j:

d k (x i , x i ) = p j=1 d j (x j i , x j i ) = p j=1 λ j k d(x j i , x j i ) ( 3 ) with d j (x j i , x j i ) = λ j k d(x j i , x j i ), λ j k > 0 and p j=1 λ j k = 1.
According to the definition of d k given in (3), the set L is written L = (G, λ) where λ = (λ 1 , . . ., λ K ), with λ k = (λ 1 i , ..., λ p i ) being the vector of weights of the fixed distance d. The adaptivity of the distance d k is expressed by the vector of weights λ k .

When using adaptive distances, the representation step is divided in two stages: a1) Stage 1 (the partition P and λ are fixed).

Find for k = 1, ..., K, the prototype y k that minimizes the adequacy criterion

i∈C k d k (x i , y k ) = i∈C k p j=1 λ j k d(x j i , y j k ) = p j=1 λ j k i∈C k d(x j i , y j k ).
a2) Stage 2 (the partition P and the set of prototypes G are fixed).

Find for k = 1, ..., K, the vector of weights λ k that minimizes the adequacy cri-

terion i∈C k d k (x i , y k ) = p j=1 λ j k i∈C k d(x j i , y j k ) = p j=1 λ j k Φ j where Φ j = i∈ C k d(x j i , y j k ).
The allocation step of the algorithm is once again:

b) Allocation step (the set of prototypes G and and the vector λ are fixed):

Find for k = 1, ..., K, C k = {i ∈ Ω | d k (x i , y k ) ≤ d m (x i , y m ) , ∀m = 1, ..., K}.
Once these two steps have been properly defined, the partitioning criterion (2) decreases at each iteration and the algorithm converges to a stationary value of this criterion under the three following conditions:

i) Unicity of the choice for the cluster affectation of each object of Ω; ii) Unicity of the choice of the prototype y k that minimizes the adequacy cri-

terion i∈C k d k (x i , y k ) iii) Unicity of the choice of the vectors of weights λ k that minimizes the ade- quacy criterion p j=1 λ j k i∈C k d(x j i , y j k ).

Dynamic clustering method with an adaptive Hausdorff distances

In this paper we are concerned with objects that are represented by a vector of intervals (we consider a point as an interval with equal lower and upper bounds). Let Ω be a set of n objects indexed by i and described by p interval variables indexed by j. An interval variable X [START_REF] Bock | Analysis of Symbolic Data, Exploratory methods for extracting statistical information from complex data[END_REF]) is a correspondence defined from Ω in such that for each i ∈ Ω,

X(i) = [a, b] ∈ ,
where is the set of closed intervals defined from .

Each object i is represented as a vector of intervals

x i = (x 1 i , • • • , x p i )
, where

x j i = [a j i , b j i ] ∈ = {[a, b] : a, b ∈ , a ≤ b}.
In this paper, an interval data table {x j i } n×p is made up of n rows that represent n objects to be clustered and p columns that represent p interval variables. Each cell of this table contains an interval

x j i = [a j i , b j i ] ∈ .
A prototype y k of cluster C k ∈ P is also represented as a vector of intervals

y k = (y 1 k , • • • , y p k ), where y j k = [α j k , β j k ] ∈ .
It is now a matter of choosing an adaptive distance between vectors of intervals and properly defining the representation step of the dynamic algorithm with adaptive distances given in the previous section. In other words, we will give an explicit formula for the prototype y k and for the vector of weights λ k that minimizes both the adequacy criterion p j=1 λ j k i∈C k d(x j i , y j k ).

The adaptive Hausdorff distances for interval data

A number of proximity measures have been introduced in the literature for interval data (as well as for other types of symbolic data). [START_REF] Gowda | Symbolic clustering using a new dissimilarity measure[END_REF] and [START_REF] Gowda | Symbolic clustering using a new similarity measure[END_REF] We have seen that the distance d k associated with the cluster C k is defined as a weighted sum of distances d j , where d j compares a pair of objects according to variable j:

d j (x j i , x j i ) = λ j k d(x j i , x j i )
Here, the two feature values x j i and x j i are, respectively, the two intervals [a j i , b j i ] and [a j i , b j i ]. The distance d (see equation 3) chosen to compare two intervals is the Hausdorff distance. The Hausdorff distance [START_REF] Nadler | Visualization of Large Data Sets: The Zoom Star Solution[END_REF], [START_REF] Rote | Computing the minimum Hausdorff distance between two point sets on a line under translation[END_REF]) is often used in image processing [START_REF] Huttenlocher | Comparing images using the Hausdorff Distance[END_REF]) and is defined to compare two sets of objects A and B. This distance depends on the distance chosen to compare two objects u and v respectively in A and B. We consider the euclidean distance and the Hausdorff distance is defined by:

d H (A, B) = max(h(A, B), h(B, A))
where:

h(A, B) = sup u∈A inf v∈B ||u -v|| At times, h is called the directed Hausdorff distance.
In this work, A and B are two intervals

x j i = [a j i , b j i ] and x j i = [a j i , b j i ].
The previous Hausdorff distance is simplified to:

d H (x j i , x j i ) = max{|a j i -a j i |, |b j i -b j i |} (4)
Finally, the adaptive distance d k associated with the cluster C k and defined in equation ( 3) is:

d k (x i , x i ) = p j=1 λ j k d H (x j i , x j i ) = p j=1 λ j k max{|a j i -a j i |, |b j i -b j i |} (5)
with λ j k > 0 and p j=1 λ j k = 1.

Definition of the best prototypes

In section 2, we saw that the representation step of the dynamic clustering algorithm with adaptive distances was divided into two stages, corresponding to two minimization problems. The first problem, when the partition P and the vector λ are fixed, is to find for k = 1, ..., K the prototype y k that minimizes the adequacy criterion i∈C k d k (x i , y k ). With d k defined in (5) and with

x i = (x 1 i , • • • , x p i ) and y k = (y 1 k , • • • , y p k )
, the adequacy criterion is:

i∈ C k d k (x i , y k ) = i∈ C k p j=1 λ j k d H (x j i , y j k ) = p j=1 λ j k i∈ C k d H (x j i , y j k ) ( 6 )
The vector of weights being fixed, the problem is now to find for j = 1, . . . , p the interval

y j k = [α j k , β j k ] that minimizes: i∈C k d H (x j i , y j k ) = i∈C k max{|a j i -α j k |, |b j i -β j k |} (7)
According to [START_REF] Chavent | Dynamical Clustering Algorithm of Interval Data: Optimization of an Adequacy Criterion Based on Hausdorff Distance[END_REF], an explicit formula for the components y j k of the best prototype is found by transforming the previous minimization problem into two well-known L 1 norm problems. Let m j i = (a j i + b j i )/2 be the midpoint of the interval x j i = [a j i , b j i ] (for j = 1, . . . , p) and l j i = (b j ia j i )/2 be half of its length. From this we have:

a j i = m j i -l j i and b j i = m j i + l j i (8)
Also, let µ j k = (α j k + β j k )/2 be the midpoint of the interval y j k = [α j k , β j k ] (for j = 1, . . . , p) and ρ j k = (β j kα j k )/2 be half of its length. We have:

α j k = µ j k -ρ j k and β j k = µ j k + ρ j k (9)
From the equations ( 7), ( 8) and ( 9), the equation ( 6) can be written as:

i∈C k d H (x j i , y j k ) = i∈C k max{|(m j i -l j i ) -(µ j k -ρ j k )|, |(m j i + l j i ) -(µ j k + ρ j k )|} = i∈C k max{|(m j i -µ j k ) -(l j i -ρ j k )|, |(m j i -µ j k ) + (l j i -ρ j k )|}
According to the following property defined for x and y in ,

max(|x -y|, |x + y|) = |x| + |y| Then: i∈C k d H (x j i , y j k ) = i∈C k (|m j i -µ j k | + |l j i -ρ j k |) = i∈C k |m j i -ρ j k | + i∈C k |l j i -ρ j k | (10)
This yields two well-known minimization problems in L 1 norm: find µ j k ∈ and ρ j k ∈ that respectively minimize:

i∈C k |m j i -µ j k | and i∈C k |l j i -ρ j k |
The solution μj k and ρj k are, respectively, the median of the set {m j i , i ∈ C k } (the midpoints of the intervals

x j i = [a j i , b j i ], i ∈ C k )
, and the median of the set

{l j i , i ∈ C k } (the half-lengths of the intervals x j i = [a j i , b j i ], i ∈ C k ). Finally, the solution ŷj k = [α j k , βj k ] is given by αj k = μj k -ρj k and βj k = μj k + ρj k (11)

Definition of the best distances

The second stage of the representation step of the dynamic clustering algorithm with adaptive distances, when the partition P and the set of prototypes G are fixed, is to find for k = 1, ..., K the vector of weights λ k that minimizes the adequacy criterion defined in (6) by:

i∈ C k d k (x i , y k ) = p j=1 λ j k Φ j where Φ j = i∈ C k d H (x j i , y j k )
Following [START_REF] Diday | Classification avec Distances Adaptatives[END_REF] and [START_REF] Govaert | Classification automatique et distances adaptatives[END_REF], the weights λ j k are calculated by the Lagrange multiplier method:

∂ ∂λ j k ( p j=1 λ j k Φ j -µ p h=1 λ h k ) = 0 for j = 1, . . . , p (12)
From equation ( 12), we have the following result:

Φ j -µ p h=1 λ h k λ j k = 0 ⇒ λ j k = µ Φ j ( p h=1 λ h k ) (13)
Remembering that Π p h=1 λ h k = 1, the parameter λ j k in equation ( 13) is given by

λ j k = µ Φ j
The restriction Π p h=1 λ h k = 1 can be written as:

1 = p h=1 µ Φ h = µ p p h=1 Φ h and then µ = ( p h=1 Φ h ) 1 p
Finally, the solution λj k to the parameter λ j k is:

λj k = µ Φ j = p h=1 ( i∈C h max{|a h i -αh k |, |b h i -βh k |}) 1 p i∈C k max{|a j i -αj k |, |b j i -βj k |} (14)

The algorithm

The algorithm schema of dynamic clustering algorithm with Hausdorff adaptive distances for interval data is as follows:

λj k max{|a j i -αj k |, |b j i -βj k |} if i ∈ C k and k * = k test ← 1 C k * ← C k * ∪ {i} C k ← C k \ {s} (4

) Stopping criterion

If test = 0 then STOP, otherwise go to (2).

For classical dynamical clustering methods, the initialization step and stopping rules can be modified. For example, the points chosen randomly at initialization can be chosen in such a way that they are as dissimilar as possible.

Concerning the stopping rule, a minimum value for the clustering criterion or a maximum number of iterations can also be given.

Another remark is when all the weights λ j k are fixed to 1, the distances are nonadaptive and the previous algorithm is equivalent to the dynamic clustering algorithm of interval data proposed in [START_REF] Chavent | Dynamical Clustering Algorithm of Interval Data: Optimization of an Adequacy Criterion Based on Hausdorff Distance[END_REF]. This remark will be used in the next section for evaluating the adaptive dynamic clustering algorithm.

Experimental results

To show the usefulness of this method, two synthetic interval data sets with linearly non-separable clusters of different shapes and sizes have been drawn.

Real applications are then considered.

Our aim is to achieve a comparison of the dynamic clustering algorithm considering different distances between vectors of intervals: adaptive Hausdorff distance proposed in this paper, non-adaptive Hausdorff distance [START_REF] Chavent | Dynamical Clustering Algorithm of Interval Data: Optimization of an Adequacy Criterion Based on Hausdorff Distance[END_REF]), one component adaptive city-block distance [START_REF] Souza | Clustering of interval data based on city-block distances[END_REF]) and non-adaptive city-block distance.

To compare the results furnished by the dynamic clustering algorithm with these different distances, an external validity index is used. For synthetic interval data sets, rectangles are built from three clusters of points drawn from three bi-variate normal distributions. Next, the a priori partition of the objects is known. For the car interval data set describing car models, it is defined a a priori partition into four groups according to a car category. For the interval data set describing species of freshwater fish, it is considered a a priori partition of the species into four groups according to diet.

The idea of external validity is simply to compare the a priori partition with the partition obtained from the clustering algorithm. In this paper, we use the corrected Rand (CR) index defined in [START_REF] Hubert | Comparing Partitions[END_REF] for comparing two partitions, the definition of which is as follows.

Let U = {u 1 , . . . , u i , . . . , u R } and V = {v 1 , . . . , v j , . . . , v C } be two partitions of the same data set having respectively R and C clusters. The corrected Rand index is:

CR = R i=1 C j=1 n ij 2 -n 2 -1 R i=1 n i. 2 C j=1 n .j 2 1 2 [ R i=1 n i. 2 + C j=1 n .j 2 ] -n 2 -1 R i=1 n i. 2 C j=1 n .j 2 (15)
where n 2 = n(n-1) 2 and n ij represents the number of objects that are in clusters u i and v i ; n i. indicates the number of objects in cluster u i ; n .j indicates the number of objects in cluster v j ; and n is the total number of objects in the data set. CR takes its values from the interval [-1,1], where the value 1 indicates perfect agreement between partitions, whereas values near 0 (or negatives) correspond to cluster agreement found by chance.

Synthetic interval data sets

In this paper, we consider the same data point configurations presented in [START_REF] Souza | Clustering of interval data based on city-block distances[END_REF]. Two data sets of 350 points in 2 were constructed. In each data set, the 350 points are drawn from three bi-variate normal distributions of independent components. There are three clusters of unequal sizes and shapes: two clusters with an ellipsoidal shape and size 150 and one cluster with a spherical shape and size 50. The mean vector and the covariance matrix of the bi-variate normal distributions are noted:

µ =    µ 1 µ 2    and Σ 11 =    σ 2 1 0 0 σ 2 2   
Data set 1 shows well-separated clusters (Figure 1). The data points of each cluster in this data set were drawn according to the following parameters: In order to build interval data sets from data sets 1 and 2, each point (z 1 , z 2 ) of these data sets is considered as the 'seed' of a rectangle. Each rectangle is therefore a vector of two intervals defined by:

([z 1 -γ 1 /2, z 1 + γ 1 /2], [z 2 -γ 2 /2, z 2 + γ 2 /2]) ( 16
)
The parameters γ 1 and γ 2 are the width and the height of the rectangle. They are drawn randomly within a given range of values. For example, the width and the height of all the rectangles can be drawn randomly within the interval [1,8]. Figure 2 shows two synthetic interval data sets built from data set 1 and data set 2 when γ 1 and γ 2 are drawn randomly from [1,8]. In the framework of a Monte Carlo experiment, 100 replications of the previous process were carried out for parameters γ 1 and γ 2 , drawn randomly 100 times from each of the following intervals: [1,8], [1,16], [1,24], [1,32], [1,40]. This process has also been repeated for seeds taken from data set 1 and data set 2.

Dynamic clustering algorithms considering different distances between vectors of intervals have been performed on these data sets. The 3-cluster partitions obtained with these clustering methods were compared with the 3-cluster partition known a priori. The comparison index used is the corrected Rand index CR given in equation ( 15). For each 100 replications, the average corrected Rand index CR is calculated.

Table 2 gives the values of the average CR index obtained with adaptive and non-adaptive distances for interval data sets 1 and 2 as well as γ 1 and γ 2 drawn from [1,8], [1,16], [1,24], [1,32], [1,40]. As expected, in each case the average CR indices are better with adaptive distances.

Concerning the data configurations presenting well separated classes, the Hausdorff (non-adaptive) distance shows better CR indices than city-block (nonadaptive) distance regardless the range of the predefined intervals in Table 2. Moreover, the Hausdorff distance is also the best option for data configuration presenting overlapping classes as long as the widest intervals are considered.

For both type of data configurations (well separated classes and overlapping classes), the average CR indices provided by the adaptive Hausdorff distance are again better than those provided by the adaptive city-block distance for those data configurations where the range of the predefined intervals are the widest. Table 3 gives the corresponding values of the standard deviation for the average CR index.

The evaluation of the performance of the dynamic clustering methods for these different distances between vectors of intervals is achieved by an independent Student's t-test with a 5% level of significance. Tables 4 and5 shows the suitable (null and alternative) hypothesis and the observed values of the test statistics following a Student's t distribution with 198 degrees of free- dom for interval data sets 1 and 2, respectively. In this table, µ 1 , µ 2 , µ 3 and µ 4 are, respectively, the average of the CR index for the dynamic clustering algorithm considering the Hausdorff distance, city-block distance, adaptive Hausdorff distance and city-block (one componente) adaptive distance, respectvely. These tables show that, in 75% of the data simulation configurations considered in this work, the dynamic clustering algorithm based on Hausdorff distance outperforms the version of this algorithm which uses the city-block distance considering both adaptive and non-adaptive cases.

H 0 : µ 1 ≤ µ 2 H 0 : µ 3 ≤ µ 4 H 1 : µ 1 > µ 2 H 1 : µ 3 > µ 4 [1,
H 0 : µ 1 ≤ µ 2 H 0 : µ 3 ≤ µ 4 H 1 : µ 1 > µ 2 H 1 : µ 3 > µ 4 [1, 8] -2.

Car data set

The car interval data set consists of a set of 33 car models described by 8 interval, 2 categorical multi-valued and one nominal variables (see Table 6).

In this application, the 8 interval variables -Price, Engine Capacity, Top Speed, Acceleration, Step, Length, Width and Height -have been considered for clustering purposes, the nominal variable Car Category has been used as a a priori classification.

Ecotoxicology data set

Several studies realized in French Guyana indicated abnormal levels of mercury contamination in some Amerindian populations. This contamination was connected to their high consumption of contaminated freshwater fish (Boudou and Ribeyre (1998)). In order to get a better knowledge of this phenomenon, a data set has been collected by researchers from the LEESA (Laboratoire d'Ecophysiologie et d'Ecotoxicologie des Systèmes Aquatiques) laboratory.

This data set concerns 12 species of freshwater fish, each species being described by 13 interval variables. These species are grouped into four a priori clusters of unequal sizes according to diet: two clusters (Carnivorous and Detritivorous) of size 4 and two clusters of size 2 (Omnivorous and Herbivorous).

Table 6 shows part of the freshwater fish data set.

Dynamic clustering algorithms considering different distances between vectors of intervals have also been performed on this data set. The 4-cluster partitions obtained with these clustering methods were compared with the 4-cluster partition known a priori. Again, the comparison index used is the corrected Rand index CR given in equation ( 15). The a priori classification, indicated by the suffix attached to the freshwater specie denomination, is as follows: performance). However, for this data sets, the version of the dynamic clustering algorithm with adaptive distances outperforms the version of this algorithm with non-adaptive distances. For the case of the dynamic clustering algorithm considering the Hausdorff adaptive distances performed on the Ecotoxicology interval data set, Tables 10 and11, respectively, give the prototype descriptions and the corresponding weight vectors of the (Hausdorff) adaptive distances associated to each class, according to the 13 interval variables. presented, especially a type of graphic called Zoom Star. In this graphical representation, each axis corresponds to an interval variable. In each axis, the lower and upper bounds of the interval value assumed by an interval variable for a given object are represented. The lower bounds (as well the upper bounds) of the intervals assumed by each interval variable are linked to form a polygon. The Zoom Star shows the area between the upper-bound and lower-bound polygons. Figure 5 gives the visualization of the prototype of each cluster from Table 10 according to the Zoom Star method.

All the interval variables for the prototype of Cluster 1 show intervals with low spread, whereas they show intervals with medium spread for the prototype of Cluster 4. For the prototype of Cluster 2, most of the interval variables that do not represent ratios show a high spread, whereas the interval variables that express ration show a low spread. Concerning the prototype of Cluster 3, the role of the ratio and non-ratio interval variables are inverted in comparison to their role in the prototype of Cluster 2.

In conclusion, for this data set, the performance of the adaptive methods measured by the CR index is superior to the non-adaptive methods. The dynamic clustering algorithm with adaptive Hausdorff distances starts from an initial partition and alternates a representation step and an allocation step until convergence when the adequacy criterion reaches a stationary value representing a local minimum. The representation step has two stages.

In the first stage, the partition and the Hausdorff distances are fixed and the algorithm looks for the best prototype of each class which minimizes the adequacy criterion. The solution for the best prototype of each class, presented in this paper, is a vector of intervals whose lower bounds, for a given variable, are the difference between the median of midpoints of the intervals computed for the objects belonging to this class and the median of their half-lengths, and whose upper bounds, for a given variable, are the sum of the median of midpoints of the intervals computed for the objects belonging to this class plus the median of their half-lengths. In the second stage, the partition and the prototype of each class are fixed and the algorithm looks for the best Hausdorff distance associated to each class which minimizes the adequacy criterion. The Hausdorff distance associated to each class is parameterized by a vector of weights and the best solution for this vector of weights provided by the clustering method is also presented in this paper. In the allocation step, the individuals are assigned to the classes according to their (minimum) adaptive Hausdorff distance to the prototypes.

Experiments with real and synthetic interval data sets showed the usefulness of this clustering method. The accuracy of the results furnished by the dynamic clustering algorithm based on adaptive Hausdorff distance is assessed by the CR index and compared with the results provided by this algorithm considering non-adaptive Hausdoff distance and adaptive and non-adaptive city-block distances.

Concerning the synthetic interval data sets, the CR index is calculated in the framework of a Monte Carlo experiment with 100 replications. For the data configurations showing well separated classes, the Hausdorff distance outperforms the city-block distance for the non-adaptive version of the dynamic clustering algorithm. Moreover, for the non-adaptive version of the dynamic clustering algorithm, the Hausdorff distance is also the best option for data configurations presenting overlapping classes as long as the widest intervals are considered. For both types of data configurations (well separated classes and overlapping classes) the adaptive Hausdorff distance outperforms the adaptive city-block distance also as long as the widest intervals are considered.

Concerning the car interval data set, the Hausdorff distance outperforms the city-block distance for the non-adaptive version of the dynamic clustering algorithm. Moreover, these distances presented the same performance when the adaptive version of the dynamic clustering algorithm is applied on this data set. Concerning the ecotoxicology data set, the Hausdorff and city-block distances presented the same performance for the non-adaptive and adaptive version of the dynamic clustering algorithm. However, for both data sets, the version of the dynamic clustering algorithm with adaptive distances outperforms the version of this algorithm with non-adaptive distances.

  introduced, respectively, dissimilarity and similarity functions with components based on position, span and content. The component based on position indicates the relative positions of two feature values on a real line. The component based on span indicates the relative sizes of the feature values without referring their common parts. The component based on content is a measure of the common parts between two features values.Ichino and Yaguchi (1994) presented the generalized Minkowski metrics for mixed feature variables. Similarity and dissimilarity measures between symbolic data restricted by dependency rules between feature values can be found in De[START_REF] De Carvalho | Proximity coefficients between Boolean symbolic objects[END_REF][START_REF] De Carvalho | Statistical proximity functions of boolean symbolic objects based on histograms[END_REF] and De Carvalho and[START_REF] De Carvalho | New metrics for Constrained Boolean Symbolic Objects[END_REF].

Fig. 1 .

 1 Fig. 1. Seed data sets 1 and 2 showing, respectively, well-separated and overlapping classes

Fig. 2 .

 2 Fig. 2. Interval data sets 1 and 2, showing, respectively, well-separated and overlapping classes

Fig. 3 .

 3 Fig. 3. Ecotoxicology data set: description of prototypes for each cluster according to Zoom Star method 5 Concluding remarks

  InitializationChoose a partition {C 1 . . . , C K } of Ω randomly or choose K distinct objects y 1 , . . . , y K belonging to Ω and assign each object i to the closest object y k * (k * = arg min k=1,...,K ∈ C k } and ρj k is the median of the set {l j i , i ∈ C k } b) For j = 1, . . . , p and k = 1, . . . , K, compute λj

	Schema of adaptive dynamic clustering algorithm
	(1) p j=1 max{|a j i -α j k |, |b j i -β j k |} to con-
	struct the initial partition {C 1 , . . . , C K }.
	(2) Representation step
	a) For i = 1 to K compute the prototype ŷk = ([α 1 k , β1 k ], . . . , αp k , βp k ]) with αj k = μj k -ρj k and βj k = μj k + ρj k where μj k is the median of the set {m j i , i k with equation (14)
	(3) Allocation step
	test ← 0
	for i = 1 to n do
	define the winning cluster C k * such that
	k * = arg min k=1,...,K	p j=1

Table 2

 2 Comparison of the methods according to the average corrected Rand index

			Interval Data Set 1	Interval Data Set 2
	Predefined	Non-Adaptive	Adaptive	Non-Adaptive	Adaptive
	Intervals	Distances	Distances	Distances	Distances
		L 1	Hausd.	L 1	Hausd.	L 1	Hausd.	L 1	Hausd.
	[1, 8]	0.684 0.691 0.935 0.923 0.379 0.378 0.470 0.448
	[1, 16]	0.664 0.706 0.931 0.931 0.375 0.374 0.432 0.434
	[1, 24]	0.636 0.700 0.892 0.909 0.369 0.377 0.406 0.418
	[1, 32]	0.622 0.685 0.773 0.912 0.361 0.385 0.389 0.412
	[1, 40]	0.618 0.702 0.701 0.886 0.348 0.378 0.373 0.393
	Table 3
	Comparison of the methods according to the standard deviation of the corrected
	Rand index
			Interval Data Set 1		Interval Data Set 2
	Predefined	Non-Adaptive	Adaptive	Non-Adaptive		Adaptive
	Intervals	Distances	Distances	Distances		Distances
		L 1	Hausd.	L 1	Hausd.	L 1	Hausd.		L 1	Hausd.
	[1, 8]	0.0127 0.0150 0.0005 0.0010 0.0013 0.0013 0.0050 0.0019
	[1, 16]	0.0118 0.0184 0.0009 0.0007 0.0010 0.0012 0.0014 0.0023
	[1, 24]	0.0041 0.0154 0.0058 0.0011 0.0014 0.0010 0.0015 0.0019
	[1, 32]	0.0019 0.0131 0.0114 0.0011 0.0013 0.0014 0.0013 0.0017
	[1, 40]	0.0014 0.0133 0.0073 0.0032 0.0010 0.0016 0.0024 0.0024

Table 6 '

 6 Car' data set with 8 interval and one nominal variables Each clustering method is run (until the convergence to a stationary value of the adequacy criterion) 100 times and the best result, according to the adequacy criterion, is selected. The corrected Rand index CR is calculated for the best result. Table7shows the clusters (individual labels) given by the nonadaptive (L 1 and Hausdorff) and adaptive (one component L 1 and Hausdorff) methods. The CR indices obtained from the results displayed in Table7are 0.35 and 0.38 for the non-adaptive L 1 and Hausdorff methods, respectively, and 0.56 for the adaptive (one component L 1 and Hausdorff) methods. Notice that, for this data set, the dynamic clustering algorithm with non-adaptive Hausdorff distance outperforms this same algorithm with non-adaptive L 1 distance. Moreover, for the case of adaptive distances, it is furnished the same partition by the dynamic clustering algorithm (the adaptive Hausdorff and one component city-block distances presented the same performance). However, for this data set, the version of the dynamic clustering algorithm with adaptive distances outperforms the version of this algorithm with non-adaptive distances.

	Price	Engine	. . .	Height	Category
		Capacity

Table 7

 7 Clustering Results for the Car data set

	Method	Cluster 1	Cluster 2	Cluster 3	Cluster 4
	L 1 (non-adaptive)	2/B 3/L 5/U	7/L 9/L 10/L	1/U 12/U 13/U	4/S 11/S 16/S
		6/B 8/B 18/L	15/S 19/S 20/S	14/B 17/U 24/U	22/L 23/L
		21/B 30/B 33/L	27/S	25/U 26/B 28/U
				29/U 31/U 32/B
	Hausdorf (non-adaptive)	1/U 12/B 13/U	7/L 9/L 10/L	2/B 3/L 5/U	4/S 11/S 16/S
		14/B 17/U 24/U	15/S 19/S 20/S	6/B 8/B 18/L
		25/U 26/B 28/U	22/L 23/L 27/S	21/B 30/B 33/S
		29/U 31/U 32/B
	L 1 (adaptive)	12/U 13/U 17/U	1/U 2/B 3/L	6/B 7/L 9/L	4/S 11/S 15/S
		24/U 25/U 28/U	5/U 8/B 14/B	10/L 22/L 23/L	16/S 19/S 20/S
		29/U 31/U	18/L 21/B 26/B		27/S
			30/B 32/B 33/L
	Hausdorf (adaptive)	1/U 2/B 3/L	12/U 13/U 17/U	4/S 11/S 15/S	6/B 7/L 9/L
		5/U 8/B 14/B	24/U 25/U 28/U	16/S 19/S 20/S	10/L 22/L 23/L
		18/L 21/B 26/B	29/U 31/U	27/S
		30/B 32/B 33/L

Table 9

 9 Clustering Results for the Ecotoxicology data set

	Method	Cluster 1	Cluster 2	Cluster 3	Cluster 4
	L 1 (non-adaptive)	1/C 4/C 7/D	2/C	3/C	5/O 6/O 9/D
		8/D 10/D			11/H 12/H
	Hausdorf (non-adaptive)	1/C 4/C 7/D	2/C	3/C	5/O 6/O 9/D
		8/D 10/D			11/H 12/H
	L 1 (adaptive)	5/O 6/O	9/D 11/H 12/H	1/C 2/C 3/C	4/C 7/D 8/D
					10/D
	Hausdorf (adaptive)	5/O 6/O	9/D 11/H 12/H	1/C 2/C 3/C	4/C 7/D 8/D
					10/D

Table 11

 11 Ecotoxicology data set: vectors of weights λ k of adaptive distance d k (k = 1, . . . , 4) according to the 13 interval variables

	Interval		Vectors of weights
	Variables	λ 1	λ 2	λ 3	λ 4
	Length	3.911076	3.784355	9.578303	5.855540
	Weight	0.099555	0.022819	0.362272	0.022301
	Muscle	0.016830	0.031758	0.375701	0.077517
	Intestine	0.062889	0.146273	0.200882	0.145701
	Stomach	0.042986	0.012046	0.319487	0.108850
	Gills	0.158115	0.097757	0.295046	0.468606
	Liver	0.032341	0.003193	0.032691	0.006864
	Kidneys	0.008428	0.004982	0.082949	0.007686
	Liver/Muscle	41.066313	33.423074	3.175321	2.804027
	Kidneys/Muscle 28.321594	39.767796	2.033658	8.990483
	Gills/Muscle	410.663114 1466.437071 12.970620 374.103973
	Intestine/Muscle 68.443852	488.812463 18.866356 184.489640
	Stomach/Muscle 40.064695	139.660709	9.453641 144.814440

. . . . . . . . . Each clustering method is run (until the convergence to a stationary value of the adequacy criterion) 100 times and the best result, according to the adequacy criterion, is selected. The corrected Rand index CR is calculated for the best result. Table 9 shows the clusters (individual labels) given by the non-adaptive (L 1 and Hausdorff) and adaptive (one component L 1 and Hausdorff) methods. The CR indices obtained from the results displayed in Table 9 are 0.488 and 0.138 for the adaptive and non-adaptive distances, respectively. Notice that, for this data set, regardless the adaptive (or the non-adaptive) distances used, the dynamic clustering algorithm furnishes the same partition (the Hausdorff and city-block distances presented the same In Noirhomme-Fraiture (2002), visualization techniques for interval data are


Alan Silva

Etienne Santos

Philibert Girard

Sylvette Leflaive

Brunet

Santos

Etienne Girard

Philibert Leflaive

Sylvette Brunet
email: sylvette.brunet@univ-poitiers.fr

Competitive adsorptions between thiophenic compounds over CoMoS/Al 2 O 3 under deep HDS of FCC gasoline

Keywords: Hydrodesulfurization, FCC gasoline, 2-methythiophene, 3-methylthiophene, benzothiophene, kinetic modeling

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Due to increasing environmental concerns, constraints relative to exhaust gas compositions have been reinforced worldwide. China V and US Tier 3 regulations have thus imposed a maximum of 10 ppmSulfur in commercial gasoline since 2017 [1]. Commercial gasoline is blended from gasoline streams mainly produced by isomerization, reforming and fluid catalytic cracking (FCC) units in refineries. In particular, catalytic cracked gasoline represents 40-80% of the total gasoline composition but accounts for almost all of the total sulfur amount present in non-hydrodesulfurized commercial gasoline [2,[START_REF] Toulhoat | Catalysis by Transition Metal Sulphides from Molecular Theory to Industrial Application[END_REF]. Consequently, there is a strong need for achieving ultra-deep HDS of FCC gasoline, which comes with developing highly active HDS catalysts. Meanwhile, these catalysts should also exhibit high selectivity to minimize olefins hydrogenation (HDO) and preserve high octane values.

The FCC gasoline is a complex blend of more than 400 compounds, mainly including paraffins, olefins, naphtenes and aromatics. The sulfur compounds in FCC gasoline are represented in various families: mercaptans, thiophenes, benzothiophenes and their alkyl derivatives. The repartition and the total amount of these sulfur compounds depend on the origin of the feedstock and the operating conditions of the FCC unit. As reported in literature [START_REF] Siddiqui | [END_REF]5], for a typical feedstock, C 1 -thiophenes -i.e. methylthiophenes-, C 2 -thiophenes and benzothiophenes correspond to more than 50% of total sulfur compounds present in the mixture. Therefore, the sulfur compounds selected in model feeds as single representative of sulfur contents in a FCC gasoline are thiophene [6], alkylthiophenes [2,7], benzothiophene and its alkylated derivatives [8,9]. However, a major drawback of such approach lies in a lack of representativeness of the reactivity of a real catalytically cracked gasoline where relative reactivities of reactants and potential competitive adsorptions between reactants and products dictate the overall feed reactivity.

The reactivity of these different sulfur molecules have been ambiguously reported in the literature with different rankings. This apparent contradiction may originate from different operating conditions. At their operating conditions -100ppmS in toluene, 1.3 MPa, 150°C, LHSV =3.5 h -1 -, Hatanaka et. al. [10,11] established the following reactivity scale: benzothiophene > thiophene > 3-methylthiophene > 2-methylthiophene > 2-ethylthiophene > 2,5-dimethylthiophene. The reaction rate of benzothiophene was approximatively twice that of thiophene. These results also confirmed the relative reactivities of 3-methylthiophene and 2methylthiophene over a CoMoS/Al 2 O 3 catalyst, as reported by Desikan et. al. [12] at atmospheric pressure and temperatures of 305 and 414ºC. However, in a review on hydrodesulfurization processes, Song et al [13] suggested a reverse reactivity order with thiophene being more reactive than benzothiophene. This assumption was based on results originally published by Nag et al. [14] on the hydrodesulfurization of benzothiophene and thiophene at 7,1 MPa and 300°C. In this study, the reactivity of thiophene was higher by a factor of 1.7 as compared to benzothiophene.

In addition, the potential competitive adsorptions over the catalyst surface between reactants and products in complex feeds such as a real catalytically cracked gasoline have never been thoroughly investigated. Indeed, most relevant articles focused on the inhibiting effect of H 2 S on the transformation of single sulfur compounds in model feeds. Thus, Fontaine et. al. [15] evidenced that the presence of H 2 S had a weak impact on the transformation of 2-methylthiophene over an unsupported CoMoS catalyst at 2 MPa and 250°C. A similar effect was observed by Dos Santos et al. [7] during the hydrodesulfurization of 3-methyl-thiophene over a CoMoS catalyst supported on alumina at 2 MPa and 185°C. A decrease in the kinetic apparent constants as a function of H 2 S partial pressure was attributed to a competitive adsorption of 3-methylthiophene and H 2 S on the active sites of the catalyst.

To develop highly efficient catalysts and bridge the gap between model feed and real catalytically cracked gasoline reactivities, there is consequently a strong need for a proper understanding of the reactivity scale of sulfur compounds under selective HDS operating conditions as well as their interplay in mixture. This paper deals with the comparison of the reactivity of sulfur model molecules representative of various families in FCC gasoline alone and in mixture in order to determine the most refractory sulfur compounds and to quantify their competitive adsorption. The transformations of 2-methylthiophene (2MT), 3-methylthiophene (3MT) and benzothiophene (BT), -the selected model molecules-over a supported CoMoS catalyst under selective HDS conditions (250°C and 2Mpa of total pressure) were thus studied alone, and then in mixture.

Following a combined experimental and theoretical approach, a kinetic modeling based on a Langmuir-Hinshelwood formalism was then developed to explain these experimental results and provide quantitative elements of reactivity of sulfur compounds in selective HDS conditions in place of the more qualitative information that is found in the available literature.

Experimental part

Catalyst and chemicals

The CoMoS/Al 2 O 3 catalyst was provided by IFPEN with the following composition: CoO (3%wt), MoO 3 (10%wt). The various characteristics are reported Table 1. The catalyst was crushed and sieved to size range between 250 and 315 m and then sulfided in situ under H 2 S/H 2 flow (10 mol% H 2 S) for 10 h at 400°C at atmospheric pressure. 2-methylthiophene (98% purity), 3-methylthiophene (98% purity), benzothiophene (95% purity) and n-heptane (>99% purity) were purchased from Sigma-Aldrich. They were used without further purification. Hydrogen sulfide (1 vol% in mixture with H 2 ) was purchased from Air Liquide.

Reaction conditions

Catalytic activity measurements were carried out in a fixed bed reactor at 250°C under a total pressure of 2 MPa with a ratio H 2 /feed of 360 NL/L [16][17][18]. The sulfur model feed (0.3wt %, 1000 ppmS), diluted in n-heptane was injected in the reactor by a HPLC Gilson pump (307 series, pump's head volume: 5 cm 3 ). The mass of catalyst used was 50 mg of CoMo/Al 2 O 3 and the contact time varied from 0 to 30 s.

Different types of feeds were used to study the transformation of the sulfur model molecules alone or in mixture:

1/ Single component feed including a sulfur model molecule alone (corresponding to 1000 ppmS) (0.3 wt% of 2MT or 3MT or 0.42wt% of BT) in n-heptane; 2/ Feed consisting of a mixture of two sulfur components including 0.3wt% of 2MT (1000 ppmS) and various amounts of BT (corresponding to 500, 1000 and 1500 ppmS) in n-heptane; 3/ Feed consisting of a mixture of either 2MT or BT (1000 ppmS) and H 2 S (from 0.2 to 1.45 kPa) in n-heptane. These H 2 S partial pressures correspond to the amount of H 2 S produced by the transformation of BT (or 2MT) in mixture with 2MT (or BT). These latter experiments were carried out in order to identify the real inhibitor of the transformation of a given sulfur compound (the other sulfur molecule present or the H 2 S produced by its transformation);

The different partial pressures of the reactants, H 2 S, H 2 introduced are reported Table 3 for the sulfur molecules alone and in mixture, respectively. n-heptane was not converted under these experimental conditions. No catalyst deactivation was observed for all the experiments.

Products analysis

The reaction products were injected on-line by means of an automatic sampling valve into a Varian gas chromatograph equipped with a PONA capillary column and a flame ionization detector as in previous works [16][17][18]. Desulfurized products, resulting from the transformation of 2-methylthiophene, 3-methylthiophene and benzothiophene are designated as HDS products. The contact time is defined as the ratio between the total amount of feed and the mass of catalyst in the oxide form2MT, 3MT and BT reactivities, are defined as the number of moles of HDS products formed per hour and per gram of catalyst and were calculated at conversion lower than 30% in a differential regime.

Kinetic modeling

The obtained experimental results were confirmed by a theoretical approach using the Langmuir-Hinshelwood model and the Arrhenius equation. Kinetic and adsorption parameters were first set from single component experiments using ReactOp software ® and were then used to fit binary experimental results without further adjustment.

Results

Transformation of single sulfur compounds

In order to draw an unambiguous reactivity ranking of 2MT, 3MT and BT, the hydrodesulfurization of single sulfur compounds was first studied. The conversion of 2MT, 3MT and BT measured separately are reported Figure 1. The distribution of the products are reported in Figure 2 for 2MT, in Figure 3 for 3MT and in Figure 4 for BT. For the three sulfur compounds, a linear dependence of conversion with contact time was observed for conversions lower than 60%, thus corresponding to a differential regime. At equivalent contact time, the conversion of BT was systematically higher than those of 2MT and 3MT, as reported in Table 4. In other words, the reactivity order observed under these operating conditions is in agreement with the results previously obtained by Hatanaka et al. [10,11] who reported that 2-methylthiophene was less reactive than 3MT and BT at a total pressure of 1.3 MPa and temperatures ranging between 130 and 170ºC.

Regarding the product distribution of the transformation of 2-MT, the main products observed are pent-1-ene, cis-pent-2-ene and trans-pent-2-ene as primary products and npentane as a secondary product (Figure 2). In agreement with the reaction scheme proposed previously (Scheme 1) [2], there are two main potential pathways: direct desulfurization (DDS) -leading to the formation of cis and trans-pent-2-ene-and the hydrogenation (HYD) pathway -leading to the formation of pent-1-ene-which both end up into the hydrogenation product (n-pentane). However, the intermediates proposed in the reaction scheme for 2MT

were not observed under selective HDS operating conditions. Similar to 2MT, the transformation of 3MT also leads to the formation of alkenes (i.e.

2-methylbut-1-ene, 2-methylbut-2-ene and 3-methylbut-1-ene) as primary products and to the formation of isopentane as secondary product resulting from the hydrogenation of alkenes (Figure 3 according to the proposed reaction scheme in the literature (Scheme 2)) [7] The transformation of BT leads to dihydrobenzothiophene (DHBT) -resulting from the hydrogenation of the thiophenic ring-which appears as a primary product. The formation of ethylbenzene (EB) is also observed as an apparent secondary product although it can be also formed directly from BT (Figure 4). Consequently, the reaction scheme established from these results can involve two pathways for the formation of EB: a direct DDS pathway from BT involving a C-S bond rupture and the hydrogenation pathway involving the formation of DHBT and a subsequent C-S bond rupture (Scheme 3).

In order to evaluate a mutual inhibiting effect between the different sulfur compounds, the transformation of a feed containing 1000 ppmS from 2MT and different amounts of BT (from 500 to 1000 ppmS) was studied. The contact time was varied in order to obtain a large range of conversions (from 15 to 80%) of 2MT and BT. A mutual inhibiting effect between 2MT and BT was noticed. Moreover, the inhibition on the transformation of 2MT resulting from BT was higher (Figure 5a) than the opposite (Figure 6a). At a contact time of 2s, the decrease in conversion for the transformation of 2MT was equal to 54% while only 15% for the transformation of BT. This inhibiting effect is due to the unconverted other sulfur compound (2MT or BT) and not to the H 2 S produced by its transformation. Indeed, if 2MT or BT were substituted by the amount of H 2 S produced by the conversion of both sulfur compounds, the inhibiting effect observed was much lower in the transformation of 2MT (i.e. lower than 5%) (Table 5) and zero in the transformation of BT (Table 6). Consequently, the decrease in the conversion of 2MT (around 60%) observed in mixture with BT corresponds to an inhibition effect of the unconverted BT. Indeed, when BT is present at high concentration (1500ppm S) in the feed, its conversion is low (i.e. 47%), which leads to an inhibition of the transformation of 2MT around 66.5% (i.e. a conversion of 10.2% of 2MT in the feed as compared to a conversion of 28.2% in the absence of BT). In the opposite, for the conversion of BT, the presence of 2MT leads to an inhibition of only around 15% whereas the concentration of unconverted 2MT remains high due to its low conversion (maximum 17%) (Table 5).

Kinetic modeling

In order to understand and quantify the strong mutual inhibiting effects observed experimentally, a kinetic model was developed to determine the kinetic and adsorption parameters of 2MT, 3MT and BT transformations. As reported previously [17,18] the following assumptions were made: a Langmuir-Hinshelwood model involving reactions between adsorbed molecules, competitive adsorptions between sulfur compounds, alkenes, alkanes and H 2 S, a unique type of site and an heterolytic dissociation of H 2 and H 2 S. The corresponding elementary steps were reported previously under similar operating conditions [17,18].

For the transformation of 2MT and 3MT, different reactions were taken into account, respectively named R 1 or R 3 for the formation of olefin intermediates and R 2 and R 4 for the formation of alkanes from olefins.

R 1 R 2 R 3 R 4
The rate equations are as follows: Regarding the transformation of BT, the proposed reaction scheme includes two pathways for its transformation: HYD of BT into dihydrobenzothiophene (DHBT) and then HDS of DHBT into ethylbenzene (EB). The second pathway corresponds to the formation of ethylbenzene directly from BT. Therefore, the consecutive reactions proposed to model the transformation of BT are given below:

2 ) 3 ( 2 ) 3 ( 2 ) 3 ( 2 ) 3 ( 2 ) 3 ( 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 .                      
P K P K e k r 2 ) 4 ( 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 .                       H S
R 5 R 6 R 7
The corresponding equation rate are as follows: To explain the order of reactivity observed experimentally as well as the mutual inhibition of sulfur compounds, a kinetic modeling of the sulfur compounds alone initially and then in a mixture was carried out. The kinetic parameters (for the hydrodesulfurization and hydrogenation reactions) and the adsorption constants calculated from single component experiments are reported in Table 7. A good fit of the transformation of 2MT, 3MT and BT is obtained. The modeling (dashed curve) and the experimental results (points) for the conversion of 2MT, BT and 3MT alone and their product distribution are respectively reported on Figure 5, Figure 6 and Figure 7.

                      H S
It can be seen that the rate constants of the desulfurization step (olefin formation) (k = 2.8 10 12 h -1 and 5.5 10 12 h -1 for 2MT and 3MT respectively) are lower than those of the corresponding saturated hydrocarbon formation (k = 6.9 10 12 h -1 and 5.8 10 12 h -1 for 2MT and 3MT respectively) (Scheme 1 and Scheme 2). In addition, modeling confirmed that the direct formation of ethylbenzene from BT was a minor pathway in selective HDS operating conditions. Indeed, the rate constant is much lower (k = 3.8 10 8 h -1 ) than that via the formation of DHBT (k = 12 10 12 h -1 ) (Scheme 3). Moreover, the scale of reactivity of the three model sulfur molecules was confirmed from the rate constants of the desulfurization steps where those of BT is higher than those of 3MT and 2MT. Similarly, the BT adsorption constant is higher than those of 3MT and 2MT (K BT = 12.4 bar -1 , K 3MT = 10.6 bar -1 , K 2MT = 7.1 bar -1 ).

In addition, when components are in mixture, the experimental results are also well modeled considering the same model without further adjustments of the parameters fitted on single component experiments. Indeed, the same parameters (kinetic and adsorption constants) determined with the single sulfur compounds were used for the modeling in mixture. The very good fit obtained confirms the validity and the robustness of the chosen model and evidences the soundness of the obtained adsorption and kinetic parameters. The mutual inhibiting effect observed when 2MT and BT were in mixture is well represented. Indeed, the BT adsorption constant was higher (K ads = 12.4 bar -1 ) than that of 2MT (7.1 bar -1 ).

The transformation of 2MT (Figure 5) was thus impacted by the presence of unconverted BT, whereas conversely BT was little impacted by 2MT (Figure 6). This confirmed competitive adsorptions of the sulfur model molecules on the catalyst surface.

Conclusion

Under selective HDS operating conditions, the experimental approach developed in this study demonstrated an unambiguous reactivity scale of sulfur compounds representative of FCC gasoline. Thus, benzothiophene is 2.4 times more reactive than 3-methylthiophene which 1.8 times more reactive than 2-methylthiophene when considering single component experiments. By kinetic modeling based on a Langmuir-Hinshelwood formalism, the reaction schemes for the transformation of 2-methylthiophene, 3-methylthiophene and benzothiophene were validated.

In mixture, the reactivity ranking of these sulfur compounds was preserved. In addition, a mutual inhibiting effect between sulfur compounds was evidenced from the experiments mixing 2-methylthiophene and benzothiophene. Besides, the influence of unconverted BT on the transformation of 2MT was much higher than that of 2MT on the transformation of BT, due to a stronger adsorption of the latter. These results were confirmed by kinetic modeling where the calculated BT adsorption constant was 40% higher than that of 2MT, thus explaining the stronger adsorption of BT on the catalytic sites. This also allowed us to identify the most refractory sulfur compounds that are present in FCC gasoline feed under deep HDS operating conditions. Notably, the transformations of sulfur model compounds alone and in mixture could be fitted from a single kinetic model based on a Langmuir-Hinshelwood formalism. This clearly validated our combined experimental and theoretical approach and also demonstrated that inhibition phenomena observed experimentally between sulfur molecules were due to competitive adsorption on the catalyst surface.

To further narrow the representativeness gap between model molecules feeds and real catalytically cracked gasoline and get a deeper understanding of complex feed reactivities, a similar approach will be applied to understand the mutual influence of sulfur compounds and olefins on their reactivity.

Figure 1 .Figure 2 .

 12 Figure 1. Conversion of 2MT, 3MT and BT as a function of contact time. (T = 250ºC, P = 2 MPa, H 2 /feed = 360 NL/L, CoMoS/Al 2 O 3 ).

Figure 3 .

 3 Figure 3. Transformation of 3MT. Distribution of products as a function of 3MT conversion. (T = 250ºC, P = 2 MPa, CoMoS/Al 2 O 3 , H 2 /feed = 360 NL/L). 2MB2N: 2methylbut-2-ene, 2MB1N: 2-methylbut-1-ene, iC5: isopentane, 3MB1N: 3-methylbut-1ene.

Figure 4 .

 4 Figure 4. Transformation of BT. Distribution of products as a function of BT conversion. (T = 250ºC, P = 2 MPa, CoMoS/Al 2 O 3 , H 2 /feed = 360 NL/L). EB: ethylbenzene, DHBT: dihydrobenzothiophene.

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Transformation of 2MT: Comparison of simulation results (dotted lines) and experimental data (points) a) for the conversions of 2MT alone and in presence of BT and b) for the distribution of HDS products (pentane and pentenes). (T = 250ºC, P = 2 MPa, CoMoS/Al 2 O 3 , H 2 /feed = 360 NL/L).

Table 1 . Amount of sulfur compounds present in a feed containing 1.05 wt% of sulfur

 1

	[17].
		Sulfur Compounds
	Sulfur Compounds
		ppmS	%
	Mercaptans	34	4.5
	Thiophene	37	4.9
	C1-thiophenes	106	14.1
	Tetrahydrothiophene	24	3.2
	C2-thiophenes	118	15.6
	C3-thiophenes/thiophenol	76	10.1
	C4-thiophenes/C1-thiophenol	83	11.0
	Benzothiophene	276	36.6
	Total Sulfur	754	100

Table 2 : Characteristic of the CoMo/Al 2 O 3 catalyst

 2

	CoO (wt%)	MoO 3 (wt%)	S BET (m2/g)	Vp (cm3/g)	Dp (

Table 3 . Partial pressures (kPa) of the different compounds for the sulfidation step and the transformation of the different feeds.

 3

	Feed	A	B	C	D	E	F
	Pressure (kPa)	Sulfidation 3MT/2MT	BT	2MT + BT 2MT + H 2 S BT + H 2 S
	P H 2 S	10	0	0	0	0.8 -1.45 0.2 -0.37
	P 2MT or P 3MT	0	2	0	2	2	-
	P BT	0	0	2	1-3	0	2
	P H 2	90	1360	1360	1360	1360	1360
	P nC 7	0	638	638	633-635	637,2-636,55	637,8-637,63
	P TOT	100	2000	2000	2000	2000	2000
	(T = 250ºC, P = 2 MPa, CoMoS/Al 2 O 3 , H 2 /feed = 360 NL/L)

Table 4 . Transformation of sulfur compounds. Comparison of the activity

 4

		2MT	3MT	BT
	Conversion (%)	28.2	23.1	33.4
	Activity (mol.h -1 .g -1 )	0.376	0.597	0.833
	(iso-conversion around 30% from feed B or C) (T = 250ºC, P = 2 MPa, CoMoS/Al 2 O 3 ,
	H 2 /feed = 360 NL/L).

Table 5 : Transformation of 2MT. Effect of the amount of H 2 S or BT.

 5

					Inhibiting effect
	BT	BT Conv.	P H2S from BT
	(ppmS)	(%)	(kPa)
				Total (%)	By H 2 S (%)	By BT (%)
	500	80.0	0.80	63.1	1.7	61.4
	1000	60.3	1.20	68.3	2.8	65.5
	1500	47.0	1.45	70.3	3.8	66.5
	(conversion of 2MT of 42% from feed E and D, T = 250°C, P = 2 MPa, H 2 /feed = 360 NL/L,
	CoMoS/Al 2 O 3 ).

Table 6 . Transformation of BT. Effect of the amount of H 2 S or 2MT.

 6

				Inhibiting effect
	2MT Conv.	P H2S from 2MT
	(%)	(bar)
			Total (%)	By H 2 S (%)	By 2MT(%)
	10.0	0.2	17.2	0	17.2
	15.4	0.31	15.5	0	15.5
	18.9	0.37	14.2	0	14.2
	(conversion of BT of 30% from Feed D and F, T = 250 °C, P = 2 MPa, H 2 /feed = 360 NL/L,
	CoMoS/Al 2 O 3 ).

Table 7 . Kinetic parameters and adsorption parameters for the reaction of sulfur compounds, H 2 S and H 2 (HDS: Hydrodesulfurization, HYD: Hydrogenation)

 7

	(T = 250°C, P = 2 MPa, H 2 /feed = 360 NL/L, CoMoS/Al 2 O 3 )

Acknowledgements

A. Silva dos Santos thanks CNPq -Brasil for a PhD grant and IFPEN Lyon for the fruitful discussions.

Scheme 3. Transformation of benzothiophene (BT): Dihydrobenzothiophene (DHBT), Ethylbenzene (EB).


Rémi Manczak
email: q1remi.manczak@xlim.fr

Sofiane Saada
email: sofiane.saada@unilim.fr

Thomas Provent
email: thomas.provent@xlim.fr

Claire Dalmay
email: claire.dalmay@xlim.fr

Barbara Bessette
email: barbara.bessette@unilim.fr

Gaëlle Bégaud
email: gaelle.begaud@unilim.fr

Serge Battu
email: serge.battu@unilim.fr

Pierre Blondy
email: pierre.blondy@xlim.fr

Marie-Odile O Jauberteau
email: o.jauberteau-marchan@unilim.fr

Canan Baristiran Kaynak
email: kaynak@ihp-microelectronics.com

Mehmet Kaynak

Cristiano Palego
email: c.palego@bangor.ac.uk

Fabrice Lalloué
email: fabrice.lalloue@unilim.fr

Arnaud Pothier
email: arnaud.pothier@xlim.fr

F Lal- Loué

UHF-Dielectrophoresis Crossover Frequency as a New Marker for Discrimination of Glioblastoma Undifferentiated Cells

Keywords: BiCMOS chip, biological cell manipulation, 25 glioblastoma cells, high frequency dielectrophoresis. 26

published or not. The documents may come   L'archive ouverte pluridisciplinaire

according to Stupp protocol. Despite recent advances in surgery, imaging, radiation therapies and chemotherapy, the median survival is less than 15 months [START_REF] Chen | Malignant glioma: lessons from genomics, mouse models, and stem cells[END_REF]. This dark prognosis of GBM is primarily due to the recurrence of tumor, which is resistant to pre-cited conventional treatments [START_REF] Cheray | Cancer stem-like cells in glioblastoma[END_REF].

Limited advances in glioblastoma treatment are closely linked to the existence of a restricted cell subpopulation also called cancer stem cells (CSCs), some very immature and undifferentiated cells, responsible for tumor cell heterogeneity [START_REF] Bradshaw | Cancer stem cell hierarchy in glioblastoma multiforme[END_REF]. However, even genetically diverse clones express undifferentiated cell markers related to cancer stem cells such as CD133 and CD44. The higher expression levels of CD133 have been correlated to poorer prognosis suggesting that this marker might play a significant role in the resistance of this type of cancer to chemotherapy and radiotherapy [START_REF] Zeppernick | Stem cell marker CD133 affects clinical outcome in glioma patients[END_REF]. Other markers such as the transcription factors OCT-4, SOX2, pSTAT3, and NANOG are considered as key players in regulating transcription of glioblastoma CSCs [START_REF] Bradshaw | Cancer stem cell hierarchy in glioblastoma multiforme[END_REF]. These CSCs are a subpopulation of undifferentiated cells, which have specific biological properties similar to normal stem cells. Currently, biologists use some immunostaining approaches to characterize CSC populations, as flow cytometry, optical microscopy or protein array analysis, based on targeting a set of undifferentiated markers previously described. These markers are required to validate the stemness lineage of CSCs from the huge heterogeneous cell population. Nevertheless, CSCs subpopulation are very rare in tumors and their isolation often requires enriching them in specific culture medium. This strategy is time consuming and delays the results. Currently, the key objective is to try to get around this problem by establishing a new way to discriminate and sort undifferentiated cell populations specifically according to theirs biological and physical characteristics.

To optimize the diagnosis and prognosis methods, the development of different approaches and techniques based on bioelectric signals of cells have been proved to carry various helpful information on cell status [START_REF] Liu | Cell-based biosensors and their application in biomedicine[END_REF]. Many sources of cell bioelectric signals, like sodium potassium channels and pumps in the plasma membrane, may affect chemical analytes homeostasis, cell patterning and cell-to-cell interactions with the extracellular matrix, which can be determined by exploiting the dielectric properties. Among these techniques, Dielectrophoresis (DEP) is a label-free, accurate, fast, and low-cost analysis method that This approach offers unique capabilities to investigate the differences on the intracellular dielectric properties of each cell population [START_REF] Salimi | Dielectric model for Chinese hamster ovary cells obtained by dielectrophoresis cytometry[END_REF]- [START_REF] Ma | Ultra-wideband impedance spectroscopy of a live biological cell[END_REF] and allows screening of the intracellular biological properties and differences within the heterogeneity of a tumor.

II. BASICS ON CELLS ELECTROMANIPULATION BY DEP

When particles presenting different polarizabilities than the surrounding medium in which they are suspended are submitted to a non-uniform electric field, a DEP force is generated inducing motion of particles [START_REF] Pohl | The motion and precipitation of suspensoids in divergent electric fields[END_REF]. There are different ways to exploit this phenomenon. In the present case, a quadrupole microelectrode system [START_REF] Hjeij | UHF dielectrophoretic handling of individual biological cells using BiCMOS microfluidic RF-sensors[END_REF] has been used as sensor and implemented in a microfluidic channel (see Fig. 1). For such electrode geometry, the DEP theory [START_REF] Pethig | Dielectrophoresis: Theory, methodology and biological applications[END_REF], [START_REF] Fiorini | Disposable microfluidic devices: Fabrication, function, and application[END_REF], [START_REF] Broche | Extraction of dielectric properties of multiple populations from dielectrophoretic collection spectrum data[END_REF] shows that considering a cell as a homogeneous spherical dielectric particle, the induced DEP force can be then computed using equation [START_REF] Cheray | Cancer stem-like cells in glioblastoma[END_REF]. where r is the particle radius, ω is the angular frequency of the 112 applied electric field, E r m s is the root mean square value of this 113 electric field, is the gradient operator and Re[K(ω)] the real 114 part of Claussius-Mossotti factor K(ω) given by (2) in which ε p * 115 and ε m * refer to the complex permittivity of the particle and the 116 suspension medium, respectively. The complex permittivity ε x * 117 are defined in [START_REF] Bradshaw | Cancer stem cell hierarchy in glioblastoma multiforme[END_REF], where ε x and σ x are the relative permittivity 118 and conductivity either of particles or immersion medium, and 119 ε 0 represent the electric constant (8.854 10 -12 F m -1 ).

F D E P = 2πε m r 3 Re [K(ω)] ∇|E r m s | 2 (1)
K(ω) = ε * p -ε * m ε * p + 2ε * m (2) ε * x = ε 0 ε x -j σ x ω (3)
120

By changing the frequency of the applied electric field [START_REF] Pohl | The motion and precipitation of suspensoids in divergent electric fields[END_REF], 121 the polarized particles would behave in various ways depending 122 on the magnitude and the sign of Re[K(ω)] which is in turn 123 determined by the effective conductivity and permittivity of the 124 particle and the dielectric properties of the surrounding medium. 125 Therefore, particles can be individually electro-manipulated ac-126 cording to their own dielectric properties.

127

Actually, the generated force is repulsive when Re[K(ω)] is 128 negative, meaning that the particle is repealed away from elec-129 trodes (see Fig. 2(a) negative DEP case -nDEP). Whereas when 130 Re[K(ω)] is positive, the force is attractive and the particle 131 moves toward the electrodes where the electric field magnitude 132 is high (see Fig. 2(b) positive DEP case -pDEP). When the 133 force becomes null just before the cell switches to negative to 134 positive DEP (or vice versa), the DEP crossover frequency is 135 then reached. This frequency can be thus considered as charac-136 teristic of a cell own properties and specificities and may differ 137 between different cells.

Depending on the type of cell properties one wants to access, the choice of the DEP frequency range is important [START_REF] Du | Validation of clausius-mossotti function in single-cell dielectrophoresis[END_REF]- [START_REF] Afshar | Change in the dielectric response of single cells induced by nutrient deprivation over a wide frequency range[END_REF].

If information about cell plasma membrane specificities are sought, conventional DEP frequencies (typically from 100 kHz to 5 MHz) are very suitable for cell analysis. At this low frequency range, the cell shape, morphology and size have strong influence on the interaction with the electric field. Conversely, Ultra High Frequencies DEP (from 50 MHz to 500 MHz) will be better to provide information about intracellular properties. Indeed, when frequency increases above several tens of MHz, the plasma membrane lets the electric field penetrate the cell and interact directly with the cell interior. As a result, the effect of DEP forces generated (i.e., attractive vs. repulsive) at high frequency regimes may be different according to the dielectric properties of the cell content. This UHF-DEP characteristic is often presented in the literature [START_REF] Pethig | Dielectrophoresis: Theory, methodology and biological applications[END_REF], [START_REF] Gimsa | Dielectrophorsis and electrorotation of neurospora slime and murine myeloma cells[END_REF], [START_REF] Chung | Intracellular potassium under osmotic stress determines the dielectrophoresis cross-over frequency of murine myeloma cells in the MHz range[END_REF] and the crossover frequency f xo2 can be written by

f 2 x02 = 1 4π 2 1 ε 2 0 (σ m -σ p ) (σ p + 2σ m ) (ε p -ε m ) (ε p + 2ε m ) . ( 4
)
For an aqueous solution, considering ε m > ε p and σ p > σ m , (4) can be simplified with the approximation:

f xo2 = σ p 2π 1 ε 0 1 (ε p -ε m ) (ε p + 2ε m ) . (5)
The value of the crossover frequency f xo2 in the UHF range is directly influenced by the intracellular properties of the cells, largely by its conductivity and to a lesser extent by its permittivity. Specifically, since undifferentiated cells exhibit different biological specificities or physiological mechanisms linked to their differentiation state for example, their crossover frequencies will be different from those of differentiated cells. Therefore, the analysis of their dielectrophoresis behavior under UHF frequencies seems very relevant for the targeted application.

III. MATERIALS AND METHODS

A. Cellular Culture

Two GBM cell lines were tested in this study, U87-MG and LN18. Both of them derived from malignant stage IV gliomas from adult patients, purchased from the American Type Culture Collection (ATCC). Two conditions were used for the analysis: Actually, under stringent culture conditions, mainly only the most resistant cells with strong aggressiveness special features can survive and grow. As illustrated in Fig. 4, it is hence expected The main purpose of this study is to characterize GBM cell 198 lines to identify their DEP crossover frequencies in the high fre-199 quency regime and establish DEP signature according to their 200 different culture conditions (normal culture medium vs. define 201 medium). Each cell population is introduced into the microflu-202 idic chip, suspended in a DEP medium, by a fluidic inlet driven 203 by a flow controller (Fluigent MFCS) and flows in a Poly-204 diméthylsiloxane (PDMS) microfluidic channel implemented 205 above dedicated sensors implemented in BiCMOS technology 206 (see Fig. 1). The experiments were done using a 40 × 40 µm 207 gap quadrupole electrodes design. This structure is based on 208 four electrodes, set at 90°, combining a pair of thick (9 µm) 209 electrodes crossing the microfluidic channel with another pair 210 of thin (0.45 µm) electrodes implemented in the middle of the 211 channel [START_REF] Hjeij | UHF dielectrophoretic handling of individual biological cells using BiCMOS microfluidic RF-sensors[END_REF].

(i) normal differentiation

212

The selected 40 µm spacing between each electrodes rep-213 resent a good compromise between an easy monitoring under 214 microscope of cell motion submitted to both positive or negative 215 DEP forces and the use of moderate RF voltage signal to bias 216 the structure and efficiently act on cells (typically the magnitude 217 of the applied voltage ranges between 2 and 4 V pp ). The same 218 frequency adjustable DEP signal has been applied to the left and 219 right electrodes whereas top and bottom ones were grounded. 220 The flow is slowed down, and when a cell arrives near the cen-221 ter of quadrupole electrodes system, the electrodes were biased 222 with a 500 MHz DEP signal expected to be much higher than 223 the f xo2 and therefore suitable to efficiently trap cells. This methodology of crossover frequency measurement has been supported and validated though electro-kinetic transient simulations using COMSOL software. Hence, the displacement of the trapped cell, from the center of quadrupole to the edge of the electrodes, according to the decrease of the DEP signals frequency, can be computed as presented in the graph below.

As Fig. 7 shows, the simulations predict that the generated DEP force, once applying a 4V pp DEP signal with a frequency set to 1.5 MHz under the cell crossover frequency, should have a sufficient influence to attract and move it 6 µm away from its initial negative DEP trapping location. Actually, such displacement value matches well with the ones observed under a microscope during experiments, as illustrated on the Fig. 6 photograph taken at 118 MHz. Consequently, we can reasonably consider that using the proposed cell characterization methodology the cell crossover frequency can be measured with a 1-2 MHz accuracy.

IV. RESULTS, ANALYSIS AND CORRELATION

A. GBM Cell Lines Phenotypic Profiles

First, control experiments were carried out to confirm the enrichment of undifferentiated cell population in define medium. Comparative analysis of the gene expression (mRNA levels) of the stemness lineage was assessed in the cells cultured 6 days in normal culture medium vs. define medium (see Fig. 8). Analyzed CSCs markers (Oct-4, Sox2 and Nanog) showed an overexpression in cells cultured in define medium (in red) compared to those cultured in normal medium (in blue).

These results were completed by analyzing protein levels with flow cytometry. They showed an enhancement of undifferentiated markers expression in both U87-MG (Fig. 9(a) and (b)) and LN18 (see Fig. 9(c) and (d)) cell lines cultured in define medium compared to normal medium culture. These biological results validate the enrichment of CSC by culture in define medium.

B. GBM Cell Crossover Frequencies

In order to support these results obtained at the biological level, the two GBM cell lines were characterized by establishing their UHF-DEP crossover frequencies according to their different culture conditions.

The set of crossover frequencies measured for both cell lines (U87-MG and LN18) cultured in the two different conditions (normal culture medium vs. define medium) is represented in the Fig. 10 (one should notice that the middle bar here represents the median value of the whole collected data). The considered crossover frequency corresponds to the frequency for which the trapped cell just starts to move away from the center of electrodes quadrupole. Based on statistical analysis of the results observed in the four populations, culturing cells in a define environment seems to have a real impact on the measured crossover frequencies, according to GBM intracellular characteristics changes. The set of statistics concerning the characterization of U87-MG and LN18 cells crossover frequencies is summarized in Table I, listing the average, median, standard deviation, standard error, minimum and maximum crossover frequency values for each cell population.

Fig. 1 .

 1 Fig. 1. Quadrupole electrodes system implemented in BiCMOS back end of line SG25H4 technology from IHP.

Fig. 2 .

 2 Fig. 2. DEP single cell manipulation principle. (a) Cell repulsion at system center with nDEP. (b) Cell attraction toward electrode with pDEP.

Fig. 3 .

 3 Fig. 3. Typical DEP spectral signature of cell with its two crossover frequencies f x o 1 and f x o 2 , respectively, at low and high frequency regimes.

Fig. 4 .

 4 Fig. 4. Representative diagram of the two different culture conditions used for the both GBM cell lines.

195B.

  Tools and Methodology for Cell Crossover 196 Frequency Measurement 197

Fig. 5 .

 5 Fig. 5. Cells suspension flowing in the microfluidic system (a). Single cell nDEP trapping at 500 MHz in the center of the quadrupole electrodes related to the generated electric field (b).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Cultured in NN -Microscope imaging sequence of GBM cell crossover frequency measurement (119 MHz) by tuning the DEP signal frequency.

Fig. 8 .

 8 Fig. 8. Comparative analysis of gene expression of three undifferentiated markers: Oct-4, Sox2 and Nanog, in U87-MG and LN18 cell lines, cultured 6 days in normal medium (NN: blue histograms) or in define medium (DN: red histograms), measured by Real Time PCR (Polymerase Chain Reaction).

Fig. 9 .

 9 Fig. 9. Comparative analysis of the undifferentiation markers expression CD44, CD133, Oct-4, Sox2 in both U87-MG and Ln18 cell lines, cultured 6 days in normal medium (NN: blue histograms) or in define medium (DN: red histograms), analyzed by multi-parameteric flow cytometry (BD Fortessa). (a) and (c) graphs represent percentage of expression and (b) and (d) histograms represent the mean fluorescent intensity (MFI) for each marker expression in U87-MG and LN18 respectively, grey graphs show the isotypes (unlabeled) negative control condition.

Fig. 10 .

 10 Fig. 10. Graphic box plots representation of U87-MG and LN18 cells crossover frequencies, cultured in two different conditions: normal medium (NN) and define medium (DN). The p value was determined using One-way ANOVA test. * * * represents p value <.0001, * * represents p value <.001.

  As shown a significant number of cells have been characterized showing statistically consolidated data. The large standard deviation, and error standard, observed for U87-MG and LN18 NN cell pool can be explained by the natural cell line heterogeneity, including a large number of different differenti-313 ated cells but also some few undifferentiated occurrence in the 314 pool. On the other hand, the DN cell pools may concentrate 315 a much higher number of undifferentiated and low differenti-316 ated cells; since DN culture conditions are not favorable for 317 differentiated cell growth. 318 Considering the measured crossover frequencies, these two 319 GBM cell lines exhibit the same behavior according to the two 320 different culture conditions. Actually, the undifferentiated en-321 riched populations (DN) show much lower crossover frequen-322 cies than the cells cultured in normal conditions, although some 323 crossover frequencies overlap exist between these two popu-324 lations -U87-MG: Average of 120 MHz for NN vs. 91 MHz 325 for DN -LN18: Average of 128 MHz for NN vs. 76 MHz 326 for DN. This decrease demonstrates a significant difference (il-327 lustrated by the run ANOVA statistical analysis tests resulting 328 in a very low p value) between these two population profiles 329 obtained by different culture conditions. This finding proves a real difference on the intracellular dielectric characteristics of 331 the undifferentiated cells enriched populations compared to dif-332 ferentiated cells, reflecting the intrinsec biological properties 333 differences. 334 This difference outlines a great potential for discrimination 335 of cell subpopulations within the whole tumor mass. Therefore 336 such technique is highly promising to achieve discrimination 337 and even isolation of undifferentiated cells allowing potential 338 cell sorting of these undifferentiated subpopulations related to 339 the CSC subpopulation.

341

  We described here a novel method of cellular subpopula-342 tions' discrimination, which completes the classical biological 343 approaches, based on the differential expression levels of a set 344 of markers. These populations, with different cellular differen-345 tiation status, are discriminated using real time measurement 346 on microfluidic lab-on-chip (LOC) platform implemented in 347 CMOS technology. Both selected GBM cell lines, showed a 348 strong correlation between the biological markers differences 349 and the measured DEP frequency signatures according to the 350 different culture conditions. Differences on crossover frequen-351 cies obtained for each subpopulation, showed a great discrimina-352 tion potential especially for the development of a novel method 353 to characterize Cancer Stem Cells. Thus, we confirmed the bi-354 ological differences analyzed by routine methods, using DEP 355 signatures differences, which complete the characterization of 356 stemness properties cells. These results correlate to the bio-357 logical differences at the functional level. The undifferentiated 358 properties of CSCs are associated to intracellular changes and re-359 flecting their high aggressiveness potential. This technic allows 360 screening of a new cell discrimination parameter, the intracel-361 lular differences and physical properties of cells, without any 362 labeling, without affecting cell integrity and viability. Hence, 363 based on the UHF-DEP spectroscopy method, we detected and 364 characterized the undifferentiated cells with unique capabilities 365 to screen biological specificities by investigating the intracellu-366 lar dielectric properties.

Finally, this method

  confirms a high potential of emerging lab-on-chip (LOC) platforms in the diagnosis and the treatment of glioblastoma. physics and material sciences from the University 441 of Toulouse, Toulouse, France, in 2016. Since 2017, Gaëlle Bégaud is a pharmacist who received the M.S. degree in molecular, cellular, and integrated physiopathology from the University of Toulouse, France, in 2003 and the Ph.D. degree in analytical chemistry and cellular biology from the University of Limoges, France, in 2007. In 2008, she joined Limoges University, where she is now an assistant professor in analytical chemistry and a researcher in the EA3842, CAPTuR team. Her research activity linked to the SUMCASTEC project is centered on the optimizaton of the use of Sedimentation Field Flow Fractionation (SdFFF) to sort Brain Cancer Stem Cells (BTSC) from complex populations such as neuroblastoma or glioblastoma cell lines. Her research interests include the implementation of new SdFFF prototypes for biological applications. Serge Battu is a pharmacist who received the M.S. degree in Chemistry from the University of Toulouse, France, in 1993, and the Ph.D. degree in Cellular and Molecular Biology from the University of Limoges, France, in 1997. In 1998 he joined Limoges University, where he is now a full professor in Analytical Chemistry and instrumentation and a researcher in the EA3842 team. His research research activity is centred on development of news SdFFF prototypes and to use SdFFF to sort cancer stem cells (BTSC) from complex populations such as neuroblastoma or glioblastoma cell line. Pierre Blondy received the Ph.D. and Habilitation degrees from the University of Limoges, Limoges, France, in 1998 and 2003, respectively. He joined the Centre National de la Recherche Scientifique as a research engineer in 1998 and the University of Limoges in 2006, where he is currently a Professor. From 2011 to 2015, he held an Institut Universitaire de France research chair. He was a visiting researcher at the University of Michigan, Ann Arbor, USA, in 1997 and at the University of California at San Diego, La Jolla, USA, in 2006 and 2008. He founded a research group in 2003 at XLIM on innovative technologies for microwave and millimeter-wave applications, and he was involved in many research projects, funded by national and european public agencies and companies. He has been working extensively on tunable filters, RF-MEMS components, and reliability. His group is also studying the properties of phase change transition materials and their applications to microwaves. Switches, power limiters, and tunable filters have been developped using innovative materials like Vanadium Dioxide and Germanium Telluride. Dr. Blondy was an Associate Editor for the IEEE MICROWAVE AND WIRELESS COMPONENTS LETTERS in 2006. He has been a member of the IEEE International Microwave Conference Technical Program Committee since 2003, and he is also serving on the technical program of the European Microwave Conference. He is the past chair of the MTT-S Technical Committee 21 on RF-MEMS. Marie-Odile Jauberteau received the Medical Doctor degree from Limoges University and the Ph.D. degree from the University of Paris VI. She has been a Medical professor of Immunology since 2000 at Limoges University, Faculty of Medicine. She is the director of the Research Team EA 3842, Control of Cell Activation, Tumor Progression and Therapeutic Resistance (CAPTuR), created in 2004. Her research is focused on oncology, especially mechanisms of cell survival and cancer stem cells. Canan Baristiran Kaynak received the B.S degree from the Chemistry department of Istanbul Technical University and the master degree in Material Science and Engineering at Sabanci University, Istanbul, Turkey, in 2007. She received the Dr.-Ing. degree from the institute of high frequency and semiconductor system technologies from Technical University of Berlin, Germany, in 2012. Since 2008, she has held a position as scientist in the Technology Department at IHP Microelectronics. She is leading projects including heterogeneous integration of silicon-based microfluidics and microbolometer devices into BiCMOS technology. Mehmet Kaynak received the B.S degree from the Electronics and Communication Engineering Department of Istanbul Technical University (ITU) in 2004, the M.S degree from the Microelectronic program of Sabanci University, Istanbul, Turkey in 2006, and the Ph.D. degree from Technical University of Berlin, Berlin, Germany, in 2014. He joined the technology group of IHP Microelectronics, Frankfurt (Oder), Germany, in 2008. From 2008 to 2015, he led the MEMS development at IHP. Since 2015, he has been the department head of Technology at IHP. Dr. Kaynak is being affiliated as Adjunct Professor at Sabanci University, Turkey. Cristiano Palego received the M.S. degree in Electrical Engineering from the University of Perugia, Italy, in 2003 and the Ph.D. degree in Microwave Engineering and Optoelectronics from the University of Limoges, France, in 2007. He was a postdoctoral associate then a research scientist at Lehigh University, Bethlehem, PA, USA, (2007-2017). In 2013, he joined Bangor University, U.K., where he is now an associate professor in smart sensors and instrumentation and a researcher in the Medical Microwave Systems group. His research interests include micro/nanotechnology, minimally invasive animal telemetry, and biomedical technology, along with energy harvesting and RF-MEMS for reconfigurable antenna arrays. Fabrice Lalloué received the Ph.D. degree in neuro-599 biology and neurodevelopment from the University 600 of Limoges, France, in 2001. He is currently Profes-601 sor for the University of Limoges in the Department 602 of Physiology and Neurosciences (Neuro-Oncology). 603 He is also assistant director of the research team 604 EA3842 CAPTuR (Control of cell Activation, Pro-605 gression of Tumor and therapeutic Resistance) at 606 the Faculty of Medicine. In 2000, he was visiting 607 the Institute of Development, Aging and Cancer in 608 Sendai, Japan. His research activities currently ad-609 dress Neuro-Oncology, Glioblastoma, Brain Tumour, Initiating stem cells, adult 610 neural stem cells, and Tumour Biomarkers. He is involved in the Scientific Pilot 611 Committee of the Canceropole Grand Sud-Ouest consortium. He has authored 612 or co-authored over 50 publications in research articles and conferences and 613 holds seven patents.

Arnaud

  Pothier received the Ph.D. degree in elec-616 trical engineering from the University of Limoges, 617 Limoges, France, in 2003. He is currently a full-618 time senior researcher with the XLIM research in-619 stitute. His current research activity is focused on RF 620 biosensors for label-free cell analysis. Dr. Pothier is a 621 member of the IEEE and EuMA society, member of 622 Technical Program Committees of IMS, EuMW, and 623 IMBIOC Conferences, and previously served as chair 624 of the IEEE MTT-S Technical Committee MTT-10 on 625 Biological Effect and Medical Applications.

TABLE I SUMMARY

 I

OF CROSSOVER FREQUENCY MEASUREMENTS (MHZ)

Rémi Manczak received the Ph.D. degree in

This work was supported by the European Union's Horizon 2020 Research and Innovation Program under Grant 737164.