Chunfeng Yang
email: chunfeng.yang@seu.edu.com

Regine Le

Bouquin Jeannes
email: regine.le-bouquin-jeannes@univ-rennes1.fr

Jean-Jacques J Bellanger
email: jean-jacques.bellanger@univ-rennes1.fr

Huazhong

Bouquin Jeannès

H Shu

A new

Keywords: Bayesian information criterion, causality, EEG signal, physiology-based model, transfer entropy

I. INTRODUCTION

PILEPSY is a n eurological disorder characterized by repetitive seizures. In 30% of cases, seizures remain drug-resistant and considerably affect all aspects of the patient's life [START_REF] Engel | Outcome with respect to epileptic seizures[END_REF]. Drug-resistant epilepsies are often partial, with an epileptogenic zone (EZ) located in a r elatively circumscribed brain area. For these partial epilepsies, surgical treatment is considered. The difficulty that arises is to determine the organization of the EZ and, thus, the part of the brain that should be excised to remove seizures without affecting normal brain function. In some patients, the pre-surgical evaluation may include recording of intracerebral ElectroEncephaloGraphic (iEEG) signals using depth electrodes. The analysis of such signals which remains a difficult task is aimed at determining which sites of the brain belong to the EZ, prior to surgery. In this context, signal processing techniques provide quantitative information that cannot be easily obtained by visual inspection. This is typically the case of correlation (in wide-sense) measures that prove useful for assessment of functional couplings between remote brain sites. Beside functional connectivity, effective connectivity provides complementary information on the directionality of couplings, i.e. the influence one neural population exerts on another one. Recently, a t heoretical information measure, namely Transfer Entropy (TE), was proposed to identify the direction of the information flow and to quantify the strength of coupling between complex systems [START_REF] Schreiber | Measuring information transfer[END_REF]. This "model-free" technique is based on t he transition probabilities between states of the considered systems from output signals. This method has been applied and tested in some nonlinear benchmark models, real EEG signals and magnetoencephalographic data in cortical and cerebellar networks with no gr ound reference [START_REF] Sabesan | Information flow in coupled nonlinear systems: application to the epileptic human brain[END_REF]- [START_REF] Wibral | Transfer entropy in magnetoencephalographic data: Quantifying information flow in cortical and cerebellar networks[END_REF]. Moreover, from these foregoing literatures, certain calibration parameters involved in the TE estimation play a significant role in acquiring the correct information flow direction between the two systems under investigation, i.e. the driven system and the driving system, for instance, the orders of the underlying Markov processes. In [START_REF] Schreiber | Measuring information transfer[END_REF], only the Markov process order of the driving system is mentioned: it is fixed to that of the driven one or set to 1 with the reason that the latter is preferable for computational reasons. In [START_REF] Sabesan | Information flow in coupled nonlinear systems: application to the epileptic human brain[END_REF], [START_REF] Sabesan | Information flow and application to epileptogenic focus localization from intracranial EEG[END_REF], Sabesan et E investigated this problem. Two measures, the delayed mutual information and the autocorrelation function [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF], [START_REF] Martinerie | Mutual information, strange attractors, and the optimal estimation of dimension[END_REF], were proposed to estimate the Markov process order of the driven system. However, the Markov process order of the driving system was still set to 1 since the authors presumed that the current state of the driving system was sufficient to arouse a considerable change in the dynamics of the driven system with one time step in the future. Furthermore, Vicente et al. [START_REF] Vicente | Transfer entropy -a model-free measure of effective connectivity for the neurosciences[END_REF] and Wibral et al. [START_REF] Wibral | Transfer entropy in magnetoencephalographic data: Quantifying information flow in cortical and cerebellar networks[END_REF] used Cao's criterion [START_REF] Cao | Practical method for determining the minimum embedding dimension of a scalar time series[END_REF] based on the false neighbor computation and autocorrelation function to estimate the Markov process orders. Two other model selection methods, the Akaike Information Criterion (AIC) [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] and the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a m odel[END_REF], have been also applied in Vectorial AutoRegressive (VAR) model order estimation and to estimate the order of Markov chains [START_REF] Katz | On some criteria for estimating the order of a Markov chain[END_REF]. To compute transfer entropy on electrophysiological signals, AIC and BIC are currently used [START_REF] Ding | Granger causality: basic theory and application to neuroscience[END_REF] to select the orders in a VAR model (including N channels 1 , ..., N XX ) under the constraint that, when regressing j X on i X , the same number = ij qq of lags is used for any ( , ij ) (where q must be determined and ij q is named hereafter the regression order from i X to j X even if no explicit regression procedure is introduced to compute TE). In this case, this procedure returns only one value for the number of past samples retained on X and on Y to compute TE. In many cases, the "true" regression orders for the driven and driving systems may be different (this is typically the case for the cross regression order when there is no influence from one system to another). Consequently, it can lead to erroneous estimation in TE when the regression orders for the driven and driving systems are chosen artificially identical. To deal with this issue, we released the constraint ( , ):

∀= ij ij q q by applying the generalized Akaike information criterion (gAIC) or the generalized Bayesian information criterion (gBIC). In Section II, after recalling the transfer entropy approach, we develop the methodology to estimate the regression orders for the driven and driving systems. We also introduce a transfer entropy based measure (called PTE) to quantify the contribution of an observation to another one and we give some details on the physiology-based model used as ground-truth. Section III is devoted to experimental results: two linear AR models were first introduced to verify the relevance of the two new order-estimated measures. Then, TE w as tested on a physiological model and on real iEEG signals involving the new criteria. Results are discussed and a conclusion is drawn in Section IV.

II. METHODS AND MATERIALS

A. Transfer Entropy

Let ( )

- -+  =  . (1)
Considering a th k -order Markov process X, the conditional probability measures (continuous or discrete) are such that:

( ) ( ) ( ) ( ) ( ) ( ) 11 // ': kk k k nnn nnn nn XX x XX x k k P dx P dx ++ ++ = = ∀> = .(2)
The relation in ( 2) can be extended by considering two random processes X and Y . The absence of information flow from Y to X is then formalized by x to the state value

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1 n x + ( 1 + n X and ( ) l n Y are independent conditionally to ( ) k n X ).
The deviation from this assumption can be quantified using the Kullback entropy, which leads to the definition of the transfer entropy where the ratio in (4) corresponds to a derivative of the conditional measure in the numerator with respect to the conditional measure in the denominator. This measure is not symmetric ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
XY Y X TE TE →→ ≠
. The values of the parameters k and l are the regression orders to predict

( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 1 , ˆ,, , 1,..
In practice, the estimations of the conditional probabilities in (6) are computed from the estimated joint probabilities in [START_REF] Vicente | Transfer entropy -a model-free measure of effective connectivity for the neurosciences[END_REF] using Bayes' rule ( ) ( ) ( ) | ,/ = PAB PAB PB . Thereby, the definition in ( 6) is rewritten as

 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1
→ + + + = × ∑ . ( 7
)
Considering the more complicated issue of real data [START_REF] Schnupp | Plasticity of temporal pattern codes for vocalization stimuli in primary auditory cortex[END_REF], we define a selective index of the flow direction, named PTE, by

→ → →→ = + YX YX XY Y X TE PTE TE TE , (8)

PTE

) reveals that no information flow exists from Y to X (resp. from X to Y ). Clearly, when estimations of PTE are substituted for theoretical PTE, some caution must be introduced. In practice, the selective index PTE is calculated only when the variables are first proven not to be independent.

B. Improved Computation of Transfer Entropy

• Selection of the parameters k and l For the two processes X and Y , the selection of the regression orders k and l plays experimentally a primary role to acquire dependable values for transfer entropy. As suggested before, to get the models' orders, AIC and BIC are often applied under the hypothesis that k and l are identical, which not only produces extra computation but also can result in erroneous estimation in TE. , are × mm coefficient matrices and q the model order. The process w n is a zero mean independent identically distributed vector process with covariance matrix Σ . It is assumed that w n and z - ni are independent for each variable 0 > i . In addition, it is assumed that z n is mean and covariance ergodic.

Suppose that the observations 12 z , z ,..., z N are generated by the VAR process given in [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. The least-squares method [START_REF] Brockwell | Time Series: Theory and Methods[END_REF] is used to fit a VAR model with order q to the observations. Suppose that the estimated model is as follows:

    12 12 z z z zw -- - =Φ +Φ + +Φ +  q n n n n nq . ( 11
)
The coefficient matrices  , 1, 2,..., Φ= i iq , of this estimated model are obtained by solving the following set of equations:

   =           TT T qq TT T qq TT T N N Nq ( ) ( ) ( ) 1 2
,:

,: A, ,:

 Φ   Φ =    Φ   T T j T q j j j and ( ) ( ) ( ) 1 2 1, z z z z + + +    =      q q qj N j j j where ( ) ,: Φ i j
is the th j row of Φ i . The least-squares solution of ( 12) is given by  ( )

1 1, A Z Z Z z ; 1, 2,..., - + = = TT j qj jm . (13)
Then, the residual covariance matrix of the VAR model is estimated as follows:

 ( )( ) 1 1 zzzz = + Σ = - - - ∑  N T nn nn nq Nq ( 14
)
where

   12 12 zz z z -- - =Φ +Φ + +Φ  
nq n n nq . To select the order parameter q, the Akaike Information Criterion is given by:

( )  ( ) ( ) 2 AIC ln det 2 = Σ+ q N mq, (15)
and the Bayesian Information Criterion is:

( )  ( ) ( ) ( ) 2 BIC ln det ln = Σ+ q N mq N . (16)
Rewriting [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]

with a 2-dimensional ( 2 = m ) VAR process ( ) ( ) [ ] z 1, 2 , = =   T T n n n nn z z xy , we have ( ) ( ) ( ) ( ) , 11 , 11 α α α α -- = = -- = =  =++     =++   ∑∑ ∑∑ yx xx xy yy q q n xx n i yx n i n x ii qq n xy n i yy n i n y ii x ix iy w y ix iy w . ( 17
)
where α are the coefficients of the VAR process. AIC and BIC can be used to determine a common summation range for the four sums in ( 17

TE

. In fact, for most of physical cases, these four orders are different, and the standard approach (AIC or BIC) only gives the maximal value among them. Moreover, if we choose the maximal value for the parameters k and l , TE computational time is drastically increased. Focusing on t his point, we propose to minimize the following criteria, namely gAIC and gBIC, with respect to the four orders

{ }

,,,

xx yx xy yy qqqq ( )  ( ) ( ) ( ) gAIC ln det 2 = Σ+ + + +  xx yx xy yy qN q q q q, (18) ( )  ( ) ( ) ( ) ( ) gBIC ln det ln = Σ+ + + +  xx yx xy yy qN q q q q N , ( 19
)
where  q is a set of estimated orders { }

,,,  xx yx xy yy qqqq

. Using gAIC and gBIC criteria is justified by considering that AIC and BIC criteria have been developed under less restrictive constraints [START_REF] Schwarz | Estimating the dimension of a m odel[END_REF], [START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF], [START_REF] Neath | Regression and time series model selection using variants of the Schwarz information criterion[END_REF] (i.e. not only for linear models such as VAR models). In this widened context, the log-likelihood function (which corresponds in ( 16) to the first term on the right hand side) must only belong to an exponential family, and the tested models do not need to be nested. When applied to VAR order selection, even if the search of candidate models may be limited to the subclass defined by xx yx xy yy qqqq = = = , as it is generally done in most applications, this constraint can be interestingly released to expect a better fitting without losing optimality relatively to the type of criterion.

Using gAIC and gBIC to get the four orders, the transfer entropy computation time is reduced, even though the computation time to evaluate the four orders is increased compared to the time needed to compute the orders with AIC and BIC. For example, given

[ ] 1, 20 ∈ q
, procedures ( 12) to [START_REF] Schnupp | Plasticity of temporal pattern codes for vocalization stimuli in primary auditory cortex[END_REF] merely need to be implemented 20 times to minimize AIC and BIC. For gAIC and gBIC, this step needs to be implemented

21 4 -1 times since { } [ ] 4 , , , 0, 20 ∈
xx yx xy yy qqqq respectively and, therefore, it becomes prohibitive. To solve this problem, we propose a "greedy" strategy [START_REF] Shu | Moment-based methods for polygonal approximation of digitized curves[END_REF] to optimize the computation. This procedure using gAIC is detailed below (let us note that this procedure would be the same using gBIC): 1) Use AIC to obtain a model order, namely the maximum order max q , and select this order to decide the range for gAIC, i.e.

[ ]

max 0, ⊂ ij qq , { } ,, ⊂ ij xy ; 2)
Considering the values of ij q for the first signal (i.e. xx q and yx q ), while setting the values of ij q for the second signal equal to max q (i.e. xy q and yy q ), implement the procedures from ( 12) to ( 14) to get the optimal orders  xx q and  yx q using the definition of (18); 3) Considering the values of ij q for the second signal (i.e.

xy q and yy q ), while setting the values of ij q for the first signal equal to the optimal order  xx q and  yx q , implement the procedures from ( 12) to [START_REF] Schnupp | Plasticity of temporal pattern codes for vocalization stimuli in primary auditory cortex[END_REF] to get the optimal orders  xy q and  yy q using the definition of [START_REF] Shu | Moment-based methods for polygonal approximation of digitized curves[END_REF]. According to this "greedy" strategy, the computation time is decreased significantly. As an example, the global procedure only needs to be implemented 2 [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF] 16 2 1 + ×-times (0.27% of the computation of AIC) when max 15 = q .

• Selection of the Radius r For the selection of the radius r , we take the similar notion as in [START_REF] Sabesan | Information flow in coupled nonlinear systems: application to the epileptic human brain[END_REF], [START_REF] Sabesan | Information flow and application to epileptogenic focus localization from intracranial EEG[END_REF], [START_REF] Grassberger | Measuring the strangeness of strange attractors[END_REF]. Firstly, the data are normalized to zero mean and unit variance in order to establish a common radius r in the state space X and Y . Then, the joint probability in log scale (

( )

ln Crobtained as an average on n over , nr C in ( 5)) is calculated and plotted vs the corresponding radius r in log scale (e.g. in Fig. 2(a)). Finally, some optimal value * r is chosen in the linear region of the curve ( ) ln Cr vs ln r so as to estimate reliable values of TE at specific values * r .

C. Model of iEEG Signals Generation

We used a physiology-based time continuous model to represent the electrical field activity of some distant -and possibly coupled -neuronal populations. Each of them generates a local field activity that is converted to an iEEG signal in a proximal electrode using a quasi-static transfer function [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. In the model, each population contains three subpopulations of neurons that mutually interact via excitatory or inhibitory feedback linking main pyramidal cells and two other types of local interneurons. Since pyramidal cells are excitatory neurons that project their axons to other areas of the brain, the model accounts for this organization by using the average pulse rate of action potentials from the main cells of one population i as an excitatory input to main cells inputs of another population j. In addition, this connection from population i to population j is represented by a parameter ij K which is proportional to the number of corresponding axonal links. Appropriate setting of this parameter allows for building systems where the neuronal populations are unidirectionally or bidirectionally coupled. Other introduced parameters are intra-population parameters. They include excitatory and inhibitory gains in feedback loops as well as coefficients related to the number of synaptic contacts between subpopulations. These parameters are adjusted to control the intrinsic activity of each population (normal background activity vs epileptic activity). This population model [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]  defined as a smooth non l inear sigmoid function f . This type of function is necessary to take into account effects of thresholding and saturation in neurons [START_REF] Wendling | Interictal to ictal transition in human temporal lobe epilepsy: insights from a computational model of intracerebral EEG[END_REF]. In [START_REF] Yang | Extracting information on flow direction in multivariate time series[END_REF], D is an 11-dimensional row vector and the populations' outputs (which are not directly measured) are obtained through the sigmoid function. The time discretization of ( 20) and ( 21), needed to get a computational model, cannot lead to a linear representation as in [START_REF] Neath | Regression and time series model selection using variants of the Schwarz information criterion[END_REF] and so the resulting model induces nonlinear dynamics. Finally, both observations x and y are obtained as a l inear combination of three of the eleven state variables.

D. Real Signals

A. Linear Models

For the first linear stochastic system we considered, the following signals with unidirectional relation were generated:

 = -+   =-+   t t tt t t tt x x xu y x yv ( 23
)
where t u and t v were independent white Gaussian noises with zero means and unit variances. To compare the performance of the 4 criteria (AIC, BIC, gAIC, and gBIC), we indicated in Table I the number of times the orders obtained with these 4 measures equal the expected optimal orders over the 200 trials.

As discussed in Section II. B, the expected optimal order values of the model described by ( 23 I). So, in this experiment, 3 parameters out of 4 a re not correctly estimated. The gAIC measure is a reasonable estimator since most of the estimated orders are the optimal ones (for example, 189 correct estimations over a total of 200 (94.5%) for xx q ). In this experiment, the gBIC measure behaves the best (100% of accuracy estimation for the 4 parameters).

Given the orders estimated by the 4 techniques, we computed the means of the corresponding TE values over the 200 simulations. Results are shown in Fig. 1 for different values of the radius r . Using gBIC,  → XY TE is the greatest in some range of ln r (around the vertical dotted line) when there is effective information flow from x to y while  → YX TE is zero when there is no information flow from y to x (see Fig. due to some erroneous model orders: the estimated orders k () = xx q and l () = yx q are always equal to 3 which are not the optimal orders. Using gAIC, the transfer entropy from x to y is close to the same measure using gBIC, while the reciprocal quantity  → YX TE is mostly close to zero (see Fig. 1(c)) thanks to only a few incorrect estimated orders. For each technique, to ensure a correct analysis of the information flow direction, a s pecific radius * r is chosen in the linear region of the curve To test the reliability of the new criteria, some parameters values of our linear models were modified and the corresponding effects summarized in Annex A. Finally, to complete these results and test the relevance of the approach, an extended range of model orders values around the best ones chosen by the algorithm was tested and subsequent impacts on TE estimation are reported in Annex B.

From this table, we observe that

For the second linear stochastic system, the following signals with bidirectional relation were generated: where t u and t v were independent white Gaussian noises with zero means and unit variances. The numbers of trials (over 200 trials) for which the model orders obtained by the 4 criteria equal the expected optimal order are presented in Table III. Given these orders, TE is computed on 200 t rials and the averaged values are displayed in Fig. 3 when the radius r in log scale is in the range [-1.5; 2]. In the same way as previously, a specific radius * r is retained, and the means and standard deviations (shown in parentheses) of the corresponding TE and PTE values are given in Table IV IV).

If we focus on the PTE index (see Tables II andIV), both PTE(AIC) and PTE(BIC) conclude to a bidirectional information flow whatever the actual situation (unidirectional or bidirectional flow) whereas the two other indices PTE(gAIC) and PTE(gBIC) allow us to illustrate the two different situations: unidirectional vs bidirectional flow. For example, when the propagation is unidirectional (see Table II

B. Physiology-based Models

The model described in Section II.C was used to simulate long duration signals (800 s) for a fixed connectivity pattern ("ground truth") among neuronal populations. Sampling rate was equal to 512 Hz. Two simulated physiology-based models with uni-and bi-directionally coupled connectivities are considered. Model parameters were such that a fast quasi-sinusoidal (25 Hz) activity (similar to that observed during seizure onset) was generated by the two populations when they were uni-or bi-directionally coupled. The coupling parameters were such as: 12 In this section, according to the results obtained on t he simulated linear models, we only tested the performance of the gBIC technique and compared it to the BIC. Given the orders estimated by gBIC, we computed the corresponding TE values and averaged them on 200 trials for different values of the radius r as shown in Fig. 4. As previously, a specific radius * r is chosen for each pair [START_REF] Engel | Outcome with respect to epileptic seizures[END_REF][START_REF] Schreiber | Measuring information transfer[END_REF], and the means and standard deviations (in parentheses) of the corresponding TE and PTE values are displayed in Table V (gBIC) are close to 0.5.

To confirm these results, we conducted a statistical test. Two independent populations were first simulated (using 12 21

= = KK

). Then, we computed TE with the order-estimated methods of BIC and gBIC when these populations are coupled or not. The results of the Mann-Whitney test are summarized in Table VI In this table, the statistical analysis is presented in two columns, the first one corresponding to the p-value and the second one to the result h of the hypothesis test (a value of h = 1 indicates that the null-hypothesis tested on the 200 trials is rejected). The expected value according to our simulation is also given into brackets. A value of h = 1 proves the presence of an information flow. For the statistical analysis on TE (gBIC), all values are identical to the expected values whatever unidirectional or bidirectional propagation flows. As for  21 → TE (BIC), in the case of unidirectional flow, the null-hypothesis is rejected which indicates that a propagation flow from population 2 t o population 1 is erroneously detected.

C. Real Signals

As discussed in Section II.D, we analyzed two real signals recorded during ictal onset. This phase is characterized by a peculiar pattern of fast activity that originates from the entorhinal/hippocampal region and, secondarily, visually propagates to the perirhinal cortex. Signal 1 was recorded in the medial entorhinal cortex (m-EC) and signal 2 in the perirhinal cortex (PRC). The time duration of the fast onset activity phase was 10 seconds. To detect the relations at different time instants, a 4-s sliding window with a step of 0.125 s was chosen for the estimation of the TE index, leading to a total of 49 TE values. In Fig. 6, we represented the evolution of this index with the orders estimated by gBIC and BIC over time during the analyzed period corresponding to this fast onset activity phase. To estimate the actual influence of signal i on signal j, and to evaluate the statistical significance of the TE measure, we tested this index against the null hypothesis of independent time series. To this end, we needed surrogate data for which we destroyed the causal dependency. We used box plots corresponding to this null hypothesis represented respectively at the bottom-left for the relation from signal 1 to signal 2, and, at the bottom-right for the relation from signal 2 to signal 1 in Fig. 6. Let us indicate that, in Fig. 6, the boxes correspond to TE estimation using gBIC knowing that results were quite comparable using BIC in this case of independence. These boxes have been obtained by randomizing the phase [START_REF] Theiler | Testing for nonlinearity in time series: the method of surrogate data[END_REF] of signal 2 without modifying signal 1, to simulate the independence hypothesis 0 H (between signal 1 and signal 2). Note that, if shuffling the phase destroyed the dependence between the two signals, it also destroyed the nonlinearity in signal 2. Even if this effect was not desired, it had only a small influence on the statistical dispersion in TE estimation. Another approach would consist in shuffling the samples time instants. In this case, to test the 0 H hypothesis, the gBIC procedure and TE estimation should be performed using the pairs of signals (

, ). Both approaches (shuffle phase based approach and shuffle time based approach) leave invariant the second-order properties of the modified signal. In our experiment, we considered the whole duration, i.e. 10 seconds of observation, and we chose 4 windows of 4 s econds each, with an overlap of 2 seconds. For each window, we performed 50 realizations by randomizing the phase so as to get 200 values ( 4 50 × ) of the index under 0 H merged in these box pl ots. From the results in Fig. 6, using gBIC, it is obvious that a strong relation exists from signal 1 to signal 2 in the former part of the interval whereas the relation from signal 2 t o signal 1 is completely non-existent. In the latter part of the sequence, the propagation from the entorhinal region to the perirhinal cortex still holds at a non negligible level, this result being consistent with currently accepted propagation pathways. Compared to BIC, the gBIC criterion allowed for getting larger dynamics on the whole sequence and helped for concluding on the a priori predominant flow direction. On the first part of the sequence, it is clear that the m-EC structure drove the PRC one whatever the criterion. For the second part of the sequence, some estimated orders were the same using either BIC or gBIC, which justifies the superposition of the two curves of TE for the relation from signal 1 to signal 2. Since some orders differed, even if the TE index remained low when estimated in the opposite direction (from signal 2 t o signal 1), its values revealed higher when using the classical criterion instead of the new one.

l

IV. DISCUSSION AND CONCLUSION

In this paper, we improved the estimation of transfer entropy by introducing a greedy strategy to derive a suitable order-estimated measure which is decisive in acquiring coherent and trustworthy results, as shown in bivariate linear AR models displaying either unidirectional or bidirectional flows. The approach including the new order-estimated measure was successfully applied to detect effective connectivity in iEEG signals simulated by physiology-based models. Moreover, the introduction of a second index, namely

→ ij

PTE

, allowed us to quantify the contribution of population i towards population j in the global information flow between the two populations under study. When tested on real signals, the proposed approach appeared robust and helpful in the determination of the organization of the epileptogenic zone. On account of computational efficiency, we just investigated bivariate signals in this study and TE was only employed to identify flow direction using pairs of signals. We are aware that the problem of direct and indirect interactions may arise when considering multivariate signals [START_REF] Yang | Extracting information on flow direction in multivariate time series[END_REF]. We therefore plan to extend this work on TE f or multivariate signals to solve this issue. Annex A When adding different regressors in the model (see equation (A.1)) or modifying the values of the weighted coefficients (equation (A.2)), or increasing the delay from one signal to the other one (equation (A.3)), the new gBIC criterion remained robust and allowed for a better estimation of the transfer entropy (Tables A.1 to A.6) compared to the standard criteria. In all of these cases, the unidirectional flow was detected with accuracy. Annex B In the model described by (23), we extended the range of model orders by changing successively the values of one and only one of the four orders xx q , yx q , xy q and yy q in the range [1; 6]. As discussed in Section III.A, the true model orders were

2 = xx q , 0 = yx q , 3 = xy q
, and 2 = yy q

. For TE calculation, when xx q or yx q is modified, only → YX TE can be theoretically impacted. In the same way, when xy q or yy q is modified, only → XY TE can be influenced. Consequently, we represented in Table B the means and standard deviations (shown in parentheses) of the corresponding TE estimated values at the specific radius to compare them with results displayed in Table II. As can be seen from this table, the unidirectional flow from signal x to signal y was perfectly detected whatever the value of xx q as  → YX TE remained null. On the other hand, the transfer entropy from signal y to signal x was erroneous as soon as the order yx q used for TE estimation differed from the true value, 0 = yx q

, not reported in the table, and leading to  0 → = YX TE

. Now, when analyzing the transfer entropy from signal x to signal y with a false value of xy q , the value of  → XY TE was not as high as expected (all the more since the order was lower than the true one). In the same manner,  → XY TE decreased drastically when the order yy q was larger than the true value. To summarize, these examples illustrate the benefits that can be expected when selecting correctly the orders.

  this absence of information flow, the Y process values ( ) l n y have no influence on the X transition probabilities from the state value ( ) k n

  ): = = = = xx yx xy yyqqqqq . Since we do not impose = kl , we estimate independently the four orders parameters so that we may have (

  The tested signals correspond to ictal epileptiform activities obtained from an animal model (epilepsy induced in guinea pig brain with bicuculline). Simultaneous signals were recorded in the medial entorhinal cortex and in the perirhinal cortex. The recorded EEG signals are plotted in Fig. 5. In the following section, a particular phase of these signals, the FOA (Fast Onset Activity), corresponding to the time interval [ ] 30 s, 40 s , is tested by the measures we developed. III. EXPERIMENTS AND RESULTS Simulations were carried out 200 times on 2048-point signals (corresponding to 4-s length EEG signals in current acquisition systems) both for linear and physiology-based models. For each computation of TE, the orders are estimated by the 4 different measures (AIC, BIC, gAIC, and gBIC) discussed in section II, and the values of r is chosen from

3 =

 3 before, the estimated order values using AIC and BIC are all the same, and are equal to the maximum order value, i.e.

  1(d)). Using AIC and BIC (see Figs. 1(a) and 1(b)), we note some non-zero values for  → YX TE

  ln r (Fig.2(a)). Then, the means and standard deviations (shown in parentheses) of the corresponding TE values at this specific radius * r are given in TableIIas well as those of the PTE index. For sake of clarity, the notation  → the technique noted M (the technique refers here to the model order criterion).

  the greatest values when the direction of information flow is from x to y , and  → YX TE (gBIC) is zero when no information flow exists from y to x since the estimated orders are always identical to the optimal orders. On the opposite,  → YX TE (AIC),  → YX TE (BIC), and  → YX TE (gAIC) still provide some non zero values due to different values of the estimated order.

Fig. 1 .

 1 Fig. 1. Means of TE values in the first linear model. The x-axis is the value of radius r (in ln scale), and the y-axis is the value of TE (bits). (a) The orders are estimated by AIC. (b) The orders are estimated by BIC. (c) The orders are estimated by gAIC. (d) The orders are estimated by gBIC. Solid line: flow direction x→y, dotted line: flow direction y→x.

Fig. 2 .Fig. 3 . 2 Fig. 4 .Fig. 5 .

 23245 Fig. 2. Results on TE with the orders estimated by gBIC between signals x and y in the first linear model. (a) Plot of the mean of C n,r with respect to n vs r (in ln scale), where C and r denote average joint probability and radius in the state space respectively. (b) Plot of TE (bits) vs radius r (in ln scale) (same as Fig. 1(d)). Solid line: flow direction x→y, dotted line: flow direction y→x.

1 Fig. 6 .

 16 Fig. 6. Results on TE with the orders estimated by gBIC and BIC for 2 real iEEG signals. Signal 1corresponds to a signal recorded in m-EC and signal 2 to a signal recorded in PRC.

  al. further Chunfeng Yang, Régine Le Bouquin Jeannès, Member, IEEE, Jean-Jacques Bellanger, and Huazhong Shu, Senior Member, IEEE A New Strategy for Model Order Identification and Its Application to Transfer Entropy for EEG Signals Analysis

  . From Table III, the gBIC order estimator is still the more accurate. This criterion leads to the highest values for  →

	(see Table	TE	XY	YX (gBIC) and  → TE	(gBIC)

  .

	TE As expected, the important values of  12 → unidirectional and bidirectional connectivities) and (gBIC) (both for TE (gBIC) (for reciprocal connectivity) indicate that an  21 → actual information flow exists from one population to the other. TE (gBIC) for unidirectional The non significant value of  21 → connectivity indicates that there is no information flow, whereas the strong value for  21 → TE (BIC) is erroneous.
	PTE Focusing on t he PTE index, the value of  21 → (0.4034) is incorrect for unidirectional connectivity since there (BIC) is no causal influence from population 2 t o population 1. In other words, PTE(BIC) fails in identifying unidirectional connectivity. On the contrary, PTE(gBIC) performs well and allows distinguishing unidirectional and bidirectional connectivities. For one-way propagation,  12 → PTE (gBIC) is 1
	and  21 → PTE	(gBIC) is zero. For two-way propagation, where
	12 = = 21 1500 KK	, the values of  12 → PTE	(gBIC) and
	 21 → PTE

TABLE A

 A

	.1 RESULTS ON THE ORDERS ESTIMATED BY THE 4 DIFFERENT MEASURES IN EQUATION (A.1)
	AIC BIC gAIC gBIC	q xx =2 0/200 0/200 190/200 200/200		q yx =0 0/200 0/200 182/200 200/200	q xy =3 200/200 200/200 200/200 200/200			q yy =2 0/200 0/200 183/200 199/200
	TABLE A.2 RESULTS ON TE AND PTE FOR THE MODEL IN EQUATION (A.1)
	AIC BIC gAIC gBIC	x→y 34.5670 (1.8184) 34.5670 (1.8184) 54.6408 (5.1603) 55.7624 (3.8396)	TE	y→x 5.9612 (0.9828) 5.9612 (0.9828) 2.0719 (5.6763) 0.0000 (0.0000)	x→y 0.8536 (0.0176) 0.8536 (0.0176) 0.9711 (0.0778) 1.0000 (0.0000)	PTE	y→x 0.1464 (0.0176) 0.1464 (0.0176) 0.0289 (0.0778) 0.0000 (0.0000)
	TABLE A.3 RESULTS ON THE ORDERS ESTIMATED BY THE 4 DIFFERENT MEASURES IN EQUATION (A.2)
	AIC BIC gAIC gBIC	q xx =2 0/200 0/200 188/200 200/200		q yx =0 0/200 0/200 183/200 200/200	q xy =3 200/200 200/200 200/200 200/200			q yy =2 0/200 0/200 193/200 200/200
		TABLE A.4 RESULTS ON TE AND PTE FOR THE MODEL IN EQUATION (A.2)
	AIC BIC gAIC gBIC	x→y 62.9057 (2.6545) 62.9057 (2.6545) 84.4181 (6.9006) 85.9249 (4.4645)	TE	y→x 19.4879 (2.1043) 19.4879 (2.1043) 5.6985 (15.6086) 0.0000 (0.0000)	x→y 0.7640 (0.0148) 0.7640 (0.0148) 0.9565 (0.1154) 1.0000 (0.0000)	PTE	y→x 0.2360 (0.0148) 0.2360 (0.0148) 0.0435 (0.1154) 0.0000 (0.0000)

TABLE A .

 A 5 RESULTS ON THE ORDERS ESTIMATED BY THE 4 DIFFERENT MEASURES INEQUATION(A.3)

	AIC BIC gAIC gBIC	q xx =2 0/200 0/200 188/200 200/200		q yx =0 0/200 0/200 183/200 200/200	q xy =5 200/200 200/200 200/200 200/200	q yy =2 0/200 0/200 186/200 200/200
	TABLE A.6 RESULTS ON TE AND PTE FOR THE MODEL IN EQUATION (A.3)
	AIC BIC gAIC gBIC	x→y 34.8387 (1.2853) 34.8387 (1.2853) 76.6072 (8.1165) 78.6284 (4.4936)	TE	y→x 3.8929 (0.8332) 3.8929 (0.8332) 0.6591 (3.3145) 0.0000 (0.0000)	x→y 0.8999 (0.0192) 0.8999 (0.0192) 0.9918 (0.0400) 1.0000 (0.0000)	PTE	y→x 0.1001 (0.0192) 0.1001 (0.0192) 0.0082 (0.0400) 0.0000 (0.0000)

This work was supported by the National Basic Research Program of China under Grant 2011CB707904, by the National Natural Science Foundation of China under Grants 61073138, 61271312, and 60911130370, by the Ministry of Education of China under Grant 20110092110023, the Key Laboratory of Computer Network and Information Integration (Southeast University), Ministry of Education, and by Natural Science Foundation of Jiangsu Province under Grant SBK200910055.

C.


Adrien Coulet
email: adrien.coulet@loria.fr

Florent Domenach
email: domenach.f@unic.ac.cy

Mehdi Kaytoue
email: mkaytoue@liris.cnrs.fr

Amedeo Napoli
email: amedeo.napoli@loria.fr

Using Pattern Structures for Analyzing Ontology-Based Annotations of Biomedical Data

Annotating data with concepts of an ontology is a common practice in the biomedical domain. Resulting annotations (data-concept relationships) are useful for data integration whereas the background ontology can guide the analysis of integrated data. Formal Concept Analysis (FCA) allows to build from a binary context a concept lattice that can be used for data analysis purposes. However annotated biomedical data are not binary and a binarization procedure is required as a preprocessing, coming with classical problems, e.g. a trade-off between expressivity and the large number of induced binary attributes. Interestingly, pattern structures offer a general method for building a concept lattice from any set of objects associated with partially ordered descriptions. In this paper, we show how to instantiate this general framework when the space of descriptions is based on an ontology. We illustrate our approach with the analysis of biomedical annotations and we show its capabilities for knowledge discovery.

Introduction

Annotating data with concepts of an ontology is a common practice in the biomedical domain. Resulting annotations define links between data and concepts that are a support for data exchange, data integration and data analysis tasks [START_REF] Daniel | Biomedical ontologies: a functional perspective[END_REF]. Annotations can be obtained in three main ways. In manual annotation, the specification of links between data and concepts is provided by human domain experts. Automated annotation consists in programs parsing data to provide such links. Semi-automated annotation is a combination of the two previous ways, where programs suggest links between data and concepts that are subsequently validated by domain experts [START_REF] Pesquita | Semantic similarity in biomedical ontologies[END_REF]. Here we are interested in the analysis of annotations of several data resources w.r.t. a "reference ontology". This allows a conjoint use of data from different biomedical domains, e.g., molecular biology and medicine. Consequently, the annotation process plays a major role in translational bioinformatics whose objective is to analyze molecular biomedical data and to discover correlations with clinical knowledge [START_REF] Atul | Viewpoint paper: Translational bioinformatics: Coming of age[END_REF]. In this way, the search for hypotheses about molecular mechanisms underlying translational bioinformatics and the discovery of connexions between molecular data and clinical observations can be performed through the analysis of annotations (as links between biomedical data and ontological concepts).

Formal Concept Analysis (FCA) is a mathematical framework for data analysis and knowledge discovery [START_REF] Ganter | Formal Concept Analysis[END_REF]. As such, it can be used for analyzing annotations of biomedical data provided that some adaptations are made. Firstly, annotations are considered as pairs < document, set of concepts > and thus cannot be directly represented within a binary context. Secondly, domain knowledge in the reference ontology is used for the annotation process and should also be taken into account as well in the analysis process. Thus annotations appear like complex data to be analyzed with FCA, requiring at least binarization.

A first solution is given by scaling, which relies on a transformation of nonbinary data into binary data. Several types of scaling are known in the FCA literature, e.g., nominal, ordinal, interordinal [START_REF] Ganter | Formal Concept Analysis[END_REF]. But it is also known that scaling leads to several problems such as arbitrary transformation of data, data loss and a potential binary attribute flooding, forbidding a comprehensive visualization of the results (see for example experiments and discussion in [START_REF] Kaytoue | Revisiting numerical pattern mining with formal concept analysis[END_REF]).

Another solution is to use pattern structures that allows to directly analyze the complex data [START_REF] Ganter | Pattern Structures and Their Projections[END_REF]. In this setting, objects may have a complex description (i.e., non-binary) but the set of descriptions must be partially ordered within a semi-lattice of descriptions. Descriptions can be of many types, i.e., numerical intervals [START_REF] Kaytoue | Mining gene expression data with pattern structures in formal concept analysis[END_REF], set of attributes [START_REF] Ganter | Pattern Structures and Their Projections[END_REF] or graphs [START_REF] Kuznetsov | Learning closed sets of labeled graphs for chemical applications[END_REF]. Pattern structures allow the application of standard FCA algorithms, e.g., for building the concept lattice, to a partially ordered set (a poset) of descriptions. The partial order on descriptions is defined thanks to a so-called similarity operator (also called meet) and an associated subsumption relation. The formalism of pattern structures was introduced in [START_REF] Ganter | Pattern Structures and Their Projections[END_REF] and gained a lot of interest in the last years due to the need for data mining and knowledge discovery associated with the availability of large volumes of web data (i.e., complex data).

In this paper, we present a first approach to analyze annotations based on a reference ontology using the formalism of pattern structures. The present approach can be seen as a materialization of what is termed as structured attribute sets in [START_REF] Ganter | Pattern Structures and Their Projections[END_REF]. A first requirement for using pattern structures is to define descriptions of objects, then a similarity operation with its associated subsumption relation (thus a partial ordering on descriptions). In the present case, annotations are based on concepts of a reference ontology, which is itself based on two posets, a poset of concepts and a poset of relations (here only the poset of concepts will be taken into account). Accordingly, we propose in this paper an original adaptation of the formalism of pattern structures to annotations based on set of concepts from a reference ontology. Here, descriptions of objects are given by sets of concepts. Then, the ordering of concepts in the reference ontology is used to define an original similarity operator and the associated subsumption relation on descriptions. This is -to the best of our knowledge-the first attempt to analyze data annotations with a pattern structure. Moreover, this shows the potential of pattern structures as an effective formalism for dealing with real-world data and providing substantial results. Actually, the resulting concept lattice can be used for guiding the document annotation process, and especially for completing annotations that are given by an automatic annotation tool. In this case, annotations may be wrong or incomplete. The work of a domain expert for correcting and completing annotations can be very time consuming when large data are considered.

The paper is organized as follows. Section 2 recalls fundamental definitions used in the paper. Section 3 presents our adaptation of pattern structures to ontology-based annotations. It introduces also a concrete example about biomedical data for illustrating the approach. Section 4 details the similarity and subsumption operations on descriptions, while Section 5 provides a discussion about the analysis of annotations of biomedical data using our approach.

2 Background definitions

Formal Concept Analysis

We recall here the standard FCA notations and we refer readers to [START_REF] Ganter | Formal Concept Analysis[END_REF] for details and proofs. A formal context (G, M, I) is defined as a set G of objects, a set M of attributes, and a binary relation I ⊆ G × M . (g, m) ∈ I means that "the object g is related with the attribute m through the relation I". Two derivation operators can be defined on sets of objects and sets of attributes as follows, ∀A ⊆ G, B ⊆ M :

A ′ = {m ∈ M ∶ ∀g ∈ A, (g, m) ∈ I} B ′ = {g ∈ G ∶ ∀m ∈ B, (g, m) ∈ I}
The two operators (⋅) ′ define a Galois connection between the power set of objects P(G) and the power set of attributes P(M ). A pair (A, B), A ⊆ G, B ⊆ M , is a formal concept iff A ′ = B and B ′ = A. A is called the extent and B the intent of the concept. The set of all formal concepts, ordered by inclusion of extents (or dually by inclusion of intents), i.e., (A 1 , B 1 ) ≤ (A 2 , B 2 ) iff A 1 ⊆ A 2 (or dually B 2 ⊆ B 1 ), forms a complete lattice [START_REF]Ordres et classification: Algèbre et combinatoire (tome II)[END_REF], called concept lattice.

Pattern structures

A pattern structure can be understood as a generalization of a formal context to describe complex data [START_REF] Ganter | Pattern Structures and Their Projections[END_REF]: An object has a description lying in a semi-lattice where an "intersection" (or meet) is defined. This intersection allows to characterize the similarity between two descriptions, i.e. what they do have in common.

Formally, let G be a set of objects, let (D, ⊓) be a meet-semi-lattice of object descriptions and let δ ∶ G → D be a mapping associating each object with its description.

A ◻ = ⊓ g∈A δ(g) f or A ⊆ G d ◻ = {g ∈ G d ⊑ δ(g)} f or d ∈ (D, ⊓).
These operators form a Galois connection between the power set of objects P(G) and (D, ⊓). Pattern concepts of (G, (D, ⊓), δ) are pairs of the form

(A, d), A ⊆ G, d ∈ (D, ⊓), such that A ◻ = d and A = d ◻ .
For a pattern concept (A, d), d is the pattern intent and is the common description to all objects in A, the pattern extent. When partially ordered by

(A 1 , d 1 ) ≤ (A 2 , d 2 ) ⇔ A 1 ⊆ A 2 (⇔ d 2 ⊑ d 1 )
, the set of all concepts forms a complete lattice called pattern concept lattice.

The operator (⋅) ◻◻ is a closure operator and pattern intents are closed patterns.

Pattern structure have been applied to numerical intervals [START_REF] Kaytoue | Mining gene expression data with pattern structures in formal concept analysis[END_REF] and to graphs [START_REF] Kuznetsov | Learning closed sets of labeled graphs for chemical applications[END_REF].

EL ontologies

Ontologies that are considered in this work are DL ontologies, i.e. are based on a set of concepts and relations represented with a Description Logic (DL) [START_REF]The Description Logic Handbook: Theory, Implementation, and Applications[END_REF].

The EL DL allows for conjunction (∧) and existential restriction (∃r.c) in definitions of concepts [START_REF] Baader | Pushing the el envelope[END_REF]. This simple DL is sufficient for our purpose, together with transitive roles and general concept inclusion axioms i.e., axioms of the form C ⩽ D where C, D can be either atomic or defined concepts. Moreover, the least common subsumer (lcs) of two concepts in EL always exists and can be computed in polynomial time, provided that their is no cycle in concept definitions, i.e., the definition of a concept c i does not include c i itself [START_REF] Baader | Computing least common subsumers in description logics with existential restrictions[END_REF].

In order to avoid any confusion and to make a clear distinction between the DL formalism and the pattern structure formalism, we use the classical logical notations 6 for the EL DL, thus ∧ for conjunction and ⩽ for subsumption, while we keep ⊓ for the similarity operator and ⊑ for the subsumption relation in pattern structures.

In the following, we consider a reference ontology denoted by O based on the EL DL. O is composed of:

-C(O) denotes a set of concepts, and R(O) denotes a set of binary relations7 , -concepts c i in C(O) are partially ordered thanks to a subsumption relation ⩽, where c 1 ⩽ c 2 means that concept c 1 is a sub-concept of c 2 and that every instance of c 1 is an instance of c 2 , -A is a set of axioms involving concepts and relations.

3 Problem statement

The UMLS Semantic Network and semantic types

The UMLS (Unified Medical Language System) is composed of two main components: a set of ontologies of various biomedical domains (such as SNOMED CT, ICD-10, MeSH) and the UMLS Semantic Network [START_REF] Bodenreider | The unified medical language system (umls): integrating biomedical terminology[END_REF]. For a sake of simplicity, we illustrate our study with annotations of a single data resource, DrugBank8 [START_REF] Knox | Drugbank 3.0: a comprehensive resource for 'omics' research on drugs[END_REF], made with a single ontology of the UMLS, the NCI (National Cancer Institute) Thesaurus [START_REF] Sioutos | Nci thesaurus: A semantic model integrating cancer-related clinical and molecular information[END_REF].

The UMLS Semantic Network is a set of broad subject categories, or semantic types, that is used as a high level categorization of concepts of UMLS ontologies [START_REF] Mccray | An upper level ontology for the biomedical domain[END_REF]. An overview of the 133 semantic types is available at http://www.nlm. nih.gov/research/umls/META3_current_semantic_types.html. The Semantic Network organizes semantic types as a simple hierarchy, i.e., a tree denoted hereafter as ST hierarchy . For example, some semantic types are broader than others such as "Organism" that is more general than "Human" or "Anatomical Structure" more general than "Tissue".

Every concept of a UMLS ontology is mapped with one or more semantic types. In addition, the hierarchy of the Semantic Network ST hierarchy can be used to map a concepts c 1 with the set of semantic types that are ancestors of the semantic types of c 1 . For example, if the concept c 1 has for semantic type "Disease or Syndrome", it can be mapped with "Pathologic Function" and "Biologic Function" too (as the laters are ancestors of the former in ST hierarchy ). Accordingly, we are using the hierarchy ST hierarchy to dispose of the full set of semantic types that can be mapped to each concept. Figure 1 illustrates the mappings of some concepts of the NCI Thesaurus with their semantic types.

In our approach, a selection of semantic types chosen by the analyst will be used as upper level categories for concepts annotating biomedical documents.

Building a pattern structure for biomedical annotations

In this work, we are interested in the discovery of associations between semantic categories (i.e., semantic types) of concepts annotating biomedical documents. This knowledge discovery method should take into account domain knowledge, i.e., the NCI Thesaurus, and semantic types. For example, an expert may be interested in a drug-disease association, e.g., Antibiotic-Inflammation, checking whether the association is frequent and searching for a potential associated molecular mechanism.

For analyzing annotations it may be worth to distinguish concepts thanks to domains of interests (kinds of points of view). For example, a domain expert may group concepts according to their membership to distinct ontologies to separate concepts from an ontology on disease and concepts from an ontology on drugs. Accordingly, we consider in this work that the domain expert defines a "scale" For example, let us consider the document DB01082 (gathering data about Streptomycin) in the DrugBank database. Figure 2 shows this document and an annotation with three concepts of the NCI Thesaurus (here the reference ontology O). Moreover, let us consider the ST scale as ST = {"Disease or Syndrome", "Bacterium", "Molecular Function", "Chemical"}. Then the annotation of DB01082 can be read as:

(DB01082, ⟨{T uberculosis}, {}, {P rotein Synthesis}, {Streptomycin}⟩)

Now we have everything for defining the pattern structure (G, (D, ⊓), δ) for analyzing annotations of biomedical documents:

-G = {g 1 , g 2 , ..., g n } is a set of annotated biomedical documents; -O is the reference ontology, i.e., the NCI Thesaurus, and C(O) is the set of concepts of O; -ST = {ST 1 , ST 2 , ..., ST k } is the set of semantic types of the UMLS Semantic Network that defines the scale ST and the dimensions of the annotation vector;

δ(g i ) = (g i , ⟨ST 1 (g i ), ST 2 (g i ), ..., ST k (g i )⟩)
where ST j (g i ) is the set of concepts of semantic type ST j annoting g i . Table 1 gives an example of this pattern structure. The fourth line of the table shows the annotation of the document DB01082 (about Streptomycin). The different columns are filled with the concepts annotating DB01082 w.r.t. the semantic type provided in the header of each column. Now, it remains to define the similarity operation ⊓ between two descriptions δ(g 1 ) and δ(g 2 ):

δ(g 1 ) = (g 1 , ⟨ST 1 (g 1 ), ST 2 (g 1 ), ..., ST k (g 1 )⟩) δ(g 2 ) = (g 2 , ⟨ST 1 (g 2 ), ST 2 (g 2 ), ..., ST k (g 2 )⟩) δ(g 1 ) ⊓ δ(g 2 ) = ⟨ST 1 (g 1 ) ⊓ ST 1 (g 2 ), ST 1 (g 1 ) ⊓ ST 1 (g 2 ), ..., ST k (g 1 ) ⊓ ST k (g 2 )⟩
where ST 1 (g 1 )⊓ST 1 (g 2 ) is the convex hull in O of all concepts in ST 1 (g 1 ) and ST 2 (g 2 ). The definition of the convex hull is made precise in the next section.

Similarity between descriptions

Given an ontology O, and two concepts c 1 and c 2 , the least common subsumer, denoted by lcs(c 1 , c 2 ), is the most specific concept subsuming both c 1 and c 2 w.r.t. the ontology O. Here O is an EL ontology where no cycle appears in concepts definitions. Thus the lcs of two concepts of O always exists [START_REF] Baader | Computing least common subsumers in description logics with existential restrictions[END_REF]. More generally, the operation can be defined (recursively) for a set of concepts C n = {c 1 , c 2 , ..., c n } as follows: ∀n ∈ N, lcs(C n ) = lcs(lcs(C n-1 ), c n )

For example, the lcs of Streptomycin and Antif ungal Agent is Anti -Inf ective Agent (see Figure 2).

The lcs itself could be used to define a simple similarity operation between two descriptions. One potential application here is to use the concept lattice based on ontological patterns to complete annotations associated with biomedical documents. In this way, taking into account the convex hull of a set of concepts allows to consider available concepts in the ontology which are linked to the initial set of concepts of the annotation. Moreover, if one concept has been missed by the annotation process and is available in the ontology, it will be potentially retrieved through the computation of the convex hull of the initial set of concepts. Now, we define the convex hull of two concepts c 1 and c 2 , denoted by CVX(c 1 , c 2 ), as the set of concepts {x 1 , x 2 , ..., x n } verifying:

-x i ⩽ lcs(c 1 , c 2 ), and -either x i ⩾ c 1 and x i ∧ c 1 ≡ c 1 or x i ⩾ c 2 and x i ∧ c 2 ≡ c 2 -x i ≢ ⊺ For example, CVX(Streptomycin, Antif ungal Agent) = {Anti-Inf ective Agent, Antibiotic, Antif ungal Agent, Streptomycin}.
As for the lcs operation, the convex hull operation can be generalized (recursively) to a set of concepts

C p = {c 1 , c 2 , ..., c p }: ∀p ∈ N, CVX(C p ) = CVX(CVX(C p-1 ), c p )
The meet operation on descriptions applies on two vectors with the same dimensions and returns a vector where the components are filled with the convex hull of the two initial sets of concepts. Formally we have:

δ(g 1 ) = (g 1 , ⟨ST 1 (g 1 ), ST 2 (g 1 ), ..., ST k (g 1 )⟩) δ(g 2 ) = (g 2 , ⟨ST 1 (g 2 ), ST 2 (g 2 ), ..., ST k (g 2 )⟩) δ(g 1 ) ⊓ δ(g 2 ) = ⟨ST 1 (g 1 ) ⊓ ST 1 (g 2 ), ST 2 (g 1 ) ⊓ ST 2 (g 2 ), ..., ST k (g 1 ) ⊓ ST k (g 2 )⟩
where

ST i (g 1 ) ⊓ ST i (g 2 ) = CVX(ST i (g 1 ) ∪ ST i (g 2 )).
The convex hull on the union of two sets of concepts is similar to the convex hull on a set of concepts as defined above.

It can be noticed that the definition of the meet operator on concepts can be likened to the the definition of the meet operator for numerical intervals as the convex hull of two intervals (see for example [START_REF] Kaytoue | Mining gene expression data with pattern structures in formal concept analysis[END_REF]). Moreover, similarly as for intervals we have the following property:

δ(g 1 ) ⊓ δ(g 2 ) = δ(g 1 ) iff δ(g 1 ) ⊑ δ(g 2 )
As an illustration let us consider the two objects "Drug1" and "DB01082" and their descriptions δ(Drug1) and δ(DB01082) given in the Table 1. Their meet is

δ(Drug1) ⊓ δ(DB01082) =
⟨{Bacterial Inf ection, M ycobacterial Inf ection, T uberculosis}, {}, {P rotein Synthesis}, {Anti -Inf ective Agent, Antibiotic, Antif ungal Agent, Streptomycin}⟩.

The meet semi-lattice of pattern elements (indeed of convex hulls) defined by the meet operation is given in Figure 3. This semi-lattice is associated with the context of Table 1 and the order defined by the NCI Thesaurus given in Figure 2.

Dually, it is also possible to define a join operation on descriptions, making (D, ⊓, ⊔) a complete lattice. This operation is not necessary for the definition of pattern structures but exists in our case because of the property of D, the space of descriptions. The join of two descriptions δ(g 1 ) and δ(g 2 ) is defined as follows:

δ(g 1 ) ⊔ δ(g 2 ) = ⟨ST 1 (g 1 ) ⊔ ST 1 (g 2 ), ST 2 (g 1 ) ⊔ ST 2 (g 2 ), ..., ST k (g 1 ) ⊔ ST k (g 2 )⟩
Fig. 3. The meet semi-lattice of convex hulls associated with the context represented in Table 1 and the NCI Thesaurus where

ST i (g 1 ) ⊔ ST i (g 2 ) = CVX(ST i (g 1 )) ∩ CVX(ST i (g 2 )).
Actually, the result of the join operation is the set of common concepts in the two convex hulls of ST i (g 1 ) and ST i (g 1 ).

For example, the join of the descriptions of "Drug1" and "DB01082" is:

δ(Drug1) ⊔ δ(DB01082) = ⟨{T uberculosis}, {}, {P rotein Synthesis}, {}⟩.
The intersection of two convex hulls may be empty as shown in the above example. However, it can be noticed that even if δ(g 1 ) and δ(g 2 ) may have no element in common, they can still have a join as illustrates the following example. Suppose that we have only one dimension and let us consider the reference ontology in Figure 2:

δ(g 1 ) = ⟨{Bacterial Inf ection, T uberculosis}⟩ δ(g 2 ) = ⟨{M ycobacterial Inf ection}⟩.
Actually, the results of the meet and join operations on these two descriptions are:

δ(g 1 )⊓δ(g 2 ) = ⟨{Bacterial Inf ection, M ycobacterial Inf ection, T uberculosis}⟩ and δ(g 1 ) ⊔ δ(g 2 ) = ⟨{M ycobacterial Inf ection}⟩.
In addition, we remark that we do not have δ(g 1 ) ⊓ δ(g 2 ) = δ(g 1 ) as δ(g 1 ) is not a convex hull and thus we do not have either δ(g 1 ) ⊑ δ(g 2 ).

Computing Pattern Structures with CloseByOne

In FCA, an efficient way of computing closed formal concepts that are the basic bricks of concept lattices is the algorithm CloseByOne [START_REF] Kuznetsov | A fast algorithm for computing all intersections of objects in a finite semi-lattice[END_REF]. To adapt Close-ByOne to the general case of pattern structures, one has to replace the original Galois connexion, usually denoted (⋅) ′ , with the derivation operator denoted (⋅) ◻ . Below, we give the basic pseudo-code of the algorithm CloseByOne (Algorithms 1 and 2) for computing ontological patterns. In addition to the new derivation operator, one must replace the intersection operation on descriptions (∩) with the meet operation on patterns (⊓, line 5 of Algorithm 2) that is adapted to the nature of the patterns.

A simple implementation of Algorithms 1 and 2 is proposed at github.com/ coulet/OntologyPatternIcfca/.

Analyzing annotations of biomedical data

We illustrate our approach with the analysis of annotations of DrugBank documents with the ontology named the NCI Thesaurus. These annotations are provided by the NCBO (National Center for Biomedical Ontology) Resource Index presented hereafter.

Alg. 1 CloseByOne.

1: L = ∅ 2: for each g ∈ G 3:

process({g}, g, (g ◻◻ , g ◻ )) 4: L is the concept set.

Alg. 2 process(A, g, (C, D)) with C = A ◻◻ and D = A ◻ and < the lexical order on object names.

if {h h ∈ C A and h < g} = ∅ then 2: L = L ∪ {(C, D)} for each f ∈ {h h ∈ G C and g < h} 4: Z = C ∪ {f } Y = D ⊓ {f ◻ } 6: X = Y ◻ process(Z, f, (X, Y )) 8: end if

A repository of annotations: The NCBO Resource Index

The NCBO Resource Index is a repository of annotations automatically populated by a Natural Language Processing tool [START_REF] Jonquet | Ncbo resource index: Ontologybased search and mining of biomedical resources[END_REF]. This tool parses the textual content of several biomedical databases (e.g., DrugBank, OMIM, ClinicalTrial.gov) searching for occurrences of terms refering to concepts of ontologies. When a concept c i is found in a document g i , an annotation i.e., a pair (g i , c i ), is created and stored. On December 18th, 2012, the NCBO Resource Index was containing annotations for 34 databases with concepts of 280 ontologies of the BioPortal [START_REF] Whetzel | Bioportal: enhanced functionality via new web services from the national center for biomedical ontology to access and use ontologies in software applications[END_REF]. The Resource Index can be queried either by a Web user interface 9 or by a REST Web service 10 . We used the second to build sets of annotations.

DrugBank annotations with the NCI Thesaurus

DrugBank is a publicly available database that contains data about drugs, their indications and their molecular targets. The database is organized into documents, or entries, every document gathering data about one drug. Data in DrugBank are for the main part made of texts is natural language. Figure 2 (b) presents the document of DrugBank that concerns Streptomycin.

As described above, the annotations we used to illustrate our approach are annotations of the DrugBank documents using concepts of the NCI Thesaurus. The NCI Thesaurus is a broad domain ontology and consequently its annotations may concern either clinics and molecular biology data that can be conjointly explored in translational bioinformatics. Moreover, the NCI Thesaurus is an EL ontology that does not contain cyclic concept definition. Thus a lcs always exists Fig. 4. The concept lattice corresponding to the pattern structure based on the context in Table 1 and on the NCI Thesaurus. The top concept has the intent with the larger descriptions and consequently its extent includes all the documents (objects). Traversing the lattice downward, the concepts present more specialized extents (with less objects) and more general intents w.r.t. the subsumption relation on ontological patterns. "Ch", "DoS" and "MF" are respectively abbreviations for the semantic types "Chemical", "Disease or Syndrome" and "Molecular Function"). and its processing is tractable. We used the version 12.04 of the NCI Thesaurus encoded in OWL and available on the NCBO Bioportal 11 .

Interpretation

We propose in Table 1 a context including annotations of five DrugBank documents based on concepts to the NCI Thesaurus. Concepts may have four distinct semantic types (the cardinality of ST is four): "Disease or Syndrome", "Bacterium", "Molecular Function" and "Chemical". The meet-semi-lattice of patterns associated with such annotations is depicted in Figure 3 and the corresponding concept lattice is given in Figure 4. Both sets of formal concepts in the semi-lattice and in the concept lattice have been obtained thanks to the implementation of CloseByOne that we adapted to ontological patterns (see subsection 3.4). Now we propose an analysis of the resulting concept lattice shown in Figure 4. Consider that one of our objectives is to repair and complete the annotations associated with biomedical documents. The top formal concept in the lattice has the "largest extent", i.e., the set of all the objects, and the "smallest intent", actually the largest convex hull for the annotations.

Let us consider the two formal concepts in the upper left part of the concept lattice, the first called c #15 has an extent containing "drug1" and "drug5" and the second called c #5 has an extent containing only "drug5". The "Chemical" semantic type (abbreviated "Ch" in Figure 4) of both concepts is {Anti -Inf ective Agent, Antibiotic, Antif ungal Agent}. The "Disease or Syndrome" dimension ("DoS") in both concepts is {Bacterial Inf ection, M ycobacterial Inf ection, T uberculosis} as in the top concept. However, the "Molecular Function" dimension ("MF") is the same for the top concept and c #15 ; i.e., {P rotein Synthesis}, while it is redefined and empty in c #5 . We propose the following interpretation:

-The value of "Chemical" in both c #15 and c #5 is completed (as a convex hull) and is the correct annotation to be associated to document "drug1" and "drug5" for the "Chemical" dimension. This shows how the final concept lattice based on the ontological pattern structure can effectively complete the original annotation process (especially when this process is automated).

-The same remark applies to the "Disease or Syndrome" dimension, which is also completed (as a convex hull). The concept lattice provides once again the complete annotation for both concepts c #15 and c #5 .

Thus, even on this small and toy example, it is possible to understand and verify the usefulness and potential of the approach: the resulting concept lattice yielded by the ontological pattern structure provides the means for completing the initial annotations in a way that respects the reference ontology.

Finally, we experimented the pattern approach on a larger real-world context. We selected 25 drugs of DrugBank out of 173 drugs returned by the query "antibiotic" and we retain the annotations provided by the NCBO Resource Index associated with 4 distinct semantic types. After 4.4 hours, we obtained 204, 801 closed concepts on a computer with two Intel Core 2 Extreme X7900 CPUs and 4GiB of memory. The resulting concept lattice is rather large and the analysis of formal concepts with a domain expert is in progress. We think that the results of the analysis will be in accordance with the analysis presented just above for the toy example.

Pattern structures provide an original and efficient approach within FCA to analyze complex data such as ontology-based annotations of biomedical documents. In this paper, we propose a framework based on pattern structures for dealing with conceptual annotations which are made of sets of concepts represented within an EL ontology. Then we propose a pattern structure providing a classification of biomedical documents according to their annotations and the semantic types of the concepts within the annotations. The resulting concept lattice can be used for analyzing and completing the original annotations.

This work shows that pattern structures are an efficient means for dealing with real-world and complex data. In the present case, more experiments remain to be done as well as a thorough study of the various pattern structures that can be associated to an annotation process depending on one or several ontologies.

  (G, (D, ⊓), δ) is called a pattern structure. Elements of D are called descriptions or patterns and are ordered by a subsumption relation ⊑ such as ∀c, d ∈ D, c ⊑ d ⇐⇒ c ⊓ d = c. A pattern structure (G, (D, ⊓), δ) gives rise to two derivation operators denoted by (⋅) ◻ :

Fig. 1 .

 1 Fig. 1. Detail of the NCI Thesaurus with associated semantic types from the UMLS. Nodes are concepts of the ontology, arrows represent subsumption relationships (⩽). Doted lines map each concept with its semantic type as defined in the UMLS Semantic Network.

Fig. 2 .

 2 Fig. 2. (a) The left part of the Figure shows the NCI Thesaurus ontology; (b) the right part is an excerpt of the document DB01082 of DrugBank related to the Streptomycin drug. Bold arrows connecting (a) and (b) represent the annotation of DB01082.

Table 1 .

 1 A context where objects are DrugBank documents and attributes are semantic types. Each document is annotated with a set of concepts of the NCI Thesaurus (our reference ontology) having distinct semantic types. The document DB01082 of DrugBank (on the fourth line) is annotated with three concepts, including the concept Tuberculosis of semantic type "Disease or Syndrome".

	❍ ❍ ❍ ❍ G ST	Disease or Syndrome	Bacterium	Molecular Function	Chemical
	Drug1	{Tuberculosis, Bacterial Infection}	{}	{Protein Synthesis}	{Antibiotic, Antifungal Agent}
	Drug2 {Bacterial Infection}	{}	{Protein Synthesis}	{}
	Drug3	{Tuberculosis, Bacterial Infection}	{}	{}	{Anti-Infective Agent}
	DB01082	{Tuberculosis}	{}	{Protein Synthesis}	{Streptomycin}
	Drug5	{Tuberculosis, Bacterial Infection}	{}	{}	{Antibiotic, Antifungal Agent}

But not classical in DL.

To avoid confusion, we will use the terms concept for DL ontologies and formal concepts or pattern concepts for FCA.

Publicly available at http://www.drugbank.ca/

Available at http://bioportal.bioontology.org/resource_index

Documented at http://www.bioontology.org/wiki/index.php/Resource_Index_ REST_Web_Service_User_Guide