Fayrouz Allam

Zaki Nossai

Hesham Gomma

Ibrahim Ibrahim

Mona Abdelsalam

Zaki Nossair

A Recurrent Neural Network Approach for Predicting Glucose Concentration in Type-1 Diabetic Patients

Keywords: glucose concentration prediction, type-1 diabetes, neural networks

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

Introduction

The current trend in research is to automatically monitor the blood glucose level and inject the insulin needed to regulate the patient glucose level. Continuous glucose sensors can be coupled with continuous insulin infusion pumps to create a closed-loop artificial pancreas. The model which is used in MPC is required to accurately predict the glucose level for long prediction horizons to compensate for the delay between the CGM readings and the blood glucose values. Since glucose-insulin interaction is a nonlinear, therefore, modeling using neural network as a nonlinear system will give better results [START_REF] Pe´rez-Gandı´a | Artificial Neural Network Algorithm for Online Glucose Prediction from Continuous Glucose Monitoring[END_REF]. However, if the model is required to predict more than one time step ahead, recurrent neural networks should be used. Neural network techniques have been used in the past for predicting glucose levels using CGM readings as input and also within glucose control systems [START_REF] Mougiakakou | A neural network approach for insulin regime and dose adjustment in type 1 diabetes[END_REF]. Oruklu et al., [START_REF] Oruklu | Adaptive control strategy for regulation of blood glucose levels in patients with type 1 diabetes[END_REF] proposed the development of the empirical models that use frequently sampled glucose data. They used autoregressive moving average (ARMA) to predict the future glucose concentration. Their model is based on the virtual subject's glucose concentration obtained from Hovorka model, and they simulate the CGM data by adding Gaussian noise to the synthetic data. The PH of their model is 30 minutes (6 steps) ahead prediction. Zainuddin et al. [START_REF] Zainuddin | A Neural Network Approach in Predicting the Blood Glucose Level for Diabetic Patients[END_REF] developed four different feed forward wavelet neural networks for four intervals: morning, afternoon, evening and night, in order to predict the blood glucose at the end of each interval. The patients need to fill in information about time of glucose measurements, blood glucose values, insulin doses, food, exercises and stress. Perez et al., [START_REF] Pe´rez-Gandı´a | Artificial Neural Network Algorithm for Online Glucose Prediction from Continuous Glucose Monitoring[END_REF] present a NNM which can predict future glucose concentration from previous 5 glucose CGM samples and the current time stamp, while the output is the glucose concentration at the chosen PH. Three different PHs are used: 15, 30, and 45 min. Their model suffers from time delay between the original and predicted glucose concentration. Gani et al., [START_REF] Gani | Universal Glucose Models for Predicting Subcutaneous Glucose Concentration in Humans[END_REF] developed a glucose prediction autoregressive (AR) model of order 30 to make short term, 30-minutes ahead prediction time, without time lag. Pappada et al., [START_REF] Pappada | Development of a Neural Network for Prediction of Glucose Concentration in Type 1 Diabetes Patients[END_REF] used a dataset of different patients obtained by CGM to construct a neural network using NeuroSolutions® software to predict glucose concentration while time varying predictive window from 50-180 minutes is used. In this paper, we present a new prediction algorithm based on a RNN. It is a fully automated prediction system which doesn't need any data to be entered from patients such as in [START_REF] Zainuddin | A Neural Network Approach in Predicting the Blood Glucose Level for Diabetic Patients[END_REF]; therefore our proposed prediction algorithm is more convenient to be used during sleeping hours.

Subjects and Dataset

The evaluation of the proposed algorithm is performed using glucose measurements from 9 type-1 diabetic patients, the average duration of glucose measurements for each patient is 2 days, 288 samples for each day. The glucose measurements for the 9 patients were obtained using Gaurdian® Real Time CGM system (Medtronic-Minimed) which provides a glucose reading every 5 minutes. Our data set consists of 4916 samples. This data is divided into two different subsets, one subset is for training and the second is for testing and validating the model. The data was smoothed using low pass filter of order 11 before using it in training and testing the neural networks. The use of smoothed version of the CGM data reduces the time lag between the predicted glucose and measured glucose values [START_REF] Gani | Universal Glucose Models for Predicting Subcutaneous Glucose Concentration in Humans[END_REF].

The Proposed Algorithm

Our proposed prediction algorithm uses CGM readings only, it is like a nonlinear AR model. Both inputs and delayed feed backed outputs are glucose concentration values. Fig. 1 shows the architecture of the RNN predictor. The input u i is the glucose reading of sample i (current glucose value), and the output y i is the predicted glucose reading at time i+PH. The input u i is entered to a tapped delay line (TDL) element to hold the previous glucose readings, u i-1 ,...,u i-N , where N is the number of inputs. The initial outputs from this TDL are u 1 to u N [means the glucose readings G 1 , G 2 , G 3 , G 4 ,…,G N ]. The network has also feedback input which is the output y i that enters to another TDL to hold y i-1 ,…,y i-M , where M is the PH divided by 5 minutes sampling period. The initial outputs from this (TDL) are y 1 to y M [means glucose readings G N+1 ,…,G N+M ]. The input u i is updated at each new sample time by the glucose reading from the CGM sensor, this means that the prediction is based on the most recent N readings of the glucose concentration. For example, if the PH is 25 minutes, and the input is u 6 [i.e. G 6 ], then the predicted glucose will be y 6 which is an estimate for G 11 . The available initial offline inputs are the first 10 samples of glucose, the first 5 samples are applied directly to the inputs and the samples from G 6 to G 10 are applied one after the other, therefore the predicted output will be G 11 to G 15 , which means 5 step ahead prediction, (each step represents 5 minutes), therefore it performs 25 minutes prediction.

Evaluation Metrics

The performance of the proposed prediction algorithm is evaluated using several metrics: root mean squared error (RMSE) in mmol/L, FIT, and normalized prediction error (NPE). These metrics have the following expressions.

N G G RMSE i i N i 2 1 ) (     (1) 100 ) 1 ( ^      G G G G FIT i i i       N i i N i i i G G G NPE 1 2 1 2 ^) (
Where Gi is the glucose reading of sample i, N is the total number of samples, i G ^ is the estimated glucose value of sample i,  G is the mean of all glucose readings. In addition to the above three metrics, the Clarke error grid approach (Clarke's EGA) [START_REF] Clarke | The original Clarke error grid analysis (EGA)[END_REF] was also used to assess the performance of the prediction algorithm with a clinically acceptable metric.

Experiments and Results

Many network architectures for the RNN are tested to optimize the predicted output. Two hidden layers and 20 neurons in the first hidden layer were found to give the best results. The prediction performance of the network is evaluated [START_REF] Mougiakakou | A neural network approach for insulin regime and dose adjustment in type 1 diabetes[END_REF] (3) (by calculating FIT) for various numbers of inputs, when the number of inputs is changed; the network architecture that optimizes the predicted output is changed, as shown in table 1. Table 2 shows the evaluation measures calculated for several PH values when using 20-13-1 architecture which gives the best predicted output (from table 1). Table 3 shows the error matrix of EGA which contains the percentage of predicted samples that are located in each zone for different PHs.

Performance Comparisons

The RNN gives more accurate predicted output than the NNM because the relationship between glucose concentration and time is highly dynamic. Many NNM architectures are tested to choose the ones that give the best predicted output, the best architectures were 15-10-1, 15-13-1, 16-13-1, 20-13-1 for PH of 15, 30, 45 and 60 minutes respectively. Fig. 2 shows the difference in prediction accuracy between the NNM and the RNN for various PHs. Table 3 shows Clarke's EGA to clinically evaluate the predicted output from the NNM and RNN. 5.98% for NNM of predicted samples located in zone D, which can lead to incorrect treatments and detections. Table 4 shows the RMSE for results that are obtained from the RNN and the SVR model used in [START_REF] Georga | Predictive Modeling of Glucose Metabolism using Free-living Data of Type 1 Diabetic Patients[END_REF].

Discussion

From all of the previous results it can be seen that, at PH= 30 min., there is no time delay in the predicted output from both the RNN and NNM. At PH=45 min., there is a wider oscillation in the predicted profile of the NNM's output. These oscillations result in a large model prediction error. The NNM failed to predict the output for PH=60 min. For short prediction horizons our proposed RNN algorithm is comparable with others that use recurrent networks such as [START_REF] Huang | A Dynamic Model with Structured Recurrent Neural Network to Predict Glucose-Insulin Regulation of Type 1 Diabetes Mellitus[END_REF] in which they used 6 RNNs to implement the Hovorka diabetic patient physiological model. Our model used real readings of glucose obtained using CGM, and trained one RNN to predict future values of glucose. This difference makes our model more simple and based on real data instead of synthetic data. The model of Huang et al. [START_REF] Huang | A Dynamic Model with Structured Recurrent Neural Network to Predict Glucose-Insulin Regulation of Type 1 Diabetes Mellitus[END_REF] gives FIT=80.9% for PH=30 min., which is less than our result for the same PH. Our results for long PH are still incomparable with that of other algorithms that are trying to predict for longer PHs such as in [START_REF] Huang | A Dynamic Model with Structured Recurrent Neural Network to Predict Glucose-Insulin Regulation of Type 1 Diabetes Mellitus[END_REF] which predicts 2 and 4 hours.

Conclusion

Training the RNN using data obtained from a CGM device, which continuously reads the glucose during the full day for different patients, gives our model some generality.

The RNN needs some fine tuning when it is used with a specific patient. This is one of the differences between our model and some others such as [START_REF] Gani | Universal Glucose Models for Predicting Subcutaneous Glucose Concentration in Humans[END_REF]. Our RNN can accurately predict the glucose values for PH=30 minutes without time delay. At PH=45 minutes the RNN can predict with very little time delay and with acceptable accuracy. But its accuracy is deteriorated for 60 minutes prediction. It is better to use the RNN than the NNM in prediction especially at longer PH, where the accuracy of the NNM is highly deteriorated. Clarke's EGA indicated that the performance of the proposed RNN prediction model is also significant from a clinical point of view. We conclude that, the RNN prediction algorithm succeeds to predict the future glucose values from CGM systems. It can be used for online glucose prediction in model prediction control systems.

Fig. 1 .

 1 Fig. 1. The architecture of the RNN which is used in the prediction algorithm.

Fig. 2 .

 2 Fig. 2. The RNN and NNM prediction results for PH of (a) 15 min., (b) 30 min.,(c) 60 min.

Table 1 .

 1 FIT values when using various number of inputs and PH=40 min. for different network architectures

Table 2 .

 2 Values of RMSE, FIT and NPE when using 40 inputs and (20-13-1) for different PHs

Table 3 .

 3 The Clarke's EGA for Feed Forward Network's Output at Different PHs.

			A	B			C	D			E
	PH NNM	RNN	NNM	RNN	NNM	RNN	NNM	RNN	NNM	RNN
	15	100	100	0	0	0	0	0	0	0	0
	30	98.5	98.6	1.47	1.3	0	0	0	0	0	0
	45	86.7	91.5	10.5	8.4	0	0	2.7	0	0	0
	60	42.09	78.7	50.42	19.3	0	0	5.98	1.95	1.49	0

Table 3

 3 shows that, for PHs less than 60 minutes, all the predicted values are located in clinically accepted zones A and B for RNN. At PH of 45 minutes, 2.7% of predicted outputs from the NNM lay in zone D. At PH of 60 minutes, there are 1.95% for RNN and

	No. of inputs	Network Architecture	FIT (%)
	30	20-7-1		72.16
	30	20-8-1		65.18
	30	20-10-1		74.5
	40	20-10-1		76.72
	40	20-13-1		77.21
	40	18-13-1		68.24
	50	18-13-1		71.5
	50	20-13-1		76.15
	PH	FIT (%)	RMSE (mmol/L)	NPE (%)
	15	95.33	0.14	1.7
	30	85.83	0.42	5.27
	45	72.3	0.84	10.28
	60	56.61	1.32	16.2

Table 4 .

 4 RMSE values for RNN prediction model and SVR based prediction model.


Pari Jahankhani

Juan A Lara

Aurora Pérez

Juan P Valente
email: jpvalente@fi.upm.es

Two different approaches of feature extraction for classifying the EEG Signals

Keywords: Wavelet, Feature extraction, EEG, Classifying 1

The electroencephalograph (EEG) signal is one of the most widely used signals in the biomedicine field due to its rich information about human tasks. This research study describes a new approach based on i) build reference models from a set of time series, based on the analysis of the events that they contain, is suitable for domains where the relevant information is concentrated in specific regions of the time series, known as events. In order to deal with events, each event is characterized by a set of attributes. ii) Discrete wavelet transform to the EEG data in order to extract temporal information in the form of changes in the frequency domain over time-that is they are able to extract non-stationary signals embedded in the noisy background of the human brain. The performance of the model was evaluated in terms of training performance and classification accuracies and the results confirmed that the proposed scheme has potential in classifying the EEG signals.

Introduction

The electroencephalograph (EEG) signal is one of the most widely signal used in the biomedicine field due to its rich information about human tasks.. In practical applications of pattern recognition, there are often diverse features extracted from raw data which needs recognising. Time series modelling has many applications like, for example, feature identification across a group of time series, or model comparison measuring the likeness among groups of time series, or the evolution of one and the same group over time. In actual fact, in many domains, like medicine, the mere observation of the model by the expert can turn out to be very useful in the decision-making process. The relation of EEG signals to the human movements and behaviour has been extensively studied in past decades [START_REF] Zoubir | Seizure detection of newborn EEG using a model approach[END_REF].

A key data mining problem is the construction of feature models from set of time series. In the field of time series data mining, there are well-established methods for comparing two time series, finding subsequence that are repeated several times throughout the same time series and techniques that try to determine whether a time series contains a particular sequence. Also there are techniques that try to generate a representative reference model from a set of time series [START_REF] Chan | Modeling multiple time series for anomaly detection[END_REF], [START_REF] Juan | Discovering similar patterns in time series[END_REF], [START_REF] Chen | A new model for multiple time series based on data mining[END_REF], [START_REF] Papadimitriou | Streaming pattern discovery in multiple timeseries[END_REF]. A powerful method was proposed in the late 1980s to perform time-scale analysis of signals: the wavelet transforms (WT). This method provides a unified framework for different techniques that have been developed for various applications.

Nevertheless, in many cases only particular regions of the series contain relevant knowledge and the data mining techniques should focus on these regions (known as events) [START_REF] Povinelli | Using genetic algorithms to find temporal patterns indicative of time series events[END_REF]. This applies to domains like seismography, the stock market or medicine. In seismography, for example, the only moments of interest are when the time series indicates an earthquake, volcanic activity leading up to the quake, or replications. The lengthy periods between these events provide hardly any information. Figure 1 shows an example of an EEG time series, highlighting an event corresponding to the electrical activity generated by the nervous system in response to a stimulus. 1) Build reference models from a set of time series, based on the analysis of the events that they contain, is suitable for domains where the relevant information is concentrated in specific regions of the time series, known as events. The method enables to define regions of interest according to the knowledge extracted from the domain experts, which is a plus compared with other methods addressing the time series as a whole without taking into account that certain regions can be irrelevant in the domain in question.

2) Discrete Wavelet Transform (DWT) has been applied for the time-frequency analysis of EEG signals and an Adaptive Fuzzy Inference Neural Network System (AFINN) [START_REF] Andrzejak | Indications of nonlinear deterministic and finitedimensional structures in time series of brain electrical activity: dependence on recording region and brain state[END_REF], [START_REF] Jahankhani | Classification Using Adaptive Fuzzy Inference Neural Network[END_REF], [START_REF] Jahankhani | Data Mining an EEG Dataset With an Emphasis on Dimensionality Reduction[END_REF] scheme for the classification using wavelet coefficients.

Method 1: Feature extraction of events

Electroencephalographic devices generate time series that record scalp electrical activity (voltage) generated by brain structures. EEG signals contain a series of waves characterised by their frequency and amplitude. In EEG time series it is possible to find certain types of special waves that are characteristic of some neurological pathologies, like epilepsy. Those waves are known as paroxysmal abnormalities and can be considered as events.

During this research we have taken into account three kinds of events:

 Spike Wave: It is a wave whose amplitude is relatively higher than the rest of waves in the signal. It has a period of between 20 and 70 millisecond.  Sharp Wave: It is a wave whose amplitude is relatively higher than the rest of waves in the signal. It has a period of between 70 and 200 millisecond. Figure 2 shows an example of a sharp wave event.  Spicule: It is a sharp wave with an abrupt change of polarity. The features characterising these events are as follows:

 Duration of the wave.  Amplitude of the wave.

To identify the EEG events and determine their features, the proposed method calculates the point where the polarity of the signal changes as shown in figure 3. The method identifies points where there is a local maximum or minimum whose distance to the polarity change value is higher than a certain threshold (∂). That distance is the amplitude of the event. The duration of the wave is then calculated by analysing the two intersections between the time series and the polarity change value line. Depending on the duration, the event is classified as a spike or a sharp wave, according to the experts" criteria. Finally, those sharp waves that have an abrupt change of polarity are classified as spicules.

Model Generation Method

The model generation method presented here is suited for domains where important information is only confined to certain regions while the remaining of the time series hardly provides any information. In order to deal with events, each event is characterized by a set of attributes.

The model generation method receives a set of time series S = {S 1 , S 2 , ..., S n }, each containing a particular number of events, and generates a reference model M that represents this set of time series. The model M is built on the basis of the most characteristic events. The most characteristic events of S are those events that appear in the highest number of timer series of S.

To find out whether a particular event in a time series S i also appears in another time series S j (j ≠ i), the event has to be characterized with an attribute vector and compared with the other events of the other series. To speed up this process, all the events present in the time series are clustered, so similar events belong to the same cluster. On the one hand, the clustering process is useful to know the different groups of events. On the other hand, it facilitates the extraction of the most characteristic events. Once we have a set of clusters, the objective is to find those clusters containing events that appear in the highest number of time series, that is, characteristic events. Having located those groups with similar events, an exhaustive cluster analysis is run in order to extract the event representative of each of these groups. This will be described later (steps 5 to 9 of the algorithm). These extracted representative events are the characteristic events of S and will be part of the final model.

Let S = {S 1 , S 2 , ..., S n } be a set of n time series and m the typical number of events that appear in the time series of S. The algorithm for generating a reference model M representing the set S is as detailed below (with the purpose of making the algorithm more legible key decisions are justified at the end of the algorithm):

1. Initialize the model. M = Ø.

Identify events.

Extract all the events E v from the series of S and use an attribute vector to characterize each event. This vector covers what the expert considers to be the key features for each type of domain event. This step is domain dependent, as the event characterization will depend on the time series type. To extract the events, the time series is examined in search of regions that meet the conditions identifying each event type defined according to the knowledge extracted from the expert.

Determine the typical number of events m.

m is the typical number of events in each time series of S. At the end of the algorithm it will be discussed how to determine this value.

Cluster events.

Cluster all the events extracted in step 2. Bottom-up hierarchical clustering techniques have been used. Taking into account that the proposal described here should be a general-purpose method and there is no a priori information for specifying the optimum number of clusters in each domain, bottom-up hierarchical clustering is a good option, as it is not necessary to specify the number of clusters k beforehand. We have used Hierarchical Clustering

Repeat steps 5 to 9 m times 5. Get the most significant cluster C k . Determine which cluster C k of all the clusters output in step 4 is the most significant. Cluster significance is measured using Equation (1).

# TS( ) SIGNF(

)

k k C C n  (1)
That is, cluster significance is given by the number of time series that have events in that cluster over the total number of time series n. Events that have already been examined (step 8 and 9) are not taken into account to calculate the numerator.

Extract the event E c that best represents the cluster.

Extract the event that is most representative of the cluster C k , that is, the event E c that minimizes the distance to the other events in the cluster. Let S j be the time series in which the event E c was found.

Add the event E c to the model.

M = M  E c .
8. Mark event E c as examined. 9. Mark the most similar events to E c as examined.

From the cluster C k obtain, for each time series S i ≠ S j , the event E p from S i that is the most similar to the representative event (E c ) output in step 6. Each E p will be represented in the model by the event E c and therefore these E p events will also be discarded in order not to be considered in later iterations.

Return M as a model of the set S.

The most significant clusters, that is, those clusters that contain events present in the highest number of time series were analysed to output the events that are part of the model. To do this, the process of identifying the most significant cluster is repeated m times, outputting a representative and marking as examined both this representative and similar events in each time series. With regard to the algorithm, note that:

a) The identification of events is domain dependent because the criteria to define events in each domain are required. The rest of the algorithm is domain independent and it can be applied to any domain without any change. Figure 4 shows the overall structure of the proposed method that receives a set of time series S and generates a model M that represents it. b) After the representative event of the most significant cluster has been output, it should not be taken into account again for the next iteration, and it is marked as an already examined event. c) A cluster may contain not just one but several events from each time series. For this reason, even if a cluster is selected as the most significant, the cluster in question is not omitted in later iterations. The events already processed are marked as examined and will not be taken into account in future iterations. Another important issue is the number of events making up the model. In this case, we have chosen the mode (m) of the number of events of the time series of S. This decision is based on the fact that if the original time series have a typical number of events m, it makes sense for the model that represents them to have the same number of events m. The typical number of events in the time series of S may not be unimodally distributed. This could happen especially if there are not many time series in the set S. For non-unimodal distributions, we have opted to take the integer value closest to the mean of the number of events.

A last point to be considered is the distance between events that has been used in the algorithm for clustering, representative event selection and discarding similar events. The city block distance is used. Given two vectors, the city block distance calculates the sum of the absolute value of the difference of each of the coordinates of the above vectors:

1 d( , ) p ii i x y x y    (2)
In Equation ( 2), x and y are the vectors (that is, the event descriptors) for comparison and p is the number of coordinates (dimension). Other distance measures have been studied, but the city block distance was finally chosen. The main reason for this choice is that the clustering algorithm uses the mean distance per attribute as the threshold for determining whether or not two elements are similar enough to belong to the same cluster. This mean distance per attribute is obtained simply by dividing the total city block distance d(x,y) by the number of attributes p. The use of the city block distance then saves time as it obviates additional transformations that would make the clustering process more complex to develop and more computationally intensive. Figure 5 shows an example of the application of the proposed method to a set S = {S 1 , S 2 , S 3 , S 4 } of 4 time series (n=4). In this case, S 1 contains 2 events (E 11 and E 12 ), S 2 contains 2 events (E 21 and E 22 ), S 3 contains 3 events (E 31 , E 32 and E 33 ) and finally S 4 contains 2 events (E 41 and E 42 ). Therefore, the typical number of events is 2 (m=2). Once the events are extracted, they are clustered into three different clusters (C 1 , C 2 and C 3 ). Then, the most significant cluster is obtained. To do that, it is necessary to calculate the significance of each cluster according to Equation (EQ). In this case, cluster C 1 have events present in 3 out of the 4 time series, cluster C 2 have events that appear in 1 out of the 4 time series and cluster = 1. Therefore, the most significant cluster is C 3 . In the next step, the event E 12 is extracted as the representative event of the cluster C 3 because E 12 is the event in C 3 that minimizes the distance to the other events in that cluster. Thus, the event E 12 is a characteristic event of S and will be part of the final model M. This process has to be repeated twice (because m=2) to build the final model that consists of the events E 12 and E 32 . The extracted wavelet coefficients provide a compact representation that shows the energy distribution of the EEG signal in time and frequency. In order to reduce the dimensionality of the extracted feature vectors, statistics over the set of wavelet coefficients were used. The following statistical features were used to represent the time-frequency distribution of the EEG signals: Extracted features for the EEG recorded class A and E shown in table 1. For each of these sub-bands, we extracted four measures of dispersion, yielding a total of 20 attributes for sample window. The complete data set consists of two sets denoted by A and E each containing 100 signal-channel EEG segments. The training and test sets of the AFINN classifier were formed by 3200 vectors (1600 vectors from each class). The extracted wavelet coefficients provided a compact representation that shows the energy distribution of the EEG signal in time and frequency. For each of wavelet sub-bands, we have extracted four measures of dispersion, yielding a total of 20 attributes per sample window. Those extracted features for two recorded classes A and E are shown in Table 1.

Table 1: The extracted features of two windows from A & E classes

Evaluation/ Conclusion

A system implementing the described method has been developed. The system has been evaluated by running a battery of experiments using a 10-fold cross validation approach. These experiments were done on time series generated by electroencephalographic devices. During this research, we have used publicly available datasets described in [START_REF] Jahankhani | EEG Signal ClassificationUsing Wavelet Feature Extraction and Neural Networks[END_REF]. The complete data set consists of five sets (denoted A-E) each containing 100 single-channel (100 electrodes) EEG recordings of 5 separate patient classes. For this study, we focused on sets labelled A (healthy patients) and E (epileptic seizure session recordings).

The ultimate aim of the evaluation is to measure how good the model generation method is. For the evaluation of the proposed method, two models were created for each class (M healthy and M epileptic ). The first model (M healthy ) was created from a training set composed of 90 of the 100 healthy patients (set A). The other 10 patients constituted the test set. The second model (M epileptic ) was generated from a training set composed of 90 of the 100 epileptic patients (set E). The other 10 patients were used as test set. The patients in the test set were chosen at random.

Once the models have been created, they have been evaluated by checking whether the M healthy model properly represents the group of healthy patients and whether the M epileptic model is representative of the group of epileptic patients. To do that, we have classified the 20 individuals in the test group according to their similarity to the two created models (that similarity was determined using the time series comparison method proposed in [START_REF] Lara | Comparing posturographic time series through events detection[END_REF]). This process was repeated ten times changing the training set and the test set.

The training data set was used to train the AFINN model for classification of the two classes of EEG signals. The proposed system was trained and tested with the extracted features using discrete wavelet transform of the EEG signals. The simulation results reveal a perfect performance compared to a classic MLP neural network.

The results of the proposed classifier, using 2 different training sets are shown in Table 2.

Fig. 1 .

 1 Fig. 1. Example of an EEG time series.

Fig. 2 .

 2 Fig. 2. Sharp wave event.

Fig. 3 .

 3 Fig. 3. Event taken from an EEG time series.

Fig. 4 .

 4 Fig. 4. Overall structure of the proposed method.

  C 3 have events present in 4 out of the 4 time series of S. Then, the significance of C 1 is SIGNF(C 1 ) = 3 4 = 0.75, the significance of C 2 is SIGNF(C 2 ) = 1 4 = 0.25 and the significance of C 3 is SIGNF(C 3 ) = 4 4

Fig. 5 .

 5 Fig. 5. Example of the application of the proposed method

1 .

 1 Maximum of the wavelet coefficients in each sub-band 2. Minimum of the wavelet coefficients in each sub-band. 3. Mean of the wavelet coefficients in each sub-band. 4. Standard deviation of the wavelet coefficients in each sub-band.

Table 2 :

 2 Comparison of the three methods

	Class	Reference Model	AFINN	MLP
	A	92%	98.12%	94.98%
	E	96%	97.96%	95.86%


Anestis Fachantidis

Ioannis Partalas
email: partalas@csd.auth.gr

Grigorios Tsoumakas

Ioannis Vlahavas
email: vlahavas@csd.auth.gr

Transferring Models in Hybrid Reinforcement Learning Agents

The main objective of transfer learning is to reuse knowledge acquired in a previous learned task, in order to enhance the learning procedure in a new and more complex task. Transfer learning comprises a suitable solution for speeding up the learning procedure in Reinforcement Learning tasks. In this work, we propose a novel method for transferring models to a hybrid reinforcement learning agent. The models of the transition and reward functions of a source task, will be transferred to a relevant but different target task. The learning algorithm of the target task's agent takes a hybrid approach, implementing both model-free and model-based learning, in order to fully exploit the presence of a model. The empirical evaluation, of the proposed approach, demonstrated significant results and performance improvements in the 3D Mountain Car task, by successfully using the models generated from the standard 2D Mountain Car.

Introduction

Reinforcement Learning (RL) is popular for its ability to deal with complex problems using limited feedback. However, RL algorithms often require a considerable amount of training, especially when they are applied to problems with very large and/or continuous state spaces.

Besides this, data efficiency is another important factor. Data efficiency refers to the ability of an agent to compute a good policy with less data. In the case of RL less data means less interaction with the environment which is specially important in real world applications such as robotics, where much interaction with the environment needs time and can be costly. Two of the most important approaches to speed up the learning process and boost data efficiency are Transfer Learning and Model -Based RL algorithms.

Whereas direct RL algorithms compute for long periods before finding good policies Model-based algorithms learn an explicit model of the environment and at the same time, the value function which determines the agent's policy. The model generated from this learning process can either be used for generating simulated experiences or for planning.

Transfer learning refers to the process of using knowledge that has been acquired in a previous learned task, the source task, in order to enhance the learning procedure in a new and more complex task, the target task. The more similar those two tasks are, the easier it is to transfer knowledge between them.

The key insight of our approach is to achieve maximum data efficiency and knowledge exploitation by combining transfer learning with model based learning. We do so by adopting a hybrid Dyna-like approach for our agent, in order to maximize the use of the model. The proposed approach is dubbed Transfer learning in Model-based Reinforcement Learning Agents (TiMRLA). The empirical evaluation, of the proposed approach, demonstrated significant results and performance improvements in the 3D Mountain Car task, by successfully using the models generated from the standard 2D Mountain Car.

Learning the model of a source task

In order to capture the dynamics of the source task environment we have to learn models of its transition and reward functions. Our main concern here is to achieve a balanced presence of the state space in the training sample and maximize the accuracy of the learned models.

To gather experiences and generate the source's task model, we use the model-based RL algorithm Model Based Policy Improvement (MBPI) [START_REF] Kalyanakrishnan | Model-based reinforcement learning in a complex domain[END_REF]. MBPI is a model-based learning algorithm that can learn, on-line, the model of a task from real, gathered, experiences. Thereafter, MBPI uses that model to generate simulated experiences and learns directly from them using Sarsa.

In continuous domains, a value approximation technique is necessary in order to approximate the value of a state using nearby state values. We use the Cerebral Model Articulation Controller (CMAC) [START_REF] Albus | A new approach to manipulator control: The cerebellar model articulation controller (cmac)[END_REF], a well-known tile coding function approximation scheme.

Our specific choice of the MBPI model-based learning algorithm in the source task, should not be considered as a limitation of the method as every other model-based RL algorithm can be equally considered, as long as it can provide us with models of the reward and transition functions of the task. The simulated experiences, in the form of a tuple < s, a, r, s ′ >, form a training set that will be used to learn the transition and reward functions of the task.

In our specific implementation, we use two cascade neural networks to learn the models of these functions. Cascade neural networks were chosen as the characteristics of this type of networks were considered suitable for the RL scheme [START_REF] Rivest | Combining td-learning with cascade-correlation networks[END_REF].

Next, the agent of the source task starts to learn and progressively creates the requested models. Every e episodes MPBI learns a new model of the transition and reward functions, based on all the gathered experiences so far. In our implementation these new models are kept only if their Validation mean square error (MSE) is less than that of the previously generated models. We choose Validation MSE and not Train MSE to avoid overfitting and test the ability of the model to predict from unseen data, such as those in the validation data set.

The TiMRLA algorithm (see Algorithm 3.3) consists of two learning mechanisms which are set up to be compatible. This hybrid-learning approach is followed, as explained later, in order to make maximum use of the source's task model.

Model-based learning

TiMRLA implements model based learning by generating mental (non-real) experiences from a model and directly learning from them with Sarsa (policy improvement) (line 30). This model can either be a source task model thus implementing transfer learning, or the new target task model being trained while learning in the target task. In the beginning of each episode the two models compete for being the current episode's "Model" (line 7). The choice is based on their error on predicting the afterstates from the current agent's state, this error is calculated as the euclidean distance between a predicted state and the one experienced (line 23).

Direct learning

Direct learning is implemented using Sarsa(λ) accompanied with eligibility traces for temporal-credit assignment. This is a model-free approach but in our implementation is accompanied by a model-based mechanism, Action Determinant, that assists the model-free learning process. This approach could be considered closer to that of Reward Shaping [START_REF] Vijaykumar Gullapalli | Shaping as a method for accelerating reinforcement learning[END_REF]. Reward Shaping is a methodology in which directly altering and adjusting the reward an agent receives, can improve its performance by biasing certain state-action values according to some prior knowledge. A simple, in form, quantity called Action Determinant is being calculated in every step (line 17). This quantity is calculated for each action α and consists of two parts. The first, CMAC.V (s ′ model ) is according to our model, the value of the state following our current state if action α was to be applied. The second part of the quantity, r model is the reward predicted from our model for the specific state-action pair (s, α). For each action α ∈ A the sum of the two quantities is summed with the respective Q(s, a) value thus creating |A| new values. These values will be used for action selection with real experiences. So, when being greedy, the agent will make a more informed choice of actions based on the CMAC weights of the predicted afterstates and on the predicted reward.

Transfer Learning

Learning in the source task provided us with the best, in terms of Validation MSE, Neural Networks. These three networks are used by TiMRLA to predict next states, whether they are terminal, and their rewards.

We describe the flow of the proposed algorithm as it executes in the target task. In lines 1-3 some basic initializations are done. The weights in the CMAC, representing our value function Q(s, a), are all set to zero and the initial model to be used is the Source Model. In the beginning of each episode we obtain the start state and we select an action (ǫ-greedy) with respect to the policy. Next, takes place the process of choosing the model that will be used (lines 7-9) in the specific episode. In lines 12-13 we obtain for each step, the current state and its reward. We then optionally calculate the quantities for the Action Determinant by predicting the next state and its value, according to our function approximator as also the reward. This is calculated for each action and passed as a parameter to the function that implements the ǫ-greedy action selection (line 19). Next model-free learning with Sarsa(λ) takes place, with its update rule using an added "CMAC" notation to emphasize the fact that these values come from the approximation of the CMAC.

Next, in lines 21-24 and for each even numbered step (not in every step, for robustness) the two models (target and source) are queried for the new state s ′ (already known) based on the preceding state-action pair (s, a). We then compute and compare two euclidean distances d(s ′ , s ′ target ) and d(s ′ , s ′ source ) where the second argument in the distances, represents the state predictions of the two models. The distances are stored in a tabular form with a window of 1000 distances and the two mean errors are calculated for that window. As mentioned before, at the start of each episode, based on this averaged error (distance), TiMRLA chooses to use the target or the source task model (the one with the less error).

In the next lines the flow continues with the MBPI algorithm for policy improvement using generated experiences from the model.

Experiments

Generalized Mountain Car

The mountain car domain was originally proposed by [START_REF] Moore | Variable resolution dynamic programming: Efficiently learning action maps in multivariate real-valued state-spaces[END_REF] and extended by [START_REF] Satinder | Reinforcement learning with replacing eligibility traces[END_REF]. In the standard 2D task an underpowered car must be driven up to a hill. The state of the environment is described by two continuous variables, the horizontal position x ∈ [-1.2, 0.6], and velocity v x ∈ [-0.007, 0.007].The actions are {Neutral, Left and Right} which modify the velocity by -0.001, 0 and 0.001 respectively. At each time step, an amount of -0.0025 * cos 3x is added to the velocity, in order to take gravity into account.

Each episode starts with the car at the bottom of the hill and ends when x becomes greater than 0.5. At each time step, the agent selects among the three available actions and receives a reward of -1. The objective is to move the car to the goal state, as fast as possible. The main difficulty in this environment is that the force of gravity is greater than the force that the engine is capable to provide. Therefore, the agent controlling the underpowered car must move back and forth, in order to exploit the effect of gravity to its favour and manage to reach the goal state.

The 3D mountain car extends the 2D task by adding an extra spatial dimension. The 3D task was originally proposed in [START_REF] Taylor | Autonomous transfer for reinforcement learning[END_REF]. The state is composed by Algorithm 1 Transferring Models in Reinforcement Learning Agents -TiMRLA four continuous variables, the coordinates in space x, and y ∈ [-1.2, 0.6], as well as the velocities v x and v y ∈ [-0.07.0.07]. The available actions are {Neutral, West, East, South, Nourth}. The Neutral action has no impact on the velocity of the car. The West and East actions add -0.001 and 0.001 to v x respectively, while South and North add -0.001 and 0.001 to v y respectively. Additionally, in order to simulate the effect of gravity, at each time step, -0.0025 * cos 3x and -0.0025 * cos 3y is added to v x and v y respectively. Each episode starts with the car at the bottom of the basin. The goal state is reached when x ≥ 0.5 and y ≥ 0.5. At each time step, the agent receives a reward of -1. In this work we use the mountain car software1 that is based on version 3.0 of the RL-Glue library2 [9].

1: Q ← Q0 D ← 0 2: M odel ←

Setup

After thorough experimentation we concluded to a set of parameters for the MBPI agent in the source task. We set ǫ to 0.5 decaying by a factor of 0.99, the step-size parameter α to 0.2 and γ, the discounting factor, to 1 as this task is an episodic one.Also, depth d of the MBPI algorithm was set to 1 as [START_REF] Kalyanakrishnan | Model-based reinforcement learning in a complex domain[END_REF] show that a greater depth decreases the overall learning performance.

As the standard MBPI algorithms specifies, after every e episodes, training of the networks takes place based on the set of experiences D gathered so far.

We set e to 40 episodes. The neural networks produced are compared with the previously used, based on their MSE, finally keeping the best of them.

As our methodology allows for different state and action variables between the source and the target task, we used manual inter-task mappings. For our experimental domain this is successfully done in [START_REF] Taylor | Transferring instances for model-based reinforcement learning[END_REF]. Manual mapping of the different actions and states should not be seen as a limitation of the method, as there exist methods that are able to automatically learn these mappings [START_REF] Taylor | Autonomous transfer for reinforcement learning[END_REF]. Table 1 presents the inter-task mapping that we used in our experiments.

For the model-based learning part of the algorithm, we got the best results by setting the model-based policy improvement period to 1, which means that policy improvement takes place in every episode. The iterations n were set to 2000, and e and the step-size parameter a, was set to 0 (greedy) and 0.07 respectively. The above values of the parameters are those who gave us the best results.

Finally, for the optional part of our algorithm that specifies the Action Determinant, we run the experiments for both versions of TiMRLA, with and without Action Determinant (TiMRLA+AD).

Both versions of the TiMRLA agent were compared with the standard Sarsa(λ) with a CMAC function approximator. In order to isolate and evaluate the effect of transferring the source model, the two agents had exactly the same parameters for learning from real experiences, with ǫ set to 0.5 , a to 0.2 and γ to 1. 1. Inter-task Mappings used by TiMRLA to construct target task instances from source task data.

3D action 2D action 3D variable 2D variable

Results and Discussion

Figure 2 presents the average cumulative reward that is obtained by each algorithm against the learning episodes.

The graphs were produced from 20 independent trials per method and with a smoothing-averaging window of 20 episodes for clarity purposes. They show a clear performance gain of the proposed algorithm TiMRLA over the standard model-free agent, which is also statistically significant at p < 0.05 (one-tail t test). Its performance demonstrates the positive effect of transferring a model from a source task and also that this model can still assist learning even after many episodes.

Clearly, TiMRLA+AD, shows a more unstable performance having a slight performance gain at the start that decreases gradually after 300 episodes. Its average reward difference with the standard model-free agent is still statistically significant at p<0.05 but for the first 400 episodes. This feature for model based assisting in the action selection, although promising, needs major adjustments such as proper decaying of its influence that will allow the method to positively influence learning in the beginning and gradually decrease its effect.

Furthermore we calculated the mean predictive accuracy of the mountain car 2D model. The mean euclidean distance of a real afterstate s ′ to the predicted one was 0.018 for the state variables x and y.

For the simulated experiences, it is important to note that the best results were obtained with policy improvement taking place in every episode. This shows us that the models succeeded in producing real-like experiences that were positively influencing the learning algorithm just like complementary steps. Although that, the tuned value for the step-size parameter was 0.07 being less than that used for the real experiences which was set to 0.2. This can be seen as a reasonable consequence of using approximate models since they are error prone and can produce at times, errors that could negatively alter the learning procedure.

Finally, an interesting finding is that the best results were obtained with the simulated actions ǫ parameter set to zero. We hypothesize that the model-based learning part, favours greedy action selection because the agent explores anyway by just being placed at random start states. This sub-optimal placing is the exploration mechanism in the agent's mental experiences, and he balances by implementing greedy action selection.

Related Work

This section presents related work in transfer learning and contrasts it with the proposed approach. For a broader review of transfer learning in reinforcement learning domains the reader can refer to a comprehensive survey of the field [START_REF] Taylor | Transfer learning for reinforcement learning domains: A survey[END_REF].

Fernadez and Veloso [START_REF] Fernández | Probabilistic policy reuse in a reinforcement learning agent[END_REF] propose an approach that reuses past policies from solved source tasks to speed up learning in the target task. A restriction of this approach is that the tasks (source and target) must have the same action and state space. Our method allows the state and action spaces to be different between the target and the source task.

Advice based methods have been proposed in [START_REF] Taylor | Cross-domain transfer for reinforcement learning[END_REF][START_REF] Torrey | Skill acquisition via transfer learning and advice taking[END_REF]. Taylor and Stone [START_REF] Taylor | Cross-domain transfer for reinforcement learning[END_REF] proposed a method that uses a list of learned rules from the source task as advice, to guide the learning procedure in the target task. The authors introduce three different utilization schemes of the advice. Furthermore, Torrey et al. [START_REF] Torrey | Skill acquisition via transfer learning and advice taking[END_REF] export rules in first-logic order and translate them into advices.

Taylor et al. [START_REF] Taylor | Transfer learning via inter-task mappings for temporal difference learning[END_REF] proposed a method, named Transfer via inter-task mappings, that initializes the weights of the target task, using the learned weights of the source task, by utilizing mapping functions. The method depends strongly on the function approximation scheme that is used in both tasks, as the function that transfers the weights is altered according to the approximation method.

The transfer of samples have also been proposed in [START_REF] Taylor | Transferring instances for model-based reinforcement learning[END_REF][START_REF] Lazaric | Transfer of samples in batch reinforcement learning[END_REF]. Lazaric et al. [START_REF] Lazaric | Transfer of samples in batch reinforcement learning[END_REF] proposed a method that uses samples gathered in the source task, in order to train a batch reinforcement learning algorithm in the target task. The state and action spaces are not allowed to differ between the two tasks. Taylor et al. [START_REF] Taylor | Transferring instances for model-based reinforcement learning[END_REF] present an algorithm for transferring samples in model-based algorithms. The procedure of transfer is accomplished on-line, that is, when the agent interacts with the environment and the state-action space can change.

Conclusions and Future Work

We examined the viability and benefits of transferring the full model of a task to a different but relevant task. The tasks are allowed to have different actions and state variables by constructing and using the right inter-task mappings. The use of cascade neural networks in an RL scheme showed promising results by providing us with accurate enough models of the transition and reward function of the the source task.

The source's task model was successfully used from a novel algorithm capable of fully exploiting its presence. The proposed algorithm is equipped with two learning mechanisms (hybrid) one for model-free learning and one for modelbased. Empirical results demonstrate a significant performance improvement when compared to the "no-transfer" case. The model of a source task, created with supervised learning, can be efficiently used in a standard model-free learning algorithm either through a model-based add-on module, or directly, with approaches such as that of Reward Shaping. Moreover concurrent learning with the two mechanisms can be compatible as also feasible and efficient.

Proposals for future work include the evaluation of the method with the use of different model-based algorithms and the improving, upon our design, of the model-free learning mechanism with the use of more sophisticated forms of Reward Shaping. Also, the development of new hybrid learning algorithms may further demonstrate the advantages of this approach. Finally, a lot of work has to be done in the field of trustworthy and informative evaluation metrics, in the context of transfer learning for Reinforcement Learning.

Fig. 1 .

 1 Fig. 1. (a) The Mountain Car 2D task. (b) The Mountain Car 3D task.

Fig. 2 .

 2 Fig. 2. Average reward in a period of 500 episodes. The curves are smoothed with a window of 20 episodes

Available at http://library.rl-community.org/

Available at http://glue.rl-community.org/


Hui Wei
email: weihui@fudan.edu.cn

Bo Lang

Qing-Song Zuo
email: qingsong.zuo@gmail.com

A Scale-Changeable Image Analysis Method

Keywords: Neurons, non-classical receptive field, image representation, feedback control circuit

The biological vision system is far more efficient than machine vision system. This is due to the former has rich neural layers for representation and process. In order to obtain a non-task-dependent image representation schema, the early phase of neural vision mechanism is worth simulating. We design a neural model to simulate non-classical receptive field of ganglion cell and its local feedback control circuit, and find it can represent image, beyond pixel level, self-adaptively and regularly. The experimental results prove this method can represent image faithfully with low cost, and can produce a compact and abstract approximation to facilitate successive image segmentation as well as integration operation. This representation schema is good at extracting spatial relationship from different components of image, thus it can be applied to formalize image semantics. Further it can be applied to object recognition or image classification tasks in future.

Introduction

The function of any representation scheme is to capture the essential features of an image and make those features accessible at higher processing layers. The target of image representation is to discover what objects are in an image, and how about their spatial relations. Human vision system generates an internal representation for the external stimulus and requires an internal data structure recording the spatial distribution of stimulus. Such a fundamental data structure is also required by computer vision for recording the image. A better representation should characterize an image by the spatial relationships between them such that a later processing can be adapted to preserver or enhance these relationships.

In order to achieve this goal, a bio-inspired image representation model based on non-classical receptive field (nCRF) is constructed, in which every nCRF denotes a unit of visual information processing, basic, structural and functional. Ganglion cell (GC) merges all stimuli occurring in its RF and reports a composite message upwards for further processing. By means of these dense and regular RFs, many GCs can realize a general representation for any outside stimuli.

Image representation is of primary importance for developing effective image processing and analysis techniques. However, due to lack of a complete understanding of the human vision system and the existence of an infinite number of probable image patterns, it is extremely difficult to construct a universal model for various natural scene images. Facing this dilemma, researchers in this field have resorted to designing separate models for different applications. Some mechanisms proposed for image representation range from color histogram to feature statistics, from spatial frequency based to region based, and from color detection to topology detection. For a more extensive review of image representation techniques, see [START_REF] Eauqueue | Region-based Image Retrieval: Fast Coarse Segmentation and Fine Color Description[END_REF], [START_REF] Deng | An Efficient Color Representation For Image Retrieval[END_REF], [START_REF] Saykol | A Histogram Based Approach For Object Based Query By Shape And Color In Image And Video Database[END_REF], [START_REF] Jeong | Image Retrieval Using Color Histograms Generated By Gauss Mixture Vector Quantization[END_REF]. Another image representation scheme using a set of block pattern models was constructed in which each model a small image block is satisfying certain intensity variation constraints. Although many representation schemes have been proposed, local image patterns was characterized by their statistical properties (i.e. pixels) in most of them. Those operations only differentiating and labeling pixels, this is far away from grasping semantic of an image. As we known, an image is an array of pixels, and this kind of discrete and separate structure cannot manifest semantic organization of image in-depth. Assembled or clustered pixels are more meaningful and effective than single pixel. So a novel representation fashion that can make pixels clustered into some sets or regions is required. Our formulation is unique in that each RF represents explicitly a basic image unit and is described by block-grained. Each block is a compact unit including semantic and every of them should be perfectly divisible units. There is a vivid metaphor describing it, holding a dozen of small tones must be much easier than holding 10-thousands grains of sands with same weight. Integrating pixels of image and use a set to pocket them is an indispensible step in image understanding. It provided a heuristic approach that integrating properly and reducing the load of memory will greatly facilitate semantic emerging.

The arrangement of this paper is as follows. In section II, we formalize the design of the image representation model based on nCRF. Some experiment results are presented in Section III. Finally, in Section IV, we conclude the paper by summarizing the main result and suggesting possible future investigations.

Formulations of The Multi-layer Computational Model Based on nCRF 2.1 A Computational Model for RF Self-adaptive Adjusting Dynamically

Neurophysiological shows that the horizontal cell (HC) with a wide range of electrical synaptic connections, it can receive inputs from large areas of receptor cells (RCs). When a disinhibitory nCRF are solely and fully stimulated by grating patches with low spatial frequency, large area horizontal cells can be simultaneously activated. It affects the activities of RCs by means of spatial summation and feedback and then elicits the (GC) responses to the stimuli within the disinhibitory nCRF through bipolar cells (BCs). Moreover, the amacrine cell (AC) connects many of the nearby GCs in the horizontal direction through its extensive dendritic branches. ACs also interconnect with each other. The breadth of AC's connection exceeds far beyond GC's CRF surround, So AC is properly related to the formation of nCRF with a wide range. Thus, HC and AC play the role of information integration in the outer and inner plexiform layer respectively. GC receive inputs from many neurons in the outer and inner plexiform layer, hence HC and AC are properly connected with the formation of nCRF of retinal GC. Figure 1 shows the details of computational model. According to the computational model show in Fig. 1. The RF's size can self-adaptive adjusting according to the color changing of region is uniform or sharp.

Experiments Results

A series of experiments is performed, from several single attributes to the whole function one-by-one, to test the efficiency of nCRF based representation model.

The Relationship between Entropy and Number of RF

Information entropy is a key concept in information theory. If this concept is applied in image description, it can measure the extent of color change in an image. If we define the quantity of information contained in an image as the number of blocks it contains, then according to the equation of entropy, we get:

  22 1 [ ] ( ) log ( ( )) log ( ) n i i i i I M p m p m E p m         
. where M is a set of classes of blocks with different sizes, M = {m1, m2, ……, mn}, and each class has multiple occurrences, meaning there can be a number of blocks of the same size. A set {p(mi)} defines the numerical likelihood that all classes occur. I[M] represents the information entropy of an image. In terms of the same strategies used in two previous experiments, we calculated information entropies of 300 images, and sorted them into a sequence based on these values.

Fig. 2. RF-based entropy and block-based entropy

Multiscale Image Representation Using nCRF

There exists hierarchical structure among biological tissue of vision, characteristics of external object and the image process. Therefore, the image feature must be described clearly by vision image with multi-scale. Utilizing wavelet transform, a unified mathematic description of the scale of vision image can be achieved with multiscale or multiresolution analysis. The wavelet ψ(x) is stretched as ψ(x/a) (a denote scale) to construct a set of primary functions. When the value of a arises, the primary function can represent the global feature of image, on the contrary, while it descends, it can be used to search for the detailed information of image.

From the pattern of nCRF, the equivalent functional form of single GC or LGN cell for performing calculation in local scope, and the exhaustive distribution of GC and LGN array on the whole vision, it can be concluded that the process of GC and LGN array for image is very similar to the wavelet transform. In fact, DOG function, which is widely used by researcher, is a wavelet function itself, too. A distinctive feature of wavelet transform is the multiscale analysis, which can use different observing scale to detect and represent different size characteristics. The RF of GC and LGN are capable of dynamically changed, it realized the effective extraction and representation for local characteristics of stimulation. For example, the large continuous areas are suitable for being represented by large-sized RF, while the details of image are suitable for being represented by small-sized RF. In this regard, the nervous system shows the adaptability and suitability of the structure and functionality. To explore the mathematical essence of the physiological mechanism, wavelet transform is a very appropriate tool for description. Fig. 3 illustrates the close association between image processing result from the nCRF mechanism and wavelet transformation

The experiment result shows that: (1)With the increasing of level of wavelet transformation, the details of image are gradually faded and the regional characteristic of image emerges simultaneously. The characteristic in these space positions can be exactly rep-resented by the lager-size RF. ( 2) The details of image can be exactly represented by the small-size RF.

(3) Multi-scale analysis come from wavelet transform can be represented by Multisize RF. It can be found that nCRF is an image analysis using neural computation actually, of which function is similar to wavelet transform. This method is credible because its foundation of mathematic is derived from wavelet transform. Although computation ability may not exceed wavelet transform, the representation based on neuron can achieve the synchronous enhance effect in a certain scope by means of synchronous oscillations, and make them prominent from the whole group of neurons. All of these may be the mechanism of neurons, of which can realize the image integration and segmentation.

Discussion

The feasibility for image representation based on nCRF is studied in this work. The real effect of nCRF is newly thought as increasing the encoding efficiency and lowering the redundancy of cells. Visual system can be regarded as a system of information processing and encoding according to information theory [START_REF] Zetzsche | Nonlinear And Extra-Classical Receptive Field Properties And The Statistics Of Natural Scenes[END_REF]. If visual system takes adaptation to natural stimulus as its evolutional goal, it can certainly employ the most effective encoding method for natural images. Since natural images have high relativity, the visual information inputting to RCs is largely redundant [START_REF] Ruderman | Statistics Of Natural Images: Scaling In The Woods[END_REF]. They found that the modulation of nCRF increases the selectivity of the responses of V1 neurons in and the sparsity of the response distribution of groups of neurons. The narrow effective bandwidth of the response curve of single neuron does not decrease the amount of information. Therefore, they thought that the modulation of nCRF help improve the encoding efficiency of visual information. For this purpose, two measures should be taken: first, increase the encoding efficiency of single cell and fully utilized the dynamic characteristics of the cell itself; second, reduce as much as possibly the redundancy in cells possibly and use as few cells as possible to transfer information, which remains the principal problem for our future work.

Fig. 1 .

 1 Fig. 1. A multi-layer neural architecture and how a GC adjusts its nCRF dynamically. (a)is a multiple layers architecture. For the sake of clarity, only one GC and its nCRF is drawn. A nCRF has three rings with positive, negative and positive weight values. (b) This is a neural circuit of nCRF adjusting dynamically, which is a small function unit of (a) in detail. Because a neuron can change its output projection direction according to changing stimulus, nCRF can dynamically change its size. This can be realized by three relay neurons and three switch neurons. A switch neuron imposes its backward control on three relay neurons, and selectively permits only one relay outputting its signal upwards to GC. The relay neuron may have a chance to join one of three different rings of nCRF(c)-(d). With the different switch is turning on ,the same neuron may exclusively participate in forming one of rings of a nCRF. Then a size-changing nCRF come into being. (c) is big one, (d) is middle one and (e) is small one.

Fig. 3 .

 3 Fig. 3. Multiscale image representation using nCRF. The red circle denotes the RF with different size. (a) is original image. (b) is an image after wavelet transform.(c) is the distribution of minimum-size RF, which able to represent the details of face in the image. (d) is the distribution of medium-size RF. (e) is the distribution of large-size RF. (f) is the distribution of full-size RF, which able to represent the region with homogeneous color in the image.

Acknowledgments. This work is supported by 973 Program (Project No. 2010CB327900), the NSFC major project (Project No. 30990263) and Shanghai Key Laboratory of Intelligent Information Processing, China (Grant No. IIPL-09-009).


E G Christodoulou

O D Røe

A Folarin

I Tsamardinos

Information-Preserving Techniques Improve Chemosensitivity Prediction of Tumours Based on Expression Profiles

Keywords: chemosensitivity prediction, variable selection, feature selection, regression, classification

Prior work has shown that the sensitivity of a tumour to a specific drug can be predicted from a molecular signature of gene expressions. This is an important finding for improving drug efficacy and personalizing drug use. In this paper, we present an analysis strategy that, compared to prior work, maintains more information and leads to improved chemosensitivity prediction. Specifically we show (a) that prediction is improved when the GI50 value of a drug is estimated by all available measurements and fitting a sigmoid curve and (b) application of regression techniques often results in more accurate models compared to classification techniques. In addition, we show that (c) modern variable selection techniques, such as MMPC result in better predictive performance than simple univariate filtering. We demonstrate the strategy on 59 tumor cell lines after treatment with 118 fully characterized drugs obtained by the National Cancer Institute (NCI 60 screening) and biologically comment on the identified molecular signatures of the best predicted drugs.

Introduction

Prior work shows that the sensitivity of a tumour to a drug can be predicted better than chance based on the gene-expressions of the tumour [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF], [START_REF] Augustin | Genomic and molecular profiling predicts response to temozolomide in melanoma[END_REF]. This finding paves the way to personalized therapy models. However, the machine learning and statistical analysis employed in prior work processes the data in a way that reduces the available information with potential detrimental effects both on the models' prediction performance as well as the identification of the molecular signatures [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF], [START_REF] Augustin | Genomic and molecular profiling predicts response to temozolomide in melanoma[END_REF], [START_REF] Staunton | Chemosensitivity prediction by transcriptional profiling[END_REF].

First the estimation of the response to a drug in prior work is sub-optimal [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF], [START_REF] Augustin | Genomic and molecular profiling predicts response to temozolomide in melanoma[END_REF]. The response of a tumour depends of course, on the dosage. The National Cancer Institute has treated a panel of 60 cancer cell lines with several thousands drugs and has created a dosage-response profile for each combination of drug and tumour. Often, this profile is summarized with a single value such as the log 10 GI50 where GI50 is the dosage (in µM) of the drug that reduces the natural tumour growth to 50% within 48 hours. NCI, in the majority of cases, estimates log 10 GI50 by piece-wise linear interpolation which are then employed by all prior work (e.g., [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF], [START_REF] Ma | An integrative genomic and proteomic approach to chemosensitivity prediction[END_REF], [START_REF] Staunton | Chemosensitivity prediction by transcriptional profiling[END_REF]). In this paper, we show that estimating the log 10 GI50 values by fitting a sigmoid to the dosage-response profile preserves more information about the effects of the drug that lead to statistically significantly improved predictive performance.

Second, prior work typically quantizes the log 10 GI50 values to create classes of tumours: [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF] and [START_REF] Augustin | Genomic and molecular profiling predicts response to temozolomide in melanoma[END_REF] categorize tumours as sensitive and resistant, while [START_REF] Staunton | Chemosensitivity prediction by transcriptional profiling[END_REF] and [START_REF] Ma | An integrative genomic and proteomic approach to chemosensitivity prediction[END_REF] as sensitive, intermediate, and resistant. This type of quantization allows the application of machine learning classification techniques, variable selection methods for classification tasks, and statistical hypothesis testing techniques for discrete outcomes. Our computational experiments demonstrate that maintaining the exact log 10 GI50 values and employing regression analysis instead of classification is often preferable as it improves chemosensitivity prediction in approximately half of the cases.

Third, prior work often employs simple methods for identifying molecular signatures such as selecting the top k genes that are mostly differentially expressed between different classes of tumours. We show that more sophisticated methods such as the Max Min Parents and Children (MMPC) algorithm for multi-variate feature selection [START_REF] Tsamardinos | The max-min hill-climbing bayesian network structure learning algorithm[END_REF] select more predictive signatures for the same parameter k. We biologically interpret these signatures for the 5 better predicted drugs.

Data and Problem Description

Data Description: Gene-expression profiles were obtained for the NCI-60 cellline panel [6] (these actually contain expressions only for 59 of the 60 cell lines) representing nine types of cancers: 5 Breast, 6 Central Nervous System (CNS), 7 Colon, 6 Leukemia, 10 Melanoma, 8 Lung, 7 Ovarian, 2 Prostate, 7 Renal. The expressions were measured on AffymtrixU133plus2 array containing 54,675 probesets that correspond to about 47,000 transcript variants which in turn represent more than 39,500 of the best characterized human genes. We denote with X i the vector of expressions for cell-line i, X v i the expression value for probeset v on cell-line i, and with X = {X i } the matrix of expressions. The raw data have been subjected to GCRMA normalization before analysis as implemented in the BioConductor platform [7]. The drug-response data for all 59 cell-lines were obtained from the CellMiner database [START_REF] Shankavaram | Cellminer: a relational database and query tool for the nci-60 cancer cell lines[END_REF] for a panel of 118 drugs that are fully characterized. Specifically, for each combination of drug and cell-line, the data contain several pairs of d, r , where d is the log 10 drug dosage and r is the percentage of tissue that survived at 48 hours after treatment. We denote with R i,j the set of such pairs for cell-line i and drug j. Problem Definition: The analysis task we address is to predict the response to a drug of a tissue with expression vector X. The response of cell-line i to a drug j is often characterized with a single number that we denote with GI50 i,j and corresponds to log 10 GI50. GI50 i,j is typically not available in the raw data, thus the value of GI50 i,j is estimated from the data in R i,j . Learning predictive models for GI50 i,j given a vector X is a regression task. Additionally, we are interested in identifying minimal molecular signatures that are optimally predictive of response and that could provide insight into the molecular mechanisms of the drug.

Improving the estimation of GI50

The GI50 i,j values in the publicly available NCI data are usually estimated as follows : The mean response r(d) for each dosage d is calculated and a piecewise linear function is interpolated through these mean values. The estimated GI50 value is the concentration that corresponds to r = 50% on this function, denoted as GI50 P LI i,j . According to the official NCI-60 site [9], this is the methodology followed for estimating the 55% of the GI50 i,j values. The remaining 45% were either approximated (manually, we presume) or chosen to be the highest concentration tested.

We now present an estimation method that employs all available measurements in R i,j . We assume the dosage-response curve to have a sigmoid shape where at 0 dosage (i.e., its logarithm approaches -∞) there is no reduction of the tumour (r=100%) and at infinity the tumour size is reduced to zero (r=-100%).

The equation of a sigmoid that ranges asymptotically between α and α + β and crosses the mid-range at γ is

r = α + β 1 + e (d-γ)δ (1)
where δ is a parameter controlling the slope of the function, r the response and d the dosage (expressed by its logarithm). Considering that asymptotically the drug has no effect at small dosages we set α = -100%; equivalently, at high dosages the tumour is eradicated completely which corresponds to -100% growth, and so we set β = 200. The remaining two parameters γ and δ were estimated using least-squares numerical optimization. Specifically, we used the function nlinfit of Matlab with initial values γ = -5 and δ = 1. This function performs a number of steps towards the steepest descend direction for the parameters γ and δ in order to converge to a good value. In the cases where the procedure would not converge with these initial values, we repeated it 100 times with different initial values for the parameters γ and δ uniformly sampled within [- 15 2] (the range of all concentrations in the data). Out of these 100 repetitions the parameter pair that led to the least mean squared error (MSE) was selected.

The estimated GI50 i,j values are found by setting r=50% and solving Eq. 1.

In certain cases, fitting a sigmoid leads to extreme values. In order to detect the outliers we applied the matlab function deletoutliers. This implements iteratively the Grubbs Test that tests one value at a time [10]. If outliers are found they are trimmed to ±2 * σ j , where σ j is the standard deviation of all currently fitted values to drug j. We denote the final estimates as GI50 Sig i,j . Figure 1(a) shows a graphical depiction of R i,j for the CCRF-CEM (Leukemia) cell-line and Carmustine (BCNU) with the fitted sigmoid superimposed. The corresponding piece-wise linear interpolation segments are also shown in the figure. We now show that this method of estimation leads to improvements in chemosensitivity prediction. The analysis includes the following steps:

Feature Selection: The most commonly applied feature (variable) selection method in the field of personalized medicine is to rank the genes according to their association with the class (equivalently the p-value) and select the top k. We call this method univariate filtering. In our work we additionally employed the Max Min Parents and Children algorithms (MMPC) [START_REF] Tsamardinos | The max-min hill-climbing bayesian network structure learning algorithm[END_REF] to select a minimalsize, optimally predictive set of probe-sets. MMPC is an algorithm that seeks to identify the neighbors of the variable-to-predict in the Bayesian Network capturing the data distribution by taking into account multivariate associations. It has been shown very effective in recent extensive experiments [START_REF] Aliferis | Local causal and markov blanket induction for causal discovery and feature selection for classification part i: Algorithms and empirical evaluation[END_REF] against an array of state-of-the art feature selection methods. In this work, the causal explorer implementation of MMPC was used [START_REF] Statnikov | Causal explorer: A matlab library of algorithms for causal discovery and variable selection for classification[END_REF] with the default values for the parameters.

Regression: We employed SVM Regression to construct the predictive models [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]. In our implementation we used the Radial Basis kernel and all other parameters set to default(we have experimented with other values but they did not lead to better results).

Estimation of Performance:

We used a leave-one-out cross validation protocol due to the small number of samples that we had (59). For each training set, the combination of MMPC and SVM regression produced a predictive model that was applied on the hold-out test sample. Metric of performance: The metric to measure prediction performance is the leave-one-out cross-validated R 2 (coefficient of determination), which is a conservative metric [START_REF] Steel | Principles and Procedures of Statistics[END_REF]. Specifically, for a given drug j, let µ \i be the mean value of GI50 in the data excluding cell line i (training data), GI50 \i the predicted GI50 by the model constructed excluding cell-line i, and GI50 i the GI50 as estimated by the experiments in the corresponding cell line i for drug j. We define:

R 2 j ≡ 1 -i (GI50 i -GI50 \i ) 2 i (GI50 i -µ \i ) 2 (2)
The interpretation of R 2 is that it corresponds to variance (uncertainty) by the models, or the reduction of variance by the use of the geneexpression models.

We have computed R 2 j for all 118 drugs both when the GI50 values are estimated using piece-wise linear interpolation as well as when fitting a sigmoid function, as described above. We denote the corresponding values as R P LI j and R Sig j . The results are shown in Figure 2(a).

The figure shows that GI50 values estimated by the sigmoid are better predicted using the protocol described above. Thus at least for the combination of MMPC and SVM Regression GI50 Sig values facilitate the induction of predictive models vs. using the GI50 P LI . The 95% confidence interval for the median R Sig j -R P LI j as estimated by a Wilcoxon signed-rank test is [0.025, 0.122]. Of course, one could argue that the results may not transfer to other feature selection or regression methods. The results however, corroborate our intuition that the sigmoid estimation better preserves information in the R i,j measurements and given no evidence to the contrary, we would suggest this method of estimation in future analyses.

A second observation on the results is that several drugs are well predicted by the models. More specifically, for the Pearson correlation r between two quantities, Cohen gives the following interpretation guidelines: small effect size, r = 0.1 -0.23; medium, r = 0.24 -0.36; large, r = 0.37 or larger. Interpreting R 2 as r 2 and translating the values we get approximately the intervals [0.01, 0.05), [0.05, 0.13), [0.13, 1]. Under this interpretation for 21 drugs out of the 118, the tumour expression profiles have a large effect in predicting response; for 12 drugs they have a medium effect, and for 19 drugs they have a small effect. On the other hand, 65 out of 118 drugs have a negative size effect, meaning that our prediction does not improve much compared to the prediction by the mean value.

We finally compared the prediction performance of our models using MMPC and univariate filtering as feature selection methodology. We computed the crossvalidated R 2 for both methods on all drugs using the same protocol as before. The k parameter is set to the number of genes returned by MMPC, so that both methods return the signatures of the same sizes. Figure 1(b) presents a histogram of the results. The figure shows that, on average, MMPC returns more informative signatures and thus, it should be preferred.

Comparison of Regression versus Classification

In all related prior work, to the best of our knowledge, classification models have been constructed for predicting GI50 values [START_REF] Potti | Genomic signatures to guide the use of chemotherapeutics[END_REF][START_REF] Staunton | Chemosensitivity prediction by transcriptional profiling[END_REF][START_REF] Ma | An integrative genomic and proteomic approach to chemosensitivity prediction[END_REF]. Given that the latter values are continuous, the authors have quantized them before applying any classifiers, as described in the previous sections. We now show that quantization is sometimes detrimental to performance and regression techniques have greater predictive power.

In this next set of computational experiments we pre-process the GI50 values of each drug to discretize them as described in [START_REF] Ma | An integrative genomic and proteomic approach to chemosensitivity prediction[END_REF]. Specifically, the class C i,j of a cell-line i and drug j is computed as sensitive, intermediate, or resistant if GI50 Sig i,j falls within (-∞, µ j -0.5σ j ], (µ j -0.5σ j , µ j +0.5σ j ], and [µ j -0.5σ j , ∞) respectively, where µ j is the average GI50 value over all cell lines for drug j and σ j the standard deviation.

To evaluate classification, we employed the same overall protocol described in Section 3 with the following modifications: we used multi-class SVM classification instead of SVM Regression [START_REF] Boser | An training algorithm for optimal margin classifiers[END_REF]. SVMs have been very popular and successful classifiers, particularly in bioinformatics [START_REF] Statnikov | A comprehensive evaluation of multicategory classification methods for microarray gene expression cancer diagnosis[END_REF]. We used the libsvm implementation of SVMs with the Radial Basis kernel and all other parameters set to default. In addition, the metric of performance for classification is accuracy, i.e., the percentage of samples whose class is correctly predicted. We denote with A j the leave-one-out cross-validated accuracy of the method on drug j.

Comparing regression vs. classification is not straightforward given that regression outputs a continuous prediction for GI50 Sig i,j while classification outputs its class. To overcome this issue we discretize the output of the regression models to the three stated classes using the same intervals as above. This allows us to compute the cross-validated accuracy of the regression for each drug j, denoted as D j . In other words, A j are computed by first discretizing the data, then using classification, and measuring the accuracy of the output, while D j is computed by using regression, then discretizing the predictions, and computing accuracy. Figure 2(b) shows the histograms of A j , D j , and their difference D j -A j . In some cases regression accuracy scores higher than classification accuracy and in other cases the reverse happens. For example, for drug Vincristine-sulfate (NSC:67574) regression accuracy D is 55,9322 while classification accuracy A is 18,6441. On the other hand, for drug Guanazole (NSC:1895) A is 57,6271 while D is 27,1186.

Biological Interpretation of the Molecular Signatures

The biological correlates of these findings are quite interesting. Below we will discuss the target systems which may connect these genes with tumour resistance and tumour response. Due to space limitation we consider only the top five drugs with the highest R 2 , i.e., the drugs for which we get the larger improvement in prediction when gene expressions are employed. These drugs and their signatures are shown in Table 1.

Table 1. The 5 best predicted drugs (highest R 2 ) along with their gene signatures.

Carmustine (BCNU) is a mustard compound and an alkylating and crosslinking agent, used in oncology since 1972, and still in use for tumours of the central nervous system, Hodgkin, non-Hodgkins disease and myelomatosis [START_REF] Egyhazi | Carmustine-induced toxicity, dna crosslinking and o6-methylguanine-dna methyltransferase activity in two human lung cancer cell lines[END_REF]. FDZ6 (Frizzled homolog 6) belongs to the family of catenin molecules and involved in various processes in developing and adult organisms, but is also implicated in carcinogenesis through the WNT/FZD signaling pathway. IKZF1 or Ikaros zinc finger 1 alteration status accurately predicts relapses [START_REF] Waanders | Integrated use of minimal residual disease classification and ikzf1 alteration status accurately predicts 79% of relapses in pediatric acute lymphoblastic leukemia[END_REF] and a nine-gene signature including deletion of IKZF1 is predictive of chemotherapy response in pediatric acute lymphoblastic leukemia [START_REF] Zuo | A pathway-based gene signature correlates with therapeutic response in adult patients with philadelphia chromosome-positive acute lymphoblastic leukemia[END_REF].

Chlorambucil and Nitrogen mustard hydrochloride compounds are derived from mustard gas and are used in chemical warfare since the first world war. They are both alkylating agents, that interfere with DNA replication and transcription of RNA, and ultimately result in the disruption of nucleic acid function. These compounds are very toxic as well as mutagenic, and thus a more tailored use by genomic signatures would be very important. The indications for Chlorambucil are chronic lymphocytic leukemia, Hodgkin's disease, and other lymphomas. Nitrogen mustard hydrochloride or chlorethamine is active on small cell, non-small cell and prostate cancer cell lines (NCI cell line database), and is indicated for the palliative treatment of Hodgkin's disease (Stages III and IV), other hematological malignancies as well as bronchogenic carcinoma. These two drugs share three genes of interest in their signatures: ARGHAP9 and zinc finger proteins (ZNF107 and ZNF764 respectively) and PTS. ARGHAP9 the RhoGAP domain and is a GTPase activator protein. It is responsible for the transduction of signals from plasma-membrane receptors and for the control of cell adhesion, motility and shape by actin cytoskeleton formation [20]. Zinc finger proteins are involved in DNA recognition, RNA packaging, transcriptional activation, regulation of apoptosis, protein folding and assembly, and lipid binding [START_REF] Laity | Zinc finger proteins: new insights into structural and functional diversity[END_REF]. PTS, encoding 6-pyruvoyltetrahydropterin synthase is crucial in folate metabolism and neutransmitter synthesis [START_REF] Meili | Disease-causing mutations improving the branch site and polypyrimidine tract: pseudoexon activation of line-2 and antisense alu lacking the poly(t)-tail[END_REF]. This has not been associated to cancer previously.

Amsacrine is an aminoacridine derivative and a potent intercalating antineoplastic agent. It is effective in the treatment of acute leukemias and malignant lymphomas, but has poor activity in the treatment of solid tumors. It is frequently used in combination with other antineoplastic agents in chemotherapy protocols. The gene CHMP2B is involved in protein sorting and transport from the endosome to the vacuole/lysosome in eukaryotic cells. It affects the MVB sorting pathway, which plays a critical role in the decision between recycling and degradation of membrane proteins [START_REF] Babst | Escrtiii: an endosome-associated heterooligomeric protein complex required for mvb sorting[END_REF]. Individuals with mutations of this gene develop neurodegenerative disease, probably due to dysfunctional autophagy [START_REF] Rusten | Escrt functions in autophagy and associated disease[END_REF]. RGS19 belongs to the RGS (Regulator of G Protein Signalling) multi-functional, GTPase-accelerating proteins. Recently it was discovered that RGS19 overexpression can enhance cell proliferation [START_REF] Tso | Rgs19 enhances cell proliferation through its c-terminal pdz motif[END_REF]. Single nucleotide polymorphism in FLJ22536 was associated to clinically aggressive neuroblastoma, a childhood malignancy [START_REF] Maris | Chromosome 6p22 locus associated with clinically aggressive neuroblastoma[END_REF]. Interestingly, hypermethylation of the autophagy related gene ATG16L2 was associated with poorer prognosis in terms of molecular response to Imatinib treatment in acute lymphoblastic leukemia, which is also the target disease of amsacrine [START_REF] Dunwell | A genome-wide screen identifies frequently methylated genes in haematological and epithelial cancers[END_REF]. C5orf13 or P311 is involved in invasion and migration of glioblastoma cells [START_REF] Mariani | Identification and validation of p311 as a glioblastoma invasion gene using laser capture microdissection[END_REF].

Dichloroallyl lawson is not an anti-cancer drug yet, and is largely inactive in most cell lines screened. Interestingly, the eIF3a mRNA is found elevated in a number of cancer type and it shown that suppression of eIF3a expression can reverse the malignant phenotype and change the sensitivity of cells to cell cycle modulators [START_REF] Saletta | The translational regulator eif3a: the tricky eif3 subunit![END_REF].

In summary, four of the five signatures predict activity in drugs with a similar mechanism of action, that is DNA damage through alkylation, and mainly active in hematological diseases as lymphomas and central nervous system tumours. Something also worth noticing is that three of our top drugs are nitrogen mustards. The signatures include several genes that are not directly related to the classical DNA damage and repair systems and may infer novel associations and actions of these genes.

Conclusion

Predicting chemosensitivity of tumours from gene expressions is important for selecting treatment, understanding the molecular mechanisms of drug response, and selecting molecular signatures. In this paper, we show that predictive performance can be improved by employing a new method for estimating the GI50 (indication of response to drug), regression algorithms instead of classification, and state-of-the-art, multivariate feature selection. The signatures identified here have several links to cancer progression and resistance to chemotherapy, but the direct relations between the genes and the respective drugs are missing. Knowledge on these relations are still expanding and the methods used to identify those signatures may be important tool for novel biological hypotheses.

Fig. 1 .

 1 Fig. 1. (Left)The drug-response measurements for cell line CCRF-CEM (Leukemia) and Carmustine (BCNU). The Ri,j values are shown, as well as the fitted sigmoid curve (red color) and the respective piece-wise linear interpolation segments (green color), (Right) Histogram of differences of leave-one-out cross validated R 2 when MMPC is used as the feature selection procedure minus the leave-one-out cross validated R 2 when univariate filtering is used as the feature selection procedure.

Fig. 2 .

 2 Fig. 2. (Left) Histogram R 2 for GI50 Sig (GI50's fitted by a sigmoid), GI50 P LI (standard estimation), and their difference, from top to bottom respectively. (Right) Histograms of cross validated classification accuracies Aj , discretized regression accuracies Dj and the Aj -Dj differences, from top to bottom, respectively

Acknowledgements

We would like to thank Matina Fragoyanni for parts of the code, Sofia Triantafyllou, Vincenzo Lagani and Angelos Armen for suggestions and feedback. Finally, we thank FORTH for funding and support of the work.