Nicoletta Nicolaou
email: nicolett@ucy.ac.cy

Julius Georgiou

Autoregressive model order estimation criteria for monitoring awareness during anaesthesia

Keywords: anaesthesia, AR model order estimation, awareness, EEG

This paper investigates the use of autoregressive (AR) model order estimation criteria for monitoring awareness during anaesthesia. The Bayesian Information Criterion (BIC) and the Akaike Information Criterion (AIC) were applied to electroencephalogram (EEG) data from 29 patients, obtained during surgery, to estimate the optimum multivariate AR model order. Maintenance of anaesthesia was achieved with propofol, desflurane or sevoflurane. The optimum orders estimated from the BIC reliably decreased during anaesthetic-induced unconsciousness, as opposed to AIC estimates, and, thus, successfully tracked the loss of awareness. This likely reflects the decrease in the complexity of the brain activity during anaesthesia. In addition, AR order estimates sharply increased for diathermy-contaminated EEG segments. Thus, the BIC could provide a simple and reliable means of identifying awareness during surgery, as well as automatic exclusion of diathermy-contaminated EEG segments.

Introduction

Intraoperative awareness is the phenomenon where patients regain consciousness during surgery. It can lead to severe psychological consequences, such as post-traumatic stress disorder, and economic consequences, such as large insurance compensations. The reported incidence of intraoperative awareness ranges from 0.1-0.8% [START_REF] Bruhn | Depth of anaesthesia monitoring: what's available, what's validated and what's next?[END_REF]; however, due to the amnesic effect of certain anaesthetics, some patients have no recollection of regaining awareness, therefore it is likely that the real incidence of awareness is much higher. The problem becomes more serious due to the routine administration of muscle relaxants with anaesthetic agents, thus placing patients who regain awareness in a state of 'awake paralysis'. Several indices have been applied to monitor the depth of hypnosis (for a review see [START_REF] Rampil | A primer for EEG signal processing in anesthesia[END_REF]) and some devices are commercially available (e.g. Bispectral Index monitor -BIS). These devices convert some combination of various EEG characteristics into a single number from 0-100 representing the level of hypnosis (100 -'fully awake', 0 -'isoelectricity'). Commercial devices suffer from a number of issues, e.g. estimation of the level of hypnosis can be affected by the type of anaesthetic or the administration of a muscle relaxant [START_REF] Dahaba | Different conditions that could result in the bispectral index indicating an incorrect hypnotic state[END_REF][START_REF] Russell | The Narcotrend "depth of anaesthesia" monitor cannot reliably detect consciousness during general anaesthesia: an investigation using the isolated forearm technique[END_REF], and the large inter-subject variability is not captured by the fixed number that represents the level of hypnosis [START_REF] Voss | Monitoring consciousness: the current status of EEG-based depth of anaesthesia monitors[END_REF].

In this work we propose a method for monitoring awareness based on tracking the anaesthetic-induced changes in EEG signal complexity. The proposed method is based on the hypothesis that anaesthetics cause a decrease in the optimum order necessary for EEG modelling with autoregressive (AR) models, which are widely used in EEG analysis (e.g. [START_REF] Anderson | Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks[END_REF][START_REF] Mcfarland | Sensorimotor rhythm-based braincomputer interface (BCI): model order selection for autoregressive spectral analysis[END_REF][START_REF] Franaszczuk | An autoregressive method for the measurement of synchronization of interictal and ictal EEG signals[END_REF]). Since estimating the optimum AR model order is non-trivial several criteria have been introduced for this purpose, with two of the most widespread being the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF] and the Akaike Information Criterion (AIC) [START_REF] Akaike | A new look at the statistical model identification[END_REF]. The ability of the BIC and the AIC to capture this anaesthetic-induced reduction in EEG complexity is investigated. The choice of these measures for this purpose is motivated by the fact that they can provide a computationally simple method of monitoring awareness, while at the same time addressing some of the issues faced by current monitors.

Methods

Dataset

The data used in the particular study were collected from 29 patients (mean age 42.6 ± 20.5; 1 female patient) undergoing routine general surgery at Nicosia General Hospital, Cyprus. The study was approved by the Cyprus National Bioethics Committee and patients gave written informed consent. A detailed description of the dataset and patient exclusion criteria can be found in previous studies (see, for example, [START_REF] Nicolaou | EEG-based automatic classification of awake versus anesthetized state in general anesthesia using granger causality[END_REF]). In summary, anaesthesia was induced with a propofol (1%, 10 mg/ml) bolus. Maintenance of anaesthesia was achieved with (a) intravenous propofol administration (23 patients), or (b) inhalational administration of sevoflurane (4 patients) or desflurane (2 patients). In most patients this was titrated with remifentanil hydrochloride. During induction and/or maintenance some patients received boluses of other drugs, such as neuromuscular blocking agents (cisatracurium, rocuronium, or atracurium) and analgesics, as per surgery requirements.

EEG data were collected using the TruScan32 system (Deymed Diagnostic) at a sampling rate of 256Hz using 19 electrodes according to the 10/20 system (FCz reference). The data were band-pass filtered at 1-100 Hz. Data recording was performed throughout the entire surgical duration (pre-induction, induction, surgery, recovery of consciousness). The point at which the patients stopped responding verbally to commands by the anaesthetist occurred less than a minute after administration of the anaesthetic bolus, depending on patient characteristics. Recovery of consciousness was defined as the point at which the patient responded to verbal commands or tactile stimuli by the anaesthetist. Patient response was expressed either as voluntary muscular movement in response to a command by the anaesthetist or a verbal response. Throughout the recording, time stamps indicating important events, such as anaesthetic induction and recovery of consciousness, were manually inserted in the digital EEG record. These markers were used in subsequent data analysis.

Multivariate AR modelling

An AR model can be considered as a linear filter, which describes a current value of a time series with T samples, X = {x(1), x(2), ..., x(T )}, using weighted information from p past samples:

x(t) = p i=1 a i x(t -i) + e i (t) (1)
where e i is normal noise. The standard AR model assumes a 1-dimensional time series. However, in many situations it is desirable that N variables are modelled together. This is particularly useful when this N -dimensional information is obtained simultaneously, for example through EEG recordings. For this purpose, the multivariate extension of the AR model (MVAR) can be used [START_REF] Penny | Multivariate autoregressive models[END_REF]. In MVAR models information from all N -dimensional observations are taken into account in a unified model of order p:

x(t) = p i=1 A(i)x(t -i) + E(t) (2)
where

x(t) = p i=1    a 11 (i) • • • a 1N (i) . . . . . . . . . a N 1 (i) • • • a N N (i)       x 1 (t -i) . . . x N (t -i)    +    e 1 (t) . . . e N (t)    (3)

AR model order estimation

A number of measures have been developed to estimate the optimum model order. The most widespread are the Bayesian Information Criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF] and the Akaike Information Criterion (AIC) [START_REF] Akaike | A new look at the statistical model identification[END_REF]. These criteria measure the relative goodness of fit of a particular AR model. They provide relative measures of the information loss when the model is applied on the observed data. Given a set of models, the preferred one minimizes the information loss (log-likelihood) estimated from the model order criteria.

The AIC is an information-theoretic measure derived by an approximate minimization of the difference between the true data distribution and the model distribution through the Kullback-Leibler divergence. One way of estimating the AIC is via the residual sum of squares after regression on the data with a model of order p,

AIC(p) = ln|σ 2 e | + 2 k T (4)
where k = N × N × p is the number of parameters estimated for a multivariate model of order p and N dimensions, σ 2 e is the variance of the model residuals, and T is the number of samples. From the definition, the first term of the AIC measures the difference between the log-likelihood of the fitted model and the log-likelihood of the model under consideration. The second term is a penalty term, which penalizes the addition of parameters. The combination of these two terms results in large values of AIC when the model has too few or too many parameters. Thus, the best model is one that provides both a good fit, and that has a minimum number of parameters. This avoids the problem of overfitting. The BIC estimation is very similar to the AIC. However, it differs from the AIC in that it is based on Bayesian formalism and is derived from approximating the Bayes Factor, i.e. the evidence ratios of models. It is defined as:

BIC(p) = ln|σ 2 e | + klnT (5)
Just as the AIC, the BIC also penalizes free parameters, but it does so more strongly. Previous studies have demonstrated that, in general, the BIC is more consistent in contrast to the AIC [START_REF] Zhang | On the convergence of model selection criteria[END_REF]. The AIC performs well in relatively small samples, but is inconsistent and its performance does not improve in large samples; this is in contrast to the BIC, which appears to perform more poorly in small samples, but remains consistent and its performance is improved with larger sample sizes.

Methodology

Dimensionality reduction

The original data space is 19-dimensional (number of electrodes). Fitting a 19-dimensional MVAR model of order p implies estimation of 19 × 19 × p parameters. The requirement for EEG stationarity does not allow the use of long segments; therefore, such a large number of parameters cannot be estimated reliably. In order to reduce the dimensionality, and subsequently the number of estimated parameters, two approaches were followed:

1. EEG aggregates: five brain areas were defined as the aggregate activity of particular electrodes. These areas were: left frontal (LF: electrodes Fp1, F7, F3, T3, C3), right frontal (RF: Fp2, F8, F4, C4, T4), left posterior (LP: T5, P3, O1), right posterior (RP: T6, P4, O2), and midline (Z: Fz, Cz, Pz). In defining these particular aggregates it was taken into account that frontoposterior interactions appear to play an important role in (un)consciousness. Electrode impedance is measured automatically; thus, electrodes with high impedance were subsequently excluded from aggregates estimation. 2. Frontal EEG: the majority of work in the literature, and also commercially available depth of anaesthesia monitors, base their analysis on mono-or bipolar frontal EEG activity. Here, we also applied the BIC and AIC measures on activity from electrode locations Fp1 and Fp2 only, in order to investigate the reliability of such features.

Optimum AR order estimation The EEG data of each patient were analysed in non-overlapping windows of 2s duration. The window size should be large enough to ensure adequate number of samples for accurate parameter estimation, but also short enough to ensure data stationarity. The choice of 2-s fulfilled both these requirements. For each window the optimum order from a range of p = 2, 3, ..., 30 was estimated using the BIC and AIC measures, for both the EEG aggregates and the frontal EEG. For each patient, the differences in the estimated AR orders at the two conditions were tested for statistical significance (ANOVA F-test, α = 0.05).

Results

Figure 1 shows the resulting BIC (a-b) and AIC (c-d) boxplots for wakefulness and anaesthesia for all patients, for both EEG aggregates (a, c) and frontal EEG (b, d). Statistical significance was estimated using ANOVA F-test (α = 0.05). All differences between wakefulness and anaesthesia were significant. Specifically, for aggregate activity, (1) BIC: F = 3535.01, p = 0, and AIC: F = 5452.03, p = 0. For frontal activity, (1) BIC: F = 616.76, p = 0, and AIC: F = 253.00, p = 0. The median AR orders from the aggregate EEG are (a) BIC: 10 (awake) and 7 (anaesthesia); (b) AIC: 30 (awake) and 14 (anaesthesia). For frontal EEG we obtain (a) BIC: 13 (awake) and 10 (anaesthesia); (b) AIC: 30 (awake) and 30 (anaesthesia). Figure 2 shows the median optimum AR orders estimated for each patient using the BIC (a-b) and AIC (c-d) for wakefulness and anaesthesia. All differences were statistically significant (ANOVA F-test, α = 0.05), except for patient S22 (indicated with an '*' on figure 2). A general observation is the overestimation of the optimum AR orders using the AIC measure.

Discussion

The optimum AR orders estimated from the BIC using the EEG aggregate activity track the transitions between wakefulness and anaesthesia well. Administration of anaesthesia causes a decrease in the estimated optimum AR orders, which remain stable throughout the entire duration of surgery, and increase again after the patient regains consciousness at the end of surgery. This effect was observed only for the BIC, with more robust estimates obtained from the aggregate EEG. Despite a similar effect observed for the AIC, this was not robust and identification between awake and anaesthetized states was difficult. Despite the large inter-subject differences in the raw EEG activity, the optimum AR orders estimated for anaesthesia and wakefulness did not display large inter-subject variability.

The more robust BIC performance could be explained through a representative example for a randomly chosen patient and EEG segment (fig. 4). The estimated BIC values are in general higher than the AIC values, which is expected if we take into consideration the second term of eqs. 4 and 5. From fig. 4 it can be seen that the changes in the AIC values become negligible for AR models AIC (c-d) estimated using EEG aggregate activity (a, c) and frontal activity (b, d), for individual patients. All differences between 'awake' and 'anaesthetized' conditions are statistically significant at the 95% significance level, except from one difference which is not statistically significant; this is indicated by an '*'. with orders p > 13. However, the minimum AIC value is obtained at p = 23. In contrast, the BIC decreases until a clear minimum value at p = 10 is reached, after which it increases again. This is a direct result of the BIC estimation, which penalises free parameters more strongly than the AIC. From these observations: (a) the AIC indicates that the goodness of fit is not substantially improved using models of orders p > 13, therefore an appropriate choice of model order without taking into consideration the absolute minimum AIC value would be p = 13; (b) the BIC indicates that the goodness of fit is not substantially improved using models of orders p > 10, as can be seen by the penalty imposed at such models by the BIC. Therefore, the combination of increasing goodness of fit and larger penalty imposed on the number of free parameters ensures that the BIC attains a clear minimum value at the smallest appropriate order.

The proposed method has several advantages over existing methods or commercial monitors: (1) the estimated features remain stable over the entire surgical Fig. 4. An example of BIC and AIC values obtained at a randomly chosen EEG segment and a randomly chosen patient. The BIC displays a distinct minimum, as opposed to the AIC, which tends to decrease asymptotically with the model order, resulting into an overestimation of the model order. For this specific example, the optimum AR orders are 10 and 23 as estimated by the BIC and the AIC respectively. duration; (2) estimations are simple and fast; (3) the features are not affected by the particular anaesthetic protocol; (4) wakefulness/anaesthesia can be identified through simple visual inspection of the estimated features, hence there is no need for conversion to a 'universal' number from 0-100 or the use of a classifier;

(5) the estimated AR orders remain stable after anaesthetic administration is switched off, indicating that the AR orders are not simply a reflection of the metabolic decrease of the anaesthetic agent; and (6) diathermy-contaminated EEG segments can be easily identified. In addition, the EEG data were obtained from actual surgery. This allows us to investigate the stability of the estimated AR order patterns and assess the clinical applicability of the proposed method.

Conclusion

The feasibility of utilizing autoregressive model order criteria for wakefulness / anaesthesia discrimination from EEG activity during surgery was investigated. The Bayesian Information Criterion (BIC) provided a reliable means of tracking the reduction in EEG complexity during anaesthesia through providing lower optimum AR order estimates. This effect can be tracked more reliably when the BIC is applied on aggregate than frontal EEG activity. The same effect was not clearly identifiable with the Akaike Information Criterion (AIC). However, both measures provided AR orders that were statistically different during wakefulness/anaesthesia. The BIC could also be used for automatic identification of diathermy-contaminated EEG segments. Future work will involve (i) investigation of other preprocessing methods for noise removal, e.g. Independent Component Analysis, and (ii) the comparison of the proposed methodology with other methods, including other methods previously proposed by the authors (e.g. [START_REF] Nicolaou | EEG-based automatic classification of awake versus anesthetized state in general anesthesia using granger causality[END_REF]) and other methods from the literature.

Fig. 1 .

 1 Fig. 1. ANOVA Boxplots over all patients for optimum AR orders estimated during awake and anaesthesia using the BIC (a)-(b), and the AIC (c)-(d). Optimum AR orders in (a) and (c) were estimated from the EEG aggregates, while in (b) and (d) from the frontal EEG. All differences between awake and 'anaesthetized conditions are statistically significant at the 95% significance level.

Fig. 2 .

 2 Fig. 2. Median optimum AR orders from BIC (a-b) and AIC (c-d) estimated using EEG aggregate activity(a, c) and frontal activity (b, d), for individual patients. All differences between 'awake' and 'anaesthetized' conditions are statistically significant at the 95% significance level, except from one difference which is not statistically significant; this is indicated by an '*'.

Fig. 3 .

 3 Fig. 3. A: Optimum AR orders for patient S3 throughout entire surgery (maintaince with propofol), estimated by: (a) BIC and (b) AIC from aggregate EEG; and (c) BIC and (d) AIC from frontal EEG. Vertical lines: dashed -induction (AI), dotted -recovery of consciousness (ROC). Sharp peaks in (a) and (c) are caused by diathermy artefacts (see panels B and C for corresponding EEG activity). B: Aggregate EEG for patient S3 throughout surgery (62.3 minutes). Average activity at: (a) left frontal location. High frequency peaks with larger amplitude than the EEG are diathermy artefacts. (b) right frontal location, zoomed in to the amplitude level of the underlying EEG. (c) left posterior at 12-20 min (high diathermy contamination). (d) Right posterior at 12-20 min, zoomed in to EEG amplitude level; and (e) midline at 58-62.3 min (patient recovers consciousness at 56.6 min). C: Corresponding frontal EEG activity at (a) Fp1, and (b) Fp2 (zoomed in).

Acknowledgement

The authors would like to thank Dr. Pandelitsa Alexandrou and Dr. Saverios Hourris (Depr. of Anaesthesia, Nicosia General Hospital, Cyprus), the hospital staff and the anonymous volunteers who participated in this study. This work was co-funded by the Republic of Cyprus and the European Regional Development Fund through the Cyprus Research Promotion Foundation ('DESMI 2008'). Grants: New Infrastructure Project/Strategic/0308/26 and DIDAKTOR/DISEK /0308/20. The authors have no conflict of interest.


Theodoros Koutsandreas
email: th_koutsandreas@hotmail.com

Eleftherios Pilalis
email: epilalis@eie.gr

Aristotelis Chatziioannou

A machine-learning approach for the prediction of enzymatic activity of proteins in metagenomic samples

Keywords: Machine-learning, Enzymes, Proteins, Metagenomics 1

In this work, a machine-learning approach was developed, which performs the prediction of the putative enzymatic function of unknown proteins, based on the PFAM protein domain database and the Enzyme Commission (EC) numbers that describe the enzymatic activities. The classifier was trained with well annotated protein datasets from the Uniprot database, in order to define the characteristic domains of each enzymatic sub-category in the class of Hydrolases. As a conclusion, the machine-learning procedure based on Hmmer3 scores against the PFAM database can accurately predict the enzymatic activity of unknown proteins as a part of metagenomic analysis workflows.

Introduction

The emerging field of Metagenomics comprises the collection and analysis of large amounts of DNA that is contained in an environmental niche [START_REF] Lorenz | Metagenomics and industrial applications[END_REF]. Due to the recent advances in high-throughput sequencing, very large amounts of nucleotide sequences can be generated in short time. Because of the increased volume of data, metagenomics is a promising way to identify novel enzymes and protein functions. However, despite the advances in high-throughput sequencing, the development of appropriate analysis tools remains challenging. Here, we developed a classifier for the prediction of protein enzymatic activity in metagenomic samples. Enzymes are proteins that are used in a wide range of applications and industries, such as Biotechnology and Biomedicine. In order to correlate unknown amino acid/nucleotide sequences with enzyme classes the PFAM database and the Enzyme Nomenclature system were used. PFAM is a database of protein families [START_REF] Finn | The PFAM protein families database[END_REF]. Each family is represented by a multiple sequence alignment which is generated by Hidden Markov Models (HMMs) with the Hmmer3 [START_REF] Finn | HMMER web server: interactive sequence similarity searching[END_REF] program. Proteins consist of one or more functional regions which are called domains, i.e. the existence of a domain in the tertiary structure of a protein, imply a specific function. Thus, proteins of the same family will include identical or similar domains. The PFAM database contains information about protein families, their domains and their architecture. Each entry, represented by a PFAM id, corresponds to a single domain. The similarity of an unknown protein with a protein domain may give great information about its function and its phylogenetic relationships.

The Enzyme Nomenclature (EC) is a numerical classification system for enzymes, based on the chemical reaction that they catalyze. It was developed under the auspices of the International Union of Biochemistry and Molecular Biology during the second half of twenty-first century. Each entry of Enzyme Nomenclature is a four-number code, the enzyme commission number (EC number), which is associated with a specific chemical reaction. Thus each enzyme receives the appropriate EC number according to its chemical activity. The first number of code specifies the major category of catalyzed chemical reaction. There are six major categories of catalysed biochemical reactions: Oxidoreductases: -. The next two numbers specify the subclasses of major class and the last one states the substrate of the reaction. For instance, the EC number 3.1.3.-refers to the hydrolysis of phosphoric mono-ester bond and 3.1.3.11 refers to the hydrolysis of fructose-bisphosphatase which contains a phosphoric mono-ester bond. The classifier developed in the current study was able to classify unknown amino acid sequences originating from metagenomic analysis to hydrolases classes pursuant to the results of Hidden Markov Model detection. The classifier consisted of separate classification models, where the classification type was binominal, i.e. is EC number or is not EC number. In order to train the classifier we used well-annotated proteins and analyzed them with Hmmer3. The result was the score of similarity between an examined sequence and a protein domain. As a result, the features of the training data were the PFAM ids and the vector of each training example included its scores to the appropriate fields.

Dataset

In order to train the classification models, we used well-annotated proteins from the UniProt database [START_REF] Apweiler | UniProt: the Universal Protein knowledgebase[END_REF]. We specifically selected all the reviewed sequences from the reference proteome set of bacteria (taxonomy: 2). The Uniprot database was used because it is a high quality, manually annotated and non-redundant protein database. A total number of 45612 sequences were collected. During the training of the classification models, we selected an amount of known proteins according to their EC numbers. A separate classification model was trained for each EC number, using binominal training data as positive and negative examples. The positive examples were sequences that belonged to a specific EC number (for example 3.1.3.1). In contrast, sequences that belonged to the same EC upper class, but differ in the last digit (i.e. 3.1.1.-), were the negative examples. In this way, we aimed at the detection of differences in features context, which separated a specific EC number protein family from all the other proteins whose EC number differed in the last digit.

Methods

Training of enzymatic classification models

In order to train separate models for each enzymatic category of the Hydrolases class we implemented a procedure which automatically constructed the corresponding training sets (Fig. 1). In the first step, all sequences that are annotated with the specific EC number were selected from the dataset as positive examples. Each EC number has its upper class number (for EC number 3.-.-.-the upper class is all the no 3.-.-.classes, i.e. all enzymes which are not hydrolases). In the second step the procedure selected, as negative example set, an equal amount of sequences that belonged to the upper class but not to the specific EC number. Consequently, the difference between positive and negative training data examples was the last EC number digit. Note that there were some conditions to be tested, especially for the EC numbers with three or four digits, during the execution of this step. The procedure stopped if the amount of positive examples was less than five or if negative examples were not found. Thereafter these two sets of examples were analyzed by Hmmer3 against the HMM profiles of the PFAM-A database in order to collect the domain scores as training features. For each enzymatic category an HMM profile library was thus constructed, which contained all the PFAM domains having a score against the sequences in the corresponding training set. For models corresponding to three-or four-digits EC numbers, a custom HMM profile was automatically constructed from the positive example set and was added to the HMM profile library as an additional feature representing the whole sequence length. The training was performed with the k-nearest neighbor (k-NN) algorithm using the Euclidean distance of 4-nearest neighbors, in a10-fold crossvalidation process (stratified sampling). The choice of the parameter k=4 was made based on the size of the smaller training sets (12-15 examples). Thus, a stable value 4 was given to k representing approximately the one third of the smallest training data.

Application of the trained models to unlabeled sequences

The classification procedure automatically performed Hmmer3 analysis on the unlabeled sequences against the PFAM database and collected their scores as features (Fig. 2). Afterwards, the trained models were applied to the unlabeled feature sets starting from the general classification in Hydrolases class (i.e. EC number 3.-.-.-).

Then the procedure continued to the subclasses with two, three and four EC number digits. The classification procedure was not hierarchical top-down as there was not any filtering method during the descending in enzyme subclasses, like the exclusion of sequences which are not annotated as hydrolases in the first classification task. -.-etc) had some false negatives that were correctly classified during the next steps in more detailed EC number classes with three or four digits. However, not all sequences passed through the classification procedures. As mentioned above, the training set of each EC number contained a particular set of features (HMM profile library). In the unlabeled sequences classifica-tion task, the sequences having no score against this feature set were filtered out as they were considered distant from the upper-level EC number category. The volume of unlabeled sequences to be tested was thus reduced as an execution time optimization of the procedure.

Software tools

The aforementioned procedures were implemented in Python scripts. All datasets and results were stored and queried in a MySQL database. The training and the application of the models was performed using RapidMiner (http://rapid-i.com/).

HMMER

Results

The procedure comprised 163 classifiers. The accuracies and F-scores in 10-fold cross validation, of a representative classifiers part, are listed in Table 1. The vast majority of classifiers had an accuracy of above 95%.

Conclusion

In conclusion, the machine-learning procedure based on Hmmer3 scores against the PFAM database performed well and accurately predicted the enzymatic activity of unknown proteins. Future developments will include the use of other protein motifs databases, like CATH and SCOP and the development of more efficient data mining algorithms. The procedure will also be extended to other enzymatic classes (here we focused on the group of Hydrolases) and will be run-time optimized for its application on very large datasets. Finally, it will be implemented as an independent tool and it will be integrated in more extended metagenomic analysis workflows.

  1.-.-.-, Transferases: 2.-.-.-, Hydrolases: 3.-.-.-, Lyases: 4.-.-.-, Isomerases: 5.-.-.-, Ligases: 6.-.-.

  Filtering was avoided because we observed that the classification procedure in main hydrolase class (3.-.-.-) and its subclasses (3.1.-.-, 3.2.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Procedure of training of the EC number classification models

Table 1 .

 1 Table 2 indicates the mean value and the standard deviation of accuracies for each level of enzyme class. Classifiers with two digits (i.e. 3.1.-.-, 3.2.-.-, 3.4.-.-, 3.5.-.-, 3.6.-.-) had a mean accuracy of 94.56%, those with three digits had a mean accuracy of 97.71% and those with four digits had a mean accuracy of 99.39%. The high performances of the classifiers show that the PFAM domain scores were able to separate the training examples accurately. This state indicates a high vector profile difference between two classes and the lack of significant noise. Classifiers accuracies in 10-fold cross-validation

		Accuracy	F_score	EC	Accuracy	F_score
	EC number	(%)	(%)	number	(%)	(%)
	3.-.-.-	90.25	88.40	3.1.1.45	100.00	100.00
	3.1.-.-	94.26	94.42	3.1.1.61	100.00	100.00
	3.5.-.-	93.44	93.76	3.1.1.85	100.00	100.00
	3.-6.-.-	94.76	94.33	3.1.21.2	100.00	100.00
	3.1.1.-	96.03	95.89	3.1.26.3	100.00	100.00
	3.1.13.-	100.00	100.00	3.1.2.6	100.00	100.00
	3.1.2.-	96.12	96.00	3.1.3.1	95.45	90.91
	3.1.21.-	94.74	94.12	3.1.3.5	99.00	98.88
	3.1.22.-	99.20	99.13	3.1.4.17	100.00	100.00
	3.1.3.-	97.21	97.18	3.2.1.52	100.00	100.00
	3.2.1.-	98.39	98.57	3.2.2.27	100.00	100.00
	3.4.11.-	95.29	95.12	3.4.11.9	100.00	100.00
	3.5.3.-	98.64	98.53	3.4.13.9	100.00	100.00
	3.6.1.-	99.31	99.33	3.4.16.4	93.33	94.12
	3.6.5.-	99.29	99.23	3.4.21.107	100.00	100.00
	3.1.11.2	100.00	100.00	3.5.1.1	100.00	100.00
	3.1.11.5	100.00	100.00	3.5.2.5	100.00	100.00
	3.1.11.6	98.29	98.86	3.5.3.6	100.00	100.00
	3.1.13.1	100.00	100.00	3.5.4.4	97.14	96.30
	3.1.1.1	95.83	93.33	3.6.1.7	100.00	100.00
	3.1.1.29	99.31	99.25	3.6.3.12	100.00	100.00
	3.1.1.3	87.50	82.35	3.6.3.30	100.00	100.00
	3.1.1.31	100.00	100.00	3.6.4.12	99.21	99.43

Table 2 .

 2 Mean and standard deviation in function with the amount of digits in EC numbers

	amount of digits	mean value	standard devia-tion
	2	94.86	0.45
	3	97.71	0.19
	4	99.39	0.09

Acknowledgements

. The presented work in this paper has been funded by the "Cooperation" program 09SYN-11-675 (DAMP), O.P. Competitiveness & Entrepreneurship (EPAN II).


Michaël Peigney
email: michael.peigney@polytechnique.org

On Melan's theorem in temperature-dependent viscoplasticity

In plasticity, Melan's theorem is a well-known result that is both of theoretical and practical importance. That theorem applies to elastic-plastic structures under time-dependent loading histories, and gives a sufficient condition for the plastic dissipation to remain bounded in time. That situation is classically referred to as shakedown. Regarding fatigue, shakedown corresponds to the most favorable case of high-cycle fatigue. The original Melan's theorem rests on the assumption that the material properties remain constant in time, independently on the applied loading. Extending Melan's theorem to time fluctuating elastic moduli is a long standing issue. The main motivation is to extend the range of applications of Melan's theorem to thermomechanical loading histories with large temperature fluctuations: In such case, the variation of the elastic properties with the temperature cannot be neglected. In this contribution, an extension of Melan's theorem to elastic-viscoplastic materials with time-periodic elastic moduli is presented. Such a time-dependence may for instance result from time-periodic temperature variations. An illustrative example is presented and supported by numerical results obtained from incremental analysis.

Introduction

For elastic-perfectly plastic structures under prescribed loading histories, the wellknown Melan's theorem [START_REF] Melan | Theorie statisch unbestimmter Systeme aus ideal-plastischen Baustoff[END_REF][START_REF] Symonds | Shakedown in continuous media[END_REF][START_REF] Koiter | General theorems for elastic-plastic solids[END_REF] gives a sufficient condition for the evolution to become elastic in the large-time limit. That situation is classically referred to as shakedown. Intuitively, shakedown means that the plastic strain tends to a limit as time tends to infinity. The Melan's theorem has the distinctive property of being path-independent, i.e. independent on the initial state of the structure. For a parameterized loading history, Melan's theorem gives bounds on the domain of load pa-rameters for which shakedown occurs. Regarding fatigue design, shakedown corresponds to the most beneficial regime of high-cycle fatigue, as opposed to the regime of low-cycle fatigue which typically occurs if the plastic strain does not converge towards a stabilized value [4]. The shakedown theory has been the object of numerous developments, regarding both extensions of the original theorem to various nonlinear behaviors [START_REF] Pham | On shakedown theory for elastic-plastic materials and extensions[END_REF][START_REF] Nguyen | On shakedown analysis in hardening plasticity[END_REF][START_REF] Ahn | Shakedown of coupled two-dimensional discrete frictional systems[END_REF][START_REF] Peigney | Shakedown theorems and asymptotic behaviour of solids in non-smooth mechanics[END_REF][START_REF] Peigney | On shakedown of shape memory alloys structures[END_REF][START_REF] Weichert | A Historical View on Shakedown Theory[END_REF] and numerical methods for assessing the shakedown limits in the space of load parameters [START_REF] Zarka | A New Approach to Inelastic Analysis of Structures[END_REF][START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement cyclique[END_REF][START_REF] Carvelli | Shakedown analysis of defective pressure vessels by a kinematic approach[END_REF][START_REF] Peigney | Approche par contrôle optimal des structures élastoviscoplastiques sous chargement cyclique[END_REF][START_REF] Peigney | An optimal control approach to the analysis of inelastic structures under cyclic loading[END_REF][START_REF] Simon | Shakedown analysis of engineering structures with limited kinematical hardening[END_REF][START_REF] Spiliopoulos | A direct method to predict cyclic steady states of elastoplastic structures[END_REF].

This chapter is concerned with extensions of the Melan's theorem to situations in which the elastic moduli are fluctuating in time, for instance as a result of imposed temperature variations. Whereas the case of temperature-dependent yield limits is well understood [START_REF] Borino | Consistent shakedown theorems for materials with temperature dependent yield functions[END_REF], the case of temperature-dependent elastic moduli remains a long standing issue and has been the object of several conjectures [START_REF] König | A shakedown theorem for temperature dependent elastic moduli[END_REF][START_REF] Halphen | Evolution des structures élastoplastiques dont les coefficients d'élasticité dépendent de la température[END_REF][START_REF] Hasbroucq | Inelastic responses of a two-bar system with temperature-dependent elastic modulus under cyclic thermomechanical loadings[END_REF][START_REF] Hasbroucq | Analytical study of the asymptotic behavior of a thin plate with temperature-dependent elastic modulus under cyclic thermomechanical loadings[END_REF]. The main difficulty is that the proof used in the original Melan's theorem -as well as in most of its knows extensions -crucially relies on some monotonicity properties that are lost when the elastic moduli are allowed to vary in time. For instance, in the case of constant elastic moduli, the distance between two solutions (as measured by the energy norm) is always decreasing with time [START_REF] Nguyen | On shakedown analysis in hardening plasticity[END_REF], which is no longer true when the elastic moduli vary in time (see. [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] for some example). A shakedown theorem has recently been proposed for elastic-perfectly plastic materials with time-periodic elastic properties [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF]. The statement and proof of that theorem differ significantly from the case of constant material properties. A salient result is that time fluctuations of the elastic moduli need to be not too large for shakedown to occur in a pathindependent fashion.

In this contribution, we aim at extending the result of [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] to elastic-viscoplastic materials with time fluctuating elastic moduli. The proof presented in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] for elastic-perfectly plastic materials makes use of the fact that the stress remains in the elasticity domain, which is not necessarily the case in viscoplasticity. In particular, the initial residual stress can be chosen as arbitrarily large, so that the stress is expected to remain outside (and possibly far away from) the elasticity domain -at least on some time interval. This chapter is organized as follows: In Sect. 2, starting from the local constitutive relations and the equilibrium equations, we derive the differential equation that governs quasistatic evolutions of the residual stress. We comment on the special (and important) case of elastic solutions. Sect. 3 is devoted to the statement and proof of a shakedown theorem for elastic-viscoplastic materials with time-periodic material properties, that is the main result of this chapter. An illustrative example is presented in Sect. 4 and supported by numerical results provided by incremental analysis. Some concluding remarks follow.

Quasistatic evolutions of an elastic-viscoplastic medium

Consider an elastic-viscoplastic body occupying a domain Ω in the reference configuration. Under the assumption of infinitesimal strains, the strain tensor ε is derived from the displacement u by

ε = 1 2 (∇u + ∇ T u).
The total strain ε, stress σ and plastic strain ε p at position x and time t satisfy the constitutive equations

ε(x,t) = L(x,t) : σ (x,t) + ε θ (x,t) + ε p (x,t), (1)
ε p (x,t) = φ (σ (x,t), x,t), (2)
where φ is the dissipation potential, taken in the form

φ (σ , x,t) = α 2 |σ -P C (x,t) σ | 2 . (3)
In [START_REF] Melan | Theorie statisch unbestimmter Systeme aus ideal-plastischen Baustoff[END_REF], L is the (symmetric positive definite) elastic moduli tensor and ε θ is the thermal strain tensor. The double product : in (1) denotes contraction with respect to the last two indexes, i.e. (L : σ ) i j = ∑ k,l L i jkl σ lk .

In (3), C (x,t) is the elasticity domain of the material (assumed to be closed and convex), α > 0 is a viscosity parameter (assumed to be independent on (x,t) for simplicity) and P C (x,t) denotes the projection of C (x,t).

The norm | • | in (3) is defined by |σ | = ∑ i, j σ 2
i j for any symmetric tensor σ . As mentioned in Sect. 1, the space and time dependence of L, ε θ and C may reflect imposed variations of the temperature. For instance, the elastic moduli L of most materials depend on the temperature θ , what can be written as L = L(θ ). For imposed variations θ (x,t) of the temperature, the elastic moduli tensor vary as L(θ (x,t)) and can thus be regarded as a function of space and time.

Assuming quasi-static evolutions, the stress field σ satisfies the equilibrium equations div σ

+ f = 0 in Ω , σ • n = T on ∂ Ω T , (4)
where f (x,t) are body forces imposed in the domain Ω and T (x,t) are tractions prescribed on a part ∂ Ω T of the boundary ∂ Ω . Prescribed displacements v(x,t) are imposed on ∂ Ω u = ∂ Ω -∂ Ω T .

Evolution equation for the residual stress

We now use Eqs. (1-4) to derive the equation governing the evolution of the stress field. The space E of stress fields is chosen as a subspace of symmetric second-order tensor fields with square-integrable components, which is known to be a Hilbert space for the scalar product

σ , σ = Ω σ (x) : σ (x)dω.
The associated norm is denoted by , i.e. σ = σ , σ . Consider the so-called fictitious elastic response (u E , σ E ), i.e. the response of the system if it were purely elastic, defined by

ε E = L : σ E + ε θ , ε E = 1 2 (∇u E + ∇ T u E ), div σ E + f = 0 in Ω , σ E • n = T on ∂ Ω T , u E = v on ∂ Ω u . (5)
The stress field σ can be written as σ = σ E + ρ where ρ is the residual stress field and belongs to the vectorial space H ⊂ E of self-equilibrated fields, defined by

H = {ρ ∈ E : div ρ = 0 in Ω , ρ • n = 0 on ∂ Ω T }. (6)
Let K 0 (t) and K(t) be the convex subsets of E defined as

K 0 (t) = {σ ∈ E : σ (x,t) ∈ C (x,t) ∀x ∈ Ω }, K(t) = K 0 (t) -σ E (t). (7)
The set K 0 (t) is the set of stress fields that are everywhere in the elasticity domain of the material. The set K(t) is the translated of K 0 by -σ E (t). Note that K 0 is independent on time t if the yield parameters are. Under suitable regularity assumptions on ( f , T , v), it can proved that the sets H and K(t) are closed in E [START_REF] Lanchon | Elastic-plastic torsion of a cylindrical shaft with simply or multiply connected cross section[END_REF].

For an arbitrary ρ ∈ K(t), it follows from (1) that

Ω (ρ -ρ) : d(ε -ε E ) dt dω = Ω (ρ -ρ) : d dt (L : ρ) dω + Ω (ρ -ρ) : φ (σ ) dω.
(8) Using (5-6) together with the principle of virtual power shows that the left-hand side of ( 8) is equal to zero. Hence

- Ω (ρ -ρ) : d dt (L(x,t) : ρ) dω = Ω (ρ -ρ) : φ (σ , x,t) dω. (9)
The function

Φ(σ ,t) = Ω φ (σ , x,t) dω
is convex, positive, and vanishes on K 0 (t). It is a classical result [START_REF] Brézis | Opérateurs maximaux monotones et semigroupes de contractions dans les espaces de Hilbert[END_REF] that

Φ(σ ,t) = α 2 σ -P K 0 (t) σ 2 , Φ (σ ,t) = α(σ -P K 0 (t) σ ) (10)
where P K 0 (t) : E → E denotes the projection on K 0 (t) (for the scalar product •, • ). It follows from [START_REF] Weichert | A Historical View on Shakedown Theory[END_REF] that

Ω (ρ -ρ) : φ (σ ) dω = ρ -ρ, Φ (σ ) = ρ -ρ, α(σ -P K 0 (t) σ ) .
Further observing that σ -P K 0 (t) σ = ρ -P K(t) ρ, we obtain from (9) that ρ satisfies

- d dt (L(x,t)ρ) ∈ α(ρ -P K(t) ρ) + H ⊥ (11)
where H ⊥ is the orthogonal of H in E.

Eq. ( 11) can be simplified by projecting it on H. To that purpose, set L(t) = πL where π : E → H is the orthogonal projector on H. Eq. ( 11) becomes

- d dt (L(t)ρ) = α(ρ -πP K(t) ρ). (12)
Since the elastic moduli tensor L(x,t) is symmetric positive definite, it can easily be verified that L(t) is self-adjoint and positive definite. Starting from a given initial state ρ(t = 0), the evolution of the stress field in H is governed by the ordinary differential equation [START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement cyclique[END_REF]. The uniqueness of the stress rate ρ has been proved in [START_REF] Halphen | Elastic perfectly plastic structures with temperature dependent elastic coefficients[END_REF].

Elastic solutions

In the following we study the asymptotic behavior of solutions to (12) as t → ∞. We only consider the case where L(t), K(t) are periodic in time (with the same period T ) and the dimension of H is finite. A central role is played by elastic solutions of (12), i.e. solutions without any plastic yielding. Such an elastic solution ρ * (t) necessarily lies in K(t) at each time t, and therefore satisfies

d dt (L(t)ρ * (t)) = 0, ρ * (t) ∈ K(t) ∩ H ∀t ∈ [0, T ] (13)
Our main objective is to examine conditions under which every solution of ( 12) converges towards an elastic solution in the large time limit. A first requisite is obviously that elastic solutions do exist, i.e. that 0≤t≤T

L(t)K(t) = / 0. ( 14
)
Even when the elastic moduli are time-independent, the condition ( 14) is not sufficient to obtain results on the asymptotic behavior. Loosely speaking, a minimal requirement is the existence of an elastic solution that remains 'strictly inside' the elastic domain. In the case of time-independent elastic moduli, that notion of 'strictly inside' is captured by the Melan's condition defined as follows

Melan's condition (standard form):

There exists an elastic solution ρ * and some

m > 1 such that ρ * (x,t) + mσ E (x,t) ∈ C (x,t) ∀(x,t) ∈ Ω × [0, T ].
Let B(ρ * , r) denotes the ball of center ρ * and radius r in E (for the norm • ). In the following, we consider a strong version of Melan's condition, defined as follows:

Melan's condition (strong version):

There exists an elastic solution ρ * and some r > 0 such that H ∩ B(ρ

* (t), r) ∈ K(t) ∀t ∈ ×[0, T ].
When the dimension of H is finite, the strong and standard versions of Melan's condition are equivalent. This is no longer true in infinite dimension (the strong version, as its names suggests, is more restrictive).

Shakedown theorem

Let us fix some notations and assumptions to be used in the remainder:

• µ(t) denotes an arbitrary non-negative differentiable function such that inf t µ(t) > 0. We set

L 0 (t) = L(t)/µ(t), M 0 (t) = L -1 0 (t) = µ(t)M(t). (15)
with M(t) = L(t) -1 . • The set πK(t) remains bounded 1 i.e. there exists a constant M such that πρ ≤ M for all t and ρ ∈ K(t).

(

• [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF]. Similarly, the operator M(t) is assumed to be bounded. • We introduce a measure γ(0, T ) of the time-fluctuations of the elastic moduli, defined as

γ(0, T ) = T 0 ||| Ṁ0 (t)||| dt = T 0 |||µ(t) Ṁ(t) + μ(t)M(t)||| dt. ( 18
)
The objective of this section is to prove Theorem 1 below.

Theorem 1. If (i) the Melan's condition is satisfied by some (ρ * , r);

(ii) the elastic moduli are such that

γ(0, T ) inf µ < r 2CM ;
then the residual stress converges towards an elastic solution, whatever the initial state is.

In that theorem, condition (ii) sets restriction on the time variations of the elastic moduli. Setting such a restriction is necessary to obtain global convergence results: One can indeed find counterexamples in which condition (i) is fulfilled and some solutions of ( 12) do not become elastic in the large time limit [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF].

The statement of Theorem 1 above is quite similar to that obtained in perfect plasticity [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF]. The proof, however, is more complicated and detailed in the following. To clarify the exposition, the proof is broken down in 3 separate steps, covered by Sect. 3.1-3.3 below. Compared to perfect plasticity, an additional difficulty is that the stress is not restricted to remain in the elasticity domain. For instance, the initial stress can be chosen as arbitrarily large. It can be proved, however, that the stress is bounded at large time: This is the object of Lemma 1, Sect. 3.1. The next step consists in proving that, for large time, the variation of some elastic energy is controlled by the plastic dissipation, in a sense that is defined in Lemma 2, Sect. 3.2. The claimed result follows from those two lemmas, as detailed in Sect. 3.3. In all that follows, ρ(t) denotes an arbitrary solution of (12).

Bound on the stress field

Let η(t) = L(t)ρ(t). The object of the following lemma is to prove that, for large time, η(t) is bounded by some constant M that can be chosen as arbitrarily close to CM.

Lemma 1. For any M > CM, there exists t 0 ≥ 0 such that η(t) ≤ M for all t ≥ t 0 .

Proof. We have

d dt η(t) 2 = -2α ρ -πP K(t) ρ, η(t) = -2α η(t), M(t)η(t) + 2α πP K(t) ρ, η(t) . (19)
Since L(t) is symmetric positive definite, we have |||L(t)||| = max i λ i where {λ i } are the eigenvalues of L(t). The relation [START_REF] Spiliopoulos | A direct method to predict cyclic steady states of elastoplastic structures[END_REF] implies that λ i ≤ C for all i. Since M(t) is positive definite with eigenvalues {1/λ i }, it follows that

η(t), M(t)η(t) ≥ min i 1 λ i η(t) 2 ≥ η(t) 2 C .
Moreover, using Cauchy-Schwarz inequality together with [START_REF] Simon | Shakedown analysis of engineering structures with limited kinematical hardening[END_REF], we find

2α πP K(t) ρ, η(t) ≤ 2α πP K(t) ρ • η(t) ≤ 2αM η(t) .
Substituting in [START_REF] König | A shakedown theorem for temperature dependent elastic moduli[END_REF] gives

d dt η(t) 2 ≤ -2 α C η(t) 2 + 2αM η(t) . (20)
Setting

G(t) = max(M 2 C 2 , η(t) 2
), Eq. ( 20) implies that G is decreasing with time t. It can indeed easily be verified from ( 20) that G is right-differentiable and satisfies G + (t) ≤ 0 for all t, where G + (t) is the right-derivative of G.

If η(t 0 ) ≤ CM for some time t 0 , then it directly follows from the monotonicity of G that η(t) ≤ CM for any t ≥ t 0 , which proves the claim. Now consider the case where η(t) > CM for all time t. Dividing (20) by η(t) , we have 1

η(t) d dt η(t) 2 = 2 d dt η(t) ≤ -2 α C η(t) + 2αM. (21)
Using the differential form of Gronwall's lemma, [START_REF] Hasbroucq | Inelastic responses of a two-bar system with temperature-dependent elastic modulus under cyclic thermomechanical loadings[END_REF] implies that

η(t) ≤ η(t 0 ) e -αt/C +CM(1 -e -αt/C ). (22)
The right-hand side of [START_REF] Hasbroucq | Analytical study of the asymptotic behavior of a thin plate with temperature-dependent elastic modulus under cyclic thermomechanical loadings[END_REF] varies between η(t 0 ) and CM in a monotonic fashion.

Choosing for instance

t 0 = C α | log η(0) -CM M -CM |,
Eq. ( 22) implies that η(t) ≤ M for all t ≥ t 0 , which proves the claim.

Variation of the elastic energy

Consider the positive function f defined as

f (t) = 1 2 τ(t), L 0 (t)τ(t) . ( 23
)
where

τ(t) = µ(t)(ρ(t) -ρ * (t)).
The function f is referred to as the elastic energy. The aim of this section is to establish Lemma 2 below, which bounds the variation f (a + T )f (a) for large time a.

Lemma 2. Let (M ,t 0 ) satisfying Lemma 1. For any a ≥ t 0 we have

f (a + T ) -f (a) ≤ 2M γ(0, T ) -r inf µ a+T a η(t) dt. (24)
Proof. We have

ḟ (t) = τ(t), L 0 (t)τ(t) + 1 2 τ(t), L0 (t)τ(t) . (25)
Recalling that d(L(t)ρ * (t))/dt = 0, Eq.( 25) can be rewritten as ḟ (t) = τ(t), η(t) + H(t) with

H(t) = - 1 2 τ(t), L0 (t)τ(t) . ( 26
)
Integrating on the time interval [a, a + T ], we obtain

f (a + T ) -f (a) = a+T a τ(t), η(t) dt + a+T a H(t) dt. (27)
In the right-hand side of ( 27), the first term is an irreversible contribution associated with the plastic dissipation whereas the second term is a reversible contribution associated with the time fluctuations of the elastic moduli. In the following, we bound those two terms separately, starting with the irreversible contribution -a+T a τ(t), η(t) dt. To that purpose, we use a reasoning that is quite similar to that used in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] for perfect plasticity.

We first note from (12) that η(t) = -πΦ (ρ(t) + σ E (t),t). Since Φ is convex, positive and vanishes in K 0 (t), we have 0 ≤ Φ(ρ(t) + σ E (t),t) ≤ η(t), ρ (t)ρ(t) for all ρ ∈ K(t) ∩ H.

(

) 28
The strong Melan's condition implies that ρ * (t)r η/ η is in K(t) ∩ H. Hence (28) gives

0 ≤ η(t), ρ * (t) -r η η -ρ(t) = η(t), ρ * (t) -ρ(t) -r η(t) , i.e. η(t), ρ(t) -ρ * (t) ≤ -r η(t) . Hence a+T a τ(t), η(t) dt ≤ -r(inf µ) a+T a η dt. (29)
Bounding the reversible contribution a+T a H(t) dt requires a little more effort. We have

τ(t) = L -1 0 (t)(η(t) -η(a) + s) (30)
where

s = L(a)(ρ(a) -ρ * (a)) (31)
is independent on t. Substituting the expression (30) in ( 26) and using the fact that L 0 (t) is self-adjoint, we find

H(t) = - 1 2 (η(t) -η(a) + s), L -1 0 (t) L0 (t)L -1 0 (t)(η(t) -η(a) + s) .
Observing that L -1 0 (t) L0 (t)L -1 0 (t) = -Ṁ0 (t) yields

H(t) = 1 2 (η(t) -η(a) + s), Ṁ0 (t)(η(t) -η(a) + s) . (32)
Since s does not depend on t and M 0 (t) is T -periodic, the integration of (32) on the time interval [a, a + T ] gives

We also have, for t ∈ [a, a + T ],

η(t) -η(a) = t a η(t )dt ≤ t a η(t ) dt ≤ a+T a η(t ) dt . (35)
Substituting (34-35) in (33) and using the definition (18) of γ(0, T ), we obtain

a+T a H(t)dt ≤ 2M γ(0, T ) a+T a η(t) dt (36)
Replacing ( 29) and (36) in [START_REF] Lanchon | Elastic-plastic torsion of a cylindrical shaft with simply or multiply connected cross section[END_REF] gives the desired result.

Proof of the theorem

We are now in a position to prove Theorem 1. By condition (i), one can pick M > CM such that γ(0,T) inf µ < r 2M . Using Lemma 1 we such M , there exists t 0 such that η(t) ≤ M for t ≥ t 0 . Let now N 0 ∈ N be such that N 0 T ≥ t 0 . For i ≥ N 0 , Lemma 2 gives

f ((i + 1)T ) -f (i) ≤ -m (i+1)T iT η(t) dt
where m = r(inf µ)-2M γ(0, T ) is non-negative by (ii). Summing over i = N 0 , • • • , N and recalling that f is positive, we find

NT N 0 T η(t) dt ≤ 1 m f (N 0 T ).
Taking the limit N → ∞ shows that the integral S 0 η(t) dt converges as S -→ +∞. Since E is a Hilbert space, it follows that η(t) also converges towards a limit η ∞ as t -→ +∞.

That limit η ∞ is in H because H is closed in E. Setting ρ ∞ (t) = M(t)η ∞
and recalling that M(t) is bounded, we have ρ(t)ρ ∞ (t) -→ 0 as t → ∞. We now check that ρ ∞ (t) is an elastic solution, i.e. satisfies [START_REF] Carvelli | Shakedown analysis of defective pressure vessels by a kinematic approach[END_REF]. The definition of ρ ∞ (t) rightly gives d(L(t)ρ ∞ (t))/dt = 0, but the fact that ρ ∞ (t) ∈ K(t) calls for some justification. Noting by [START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement cyclique[END_REF]

that η(t) = -α(ρ(t) -πP K(t) ρ(t)), we have ρ ∞ (t)-πP K(t) ρ ∞ (t) ≤ 1 α η(t) + ρ(t)-ρ ∞ (t) + πP K(t) ρ(t)-πP K(t) ρ ∞ (t) .
It is a classical result that the projection operator on a closed convex set is a contractive mapping [START_REF] Brézis | Opérateurs maximaux monotones et semigroupes de contractions dans les espaces de Hilbert[END_REF], hence πP

K(t) ρ(t) -πP K(t) ρ ∞ ≤ ρ(t) -ρ ∞ (t) and ρ ∞ (t) -πP K(t) ρ ∞ (t) ≤ 1 α η(t) + 2 ρ(t) -ρ ∞ (t) . Recall that ρ(t) -ρ ∞ (t) = M(t)(η(t) -η ∞ ). Since M(t) is bounded, there is a constant K such that ρ(t) -ρ ∞ (t) ≤ K η(t) -η ∞ . Therefore, ρ ∞ (t) -πP K(t) ρ ∞ (t) ≤ 1 α η(t) + 2K η(t) -η ∞ (t) .
Integrating on the time interval [iT, (i + 1)T ] and observing that ρ ∞ (t) and K(t) are T -periodic, we obtain

T 0 ρ ∞ (t) -πP K(t) ρ ∞ (t) dt ≤ 1 α (i+1)T iT η(t) dt + 2K (i+1)T iT η(t) -η ∞ dt. Since ∞ 0 η(t) dt < ∞, we have (i+1)T iT η(t) dt -→ 0 as i -→ ∞.
Moreover, since η(t) converges towards η ∞ , we also have

(i+1)T iT η(t) -η ∞ dt -→ 0 as i -→ ∞. It follows that ρ ∞ (t) -πP K(t) ρ ∞ (t) = 0 on [0, T ], i.e. ρ ∞ (t) = πP K(t) ρ ∞ (t) on [0, T ]. Therefore, we have P K(t) ρ ∞ (t) = ρ ∞ (t) + q with q ∈ H ⊥ .
By the definition of the projection it follows that

ρ ∞ (t) -P K(t) ρ ∞ (t), ρ -P K(t) ρ ∞ ≤ 0 for any ρ ∈ K(t). In particular, choosing ρ = ρ * (t) ∈ K(t) ∩ H we find -q, ρ * (t) -ρ ∞ (t) -q = q 2 ≤ 0 hence q = 0, i.e. P K(t) ρ ∞ (t) = ρ ∞ (t) or equivalently ρ ∞ (t) ∈ This completes the proof that ρ ∞ (t) is an elastic solution. Since ρ ∞ (t) -ρ(t) -→ 0 as t → ∞, the convergence of ρ(t)
towards an elastic solution is obtained.

Remark: It is interesting to compare Theorem 1 with the analog result obtained in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] for elastic-perfectly plastic materials. In [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF], a condition analog to (ii) was formulated, with the factor 2 in (ii) replaced by a factor 3. That factor 3 results from Eq. ( 46) in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF], which -using the present notations -states that η(t)η(a) + 2s ≤ 6CM. Using a reasoning similar to that used in Eq. (34), one can actually observe that η(t)-η(a)+2s ≤ 4CM, which leads to an improved factor 2 instead of a factor 3.

There is also a more subtle difference between Theorem 1 above and the results in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF], regarding the definition of the constant M. In Theorem 1, M is a bound on πK(t) whereas in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF] M is a bound on K(t) ∩ H. We have K(t) ∩ H ⊂ πK(t) but the inclusion is generally strict.

Illustrative example

As an illustration of Theorem 1, consider the problem of an elastic-viscoplastic plate under cyclic thermal loading. The plate is stress-free in the e 3 direction and is in frictionless contact with rigid walls in the (e 1 , e 2 ) directions (Fig. 1). Length units are chosen such that the plate has unit volume. The constitutive material is homogeneous and the dissipation potential φ of the Perzyna type

φ (σ ) = 1 α J 2 -σ y 2 + ( 37
)
where α is a viscosity parameter, σ y is the yield strength, J 2 = 3/2|σ -(tr σ /3)I| is the second invariant of the deviatoric stress, and x + = max(0, x) denotes the positive part of a scalar x. The potential (37) can be put in the format (3) by setting C (x,t) = {σ : J 2 ≤ σ y }.

Because of the rigid walls in the e 1 and e 2 directions, thermal dilatation may result in high compressive stress and plastic flow. In the following we are interested in bounding the temperature fluctuations θ for which shakedown occurs.

Temperature-independent elastic moduli

Using the plane stress assumption, we consider a 2-dimensional model of the problem. If the imposed temperature field θ as well as the initial state ε p (t = 0) are uniform -which is assumed in the following -then the fields ε p , ε and σ remain uniform at all time. The space E is thus chosen as the 3-dimensional space of tensors σ with a matrix representation of the form

σ =   σ 11 σ 12 0 σ 12 σ 22 0 0 0 0   (38)
in the basis (e 1 , e 2 , e 3 ). The space H of residual stresses is the subspace of E constituted by uniform fields ρ of the form

ρ =   ρ 1 0 0 0 ρ 2 0 0 0 0   . (39)
The purely elastic stress response σ E (t) of the plate is given by

σ E =   f E (t) 0 0 0 f E (t) 0 0 0 0   (40)
where

f E (t) = -E(t)ε θ (t)/(1 -ν(t))
. It can easily be verified that

K 0 = {σ ∈ E : σ 2 11 + σ 2 22 -σ 11 σ 22 + 6σ 2 12 ≤ σ 2 y }, (41)
πK 0 = {ρ ∈ H : ρ 2 1 + ρ 2 2 -ρ 1 ρ 2 ≤ σ 2 y }. (42)
Let (E 1 , E 2 ) be the orthonormal basis of H defined by E 1 = diag(1, 0, 0) and E 2 = diag(0, 1, 0). The set πK 0 in (42) is a solid ellipsoid with axes

E 1 + E 2 and E 1 -E 2 .
The set πK(t) = π(K 0σ E (t)) is obtained by a time-dependent translation of πK 0 in the E 1 + E 2 direction (Fig. 2). Setting σ θ = sup t |σ E (t)|, it can easily be verified that ∩ t πK(t) is non empty as long as σ θ ≤ σ y . More precisely, ∩ t πK(t) contains B(0, r) ∩ H where B(0, r) is the ball centered at the origin with a radius r given by

r =          σ y 2 3 - σ θ σ y 2 if 0 ≤ σ θ σ y ≤ 2 3 , σ y √ 2(1 - σ θ σ y ) if 2 3 ≤ σ θ σ y ≤ 1. ( 43
)
Since ρ = 0 is an elastic solution to the evolution problem, it can be deduced from Theorem 1 (or from the standard form of Melan's theorem) that shakedown occurs if σ θ ≤ σ y .

ρ 1 ρ 2
Fig. 2 The set πK(t) is a solid ellipsoid obtained by translation of the ellipsoid πK 0 (shown in dotted lines) in the (1,1) direction. Shown in solid lines are the extreme locations of πK(t), corresponding to f E (t) = ±σ θ . The filled ball centered at the origin is included in πK(t) for all t.

Temperature-dependent elastic moduli

Let us now use Theorem 1 to estimate the shakedown limit in the case of non constant elastic moduli. Applying Theorem 1 requires to evaluate the constants M, C in [START_REF] Simon | Shakedown analysis of engineering structures with limited kinematical hardening[END_REF][START_REF] Spiliopoulos | A direct method to predict cyclic steady states of elastoplastic structures[END_REF] as well as the scalar γ(0, T ) defined in [START_REF] Borino | Consistent shakedown theorems for materials with temperature dependent yield functions[END_REF]. It can be verified from (42) that any σ in πK 0 satisfies σ ≤ √ 2σ y . Hence any σ in πK(t) satisfies

σ ≤ √ 2σ y + σ E (t) ≤ √ 2(σ y + σ θ ).
The constant M in ( 16) can thus be chosen as

M = √ 2(σ y + σ θ ). (44)
Recall that L(t)ρ is defined as the projection of L(t) : ρ on H. In the present case, we have

L(t)ρ = 1 E(t)   ρ 1 -νρ 2 0 0 0 ρ 2 -νρ 1 0 0 0 0   .
Hence the matrix representation of L in the basis (E 1 , E 2 ) is given by

L(t) = 1 E(t) 1 -ν -ν 1 .
The

In order to calculate γ(0, T ), we note that

M(t) = L -1 (t) = E(t) 1 -ν 2 (t) 1 ν ν 1 .
Choosing µ(t) = 1/ tr M(t) as suggested in [START_REF] Peigney | Shakedown of elastic-perfectly plastic materials with temperaturedependent elastic moduli[END_REF], we get

µ(t) = 1 -ν 2 (t) 2E(t) , M 0 (t) = 1 2 1 ν ν 1 , Ṁ0 (t) = 1 2 0 ν ν 0 . It follows that ||| Ṁ0 (t)||| = 1 2 | ν| and γ(0, T ) = 1 2 T 0 | ν|dt.
Using the values of M and C defined in (44) and ( 46) respectively, we obtain from Theorem 1 that shakedown occurs if

1 + ν max E min T 0 | ν|dt inf t 1-ν 2 (t) E(t) ≤ r √ 2(σ y + σ θ ) ( 47
)
where r is the function of σ θ defined in (43). The left-hand side of ( 47) is a function of the fluctuations of the elastic moduli, whereas the right-hand side is a function of the loading σ θ . Observing that (1ν 2 (t))/E(t) ≥ (1ν 2 max )/E max , it can be seen that a sufficient condition for Eq. ( 47) to be satisfied is

1 1 -ν max E max E min T 0 | ν|dt ≤ r √ 2(σ y + σ θ ) . (48)
Assume that ν vanishes only when ν(t) = ν min or ν(t) = ν max , i.e. ν(t) grows monotonically from ν min to ν max and then decreases monotonically from ν max to ν min , in a periodic fashion. In such case, we have

T 0 | ν|dt = 2(ν max -ν min )
and the inequality (48) reduces to

(ν max -ν min ) 1 -ν max E max E min ≤ r 2 √ 2(σ y + σ θ ) . (49)
For fixed values of E max /E min , the shakedown limit σ θ SD is defined at the largest value of σ θ that satisfies (49). The shakedown limit σ θ SD is plotted in Fig. 3 as a function of E max /E min , for several values of ν maxν min . The maximum value ν max is set to 0.3 for all the curves in Fig. 3. First consider the case ν max -ν min = 0 : Whatever the value of E max /E min , the shakedown limit is equal to σ y and thus coincides with the value obtained in Sect. 4.1 for temperature-independent elastic moduli. When ν maxν min increases, the shakedown limit decreases rapidily. For instance, the shakedown limit is approximatively equal to 0.6σ y in the case of a 20% variation of E and ν around nominal values (E 0 , ν 0 = 0.3). For common metals, the Young modulus E decreases with the temperature, while the Poisson ratio increases [START_REF] Ledbetter | Elastic Properties of Metals and Alloys. II. Copper[END_REF]. The Young modulus is more sensitive to temperature variations than the Poisson ratio. To put things in perspective, a 20% variation of the Young modulus E in steels would typically correspond to temperatures fluctuating between 0 ˚C and 500 ˚C [START_REF]Eurocode 3: Design of Steel Structures, ENV 1993-1-2: General Rules -Structural Fire Design[END_REF].

Incremental analysis

As a further illustration of Theorem 1, we use incremental analysis to solve [START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement cyclique[END_REF] for given loading history and fluctuations of the elastic moduli. The functions f E (t) , E(t) and ν(t) are taken as f E (t) = σ θ sin ωt, E(t) = E 0 (1 -0.2 sin 2 ωt), ν(t) = ν 0 (1 + 0.2 sin 2 ωt) with E 0 /σ y = 10 and ν 0 = 0.3. For such parameters, the shakedown limit σ θ SD as provided by Eq. ( 49) is approximatively equal to 0.6σ y . The differential equation ( 12) is solved numerically using a Runge-Kutta scheme with variable step size. In Fig. 4 is shown the evolution of the plastic strain in the case σ θ = 0.5σ y , i.e. below the shakedown limit. The initial state is taken as ρ 1 (t = 0) = σ y , ρ 2 (t = 0) = 2σ y . In accordance with the results of Sect. 4.2, the plastic strain converges towards a limit as t → ∞. Theorem 1 ensures that such behavior occurs for all inital states ρ 1 (t = 0), ρ 2 (t = 0). In Fig. 5 is plotted the evolution of ε p (t) (solid line). A fast decrease is observed for small time: This is the meaning of Lemma 1 introduced in Sect. 3.1. In dashed line is plotted the exponental upper bound on ε p (t) that is deduced from Eq. ( 22) .

Concluding remarks

The shakedown theorem presented in this chapter gives a sufficient condition for the evolution to become elastic in the large time limit, whatever the initial state is. Loosely speaking, that theorem states that if the loading is not too large (in the sense of condition (i)) and the time-fluctuations of the elastic moduli are not too large (in the sense of condition (ii)), then shakedown occurs immaterial of the initial state. We emphasize again that setting a restriction on the time-fluctuations of the elastic moduli is essential to ensure that shakedown occurs in a path-independent fashion. It can be observed that Theorem 1 is independent on the viscosity parameter α. This is consistent with the case of constant elastic moduli: the shakedown behavior is essentially determined by the elasticity domain of the material (as long as the dissipation φ satisfies standard assumptions such as convexity).

For the sake of simplicity, the viscosity parameter α has been assumed to be independent on (x,t) but there is no difficulty in extending Theorem 1 to non constant viscosities. It could be interesting to investigate whether Theorem 1 could be extended to other viscoplastic potentials than those of the form (3).

In practice, the shakedown theorem that has been presented could be useful for the fatigue design of structures submitted to large temperature variations. As illustrated in Sect. 4, that theorem gives lower bounds on the shakedown limit. It would be interesting to investigate whether the kinematic shakedown theorem of [START_REF] Koiter | General theorems for elastic-plastic solids[END_REF] could be extended to the case of temperature-dependent elastic moduli, so as to obtain upper bounds.

  )η(a), Ṁ0 (t)(η(t)η(a) + 2s) dt. (33) Using now the Cauchy-Schwarz inequality and the definition of the norm operator, we obtain η(t)-η(a), Ṁ0 (t)(η(t)-η(a)+2s) ≤ η(t)-η(a) •||| Ṁ0 (t)|||• η(t)-η(a)+2s . Since η(t)η(a) + 2s = η(t) + η(a) -2L(a)ρ * (a), we have η(t)η(a) + 2s ≤ η(t) + η(a) + 2 L(a)ρ * (a) .

Lemma 1

 1 gives η(t) ≤ M and η(a) ≤ M . Since ρ * (a) ∈ K(a) ∩ H, we have L(a)ρ * (a) ≤ CM as a consequence of (16-17). Hence η(t)η(a) + 2s ≤ 2M + 2CM ≤ 4M .

Fig. 1

 1 Fig.1An elastic-viscoplastic plate under a cyclic temperature θ (t). The plate is constrained in the e 1 and e 2 directions.

Fig. 3

 3 Fig. 3 Shakedown limit as a function of E max /E min , plotted for several values of ν maxν min . Case ν max = 0.3.

Fig. 4 22 Fig. 5

 4225 Fig.4Evolution of the plastic strain for a loading below the shakedown limit. Case ρ 1 (t = 0) = σ y , ρ 2 (t = 0) = 2σ y .

  operator L(t) being symmetric, its norm operator |||L(t)||| is equal to max i |λ i | where {λ i } are the eigenvalues of L(t). A simple calculation gives

		|||L(t)||| =	1 + ν(t) E(t)	.	(45)
	Set
	ν min = inf t	ν(t), ν max = sup t	ν(t), E min = inf t	t E(t), E max = sup	E(t).
	From (45) we have |||L(t)||| ≤ 1+ν max E min chosen as C = for all t. Hence the constant C in (17) can be 1 + ν max E min .

It can observed that πK(t) is bounded if the C (x,t) is.


O'tega Ejofodomi

Michael Olawuyi

Don Uche Onyishi

Godswill Ofualagba

Detecting Architectural Distortion In Mammograms Using A Gabor Filtered Probability Map Algorithm

Keywords: Breast Cancer, Architectural Distortion, Gabor Filters 1

Breast Cancer is a disease that is prevalent in many countries. Computer-Aided detection (CAD) systems have been developed to assist radiologists in detecting breast cancer. This paper discusses an algorithm for architectural distortion (AD) detection with a better sensitivity than the current CAD systems.

19 images containing ADs were preprocessed with a median filter and Gabor filters to extract texture information. AD probability maps were generated using a maximum amplitude map and histogram analysis on the orientation map of the Gabor filter response. AD maps were analyzed to select ROIs as potential AD sites.

AD map analysis yielded a sensitivity of 79% (15 out of 19 cases of AD were detected) with a false positive per image (FPI) of 18. Future work involves the development of a second stage in the algorithm to reduce the FPI value and application of the algorithm to a different set of database images.

INTRODUCTION

Breast cancer is a disease that is prevalent in many countries. In the United States of America, breast cancer is the most common type of cancer, side from skin cancer. The American Cancer Society estimates that there will be 234, 580 new cases of breast cancer (232,340 female, 2,240 male) and that 40,030 people (39,620 female, 410 male) will die from breast cancer in 2013 [START_REF]Cancer Facts and Figures[END_REF]. In the United Kingdom, it is estimated that there were 49,961 new cases (49,564 female, 397 male) of breast cancer, and 11, 633 deaths [START_REF] Sampat | Computer-aided detection and diagnosis in mammography[END_REF]556 female, 77 male) related to breast cancer in 2010 [2]. In some countries, the statistics for breast cancer is unreliable.

The current gold standard for detecting breast cancer is screening mammography, which is X-ray imaging of the breast. Mammography can often detects breast cancer at an early stage, when treatment is more effective and a cure is more likely, especially for women over 50 years. This reduces the mortality rate of the disease. Numerous studies have shown that early detection with mammography saves lives and increases treatment options [START_REF]Cancer Facts and Figures[END_REF][START_REF] Magnus | Effectiveness of mammography screening in reducing breast cancer mortality in women aged 39-49 years: a meta-analysis[END_REF][START_REF] Mandelblatt | Effects of mammography screening under different screening schedules: model estimates of potential benefits and harms[END_REF]. In mammography, two types of X-ray images are usually obtained: the Mediolateral Oblique (MLO) view, which is a sideways view of the breast, and the Craniocaudal view (CC), which is a top to bottom view of the breast.

Radiologists search each image for signs of lesions and abnormalities that may represent breast cancer. Typically, they search for masses, calcifications and architectural distortions. Breast lesions are described and reported according to the Breast Imaging Reporting and Data System (BIRADS ™), which is a terminology developed by the American College of Radiology (ACR), for the description of mammographic lesions [START_REF]ACR BI-RADS -Mammography, Ultrasound & Magnetic Resonance Imaging[END_REF]. A mass is defined as a space-occupying lesion seen in at least two different projections. Calcifications are calcium deposits in the breasts that show up as white, bright spots in mammograms. Architectural distortions (AD) are abnormal breast lesions in which the breast parenchyma is distorted with no visible mass. The distortion includes spiculations radiating from a point, and focal retraction or distortion of the edge of the parenchyma [START_REF]ACR BI-RADS -Mammography, Ultrasound & Magnetic Resonance Imaging[END_REF].

Breast cancer detection using mammograms is far from perfect. It is estimated that radiologists fail to detect about 10-30% of breast cancers [START_REF] Kerlikowske | Performance of screening mammography among women with and without a first-degree relative with breast cancer[END_REF][START_REF] Kolb | Comparison of the performance of screening mammography, physical examination, and breast US and evaluation of factors that influence them: an analysis of 27,825 patient evaluations[END_REF][START_REF] Bird | Analysis of cancers missed at screening mammography[END_REF]. Computer-Aided Detection (CAD) systems has been developed to enable radiologists detect lesions that may be indicative of cancer. These systems serve as support to the radiologist and are not used by themselves as a breast cancer detection tool. Commercial CAD systems accurately identified 90% cases of masses and microcalcifications [START_REF]The CADx Second Look system[END_REF][START_REF]The R2 Technology's Image Checker[END_REF]. However, the sensitivity for AD is low. The R2 Image Checker system successfully identified 49% cases containing ADs with 0.7 false positives per image. The CADx Second Look system [START_REF]The CADx Second Look system[END_REF] successfully detected 33% cases of ADs with 1.27 false positives per image [START_REF] Sampat | Computer-aided detection and diagnosis in mammography[END_REF].

There is need to improve the detection of ADs in CAD systems. The goal of this work is to develop an algorithm that can detect ADs with higher sensitivity than the currently existing CAD systems. Details of the algorithm are presented in the methods section. The results of the algorithm are shown in the results section. Analysis of the performance of the algorithm is provided in the discussion section. Future improvements to the algorithm are highlighted in the conclusion section.

MATERIALS AND METHODS

MATERIALS

All of the images used in this study were obtained from the Mammographic Image Analysis Society (MIAS) Digital Mammogram Database [START_REF] Suckling | The Mammographic Image Analysis Society Digital Mammogram Database[END_REF]. All the images were digitized at 50 micron pixel edge. The sizes of all images were 1024 by 1024 pixels. There were a total of 19 images containing ADs. The pixel coordinates of the AD sites, the pathology of the AD, and the nature of the background breast tissue are also provided. Image views were both MLO and CC views of the breast. These 19 images were used in this pilot study. Figure 1 shows a MIAS image with a circle marking the site of the AD. The ground truth provided by the database was used in determining the effectiveness of the algorithm.

METHODS

Most of the reported automated detection of AD in literature consists of several steps. There may be some preprocessing to remove noise from the mammograms and to select the breast region within the mammogram. The next stage involves the extraction of the texture orientation (line structures oriented at different angles within the image). The texture orientation is then analyzed to produce regions of interest (ROI) which are potential AD sites. Estimation of the characteristic features of ADs are then used to narrow down the number of potential AD sites, thereby reducing the number of false positives in automated detection. The methodology employed in this work consists of four major steps: preprocessing, Gabor filtering to extract texture information, generation of AD maps, and finally, analysis of AD maps to select ROIs. The flow chart for the algorithm is given in Figure 2.

Image Preprocessing

The images were filtered with a median filter to reduce the noise in the image. No region growing methods were employed to segment the breast region from the background.

Fig. 2. Flow Chart for AD probability map algorithm

Gabor Filtering

There are many different techniques employed to extract the texture orientation from an image. The method used in this study is the Gabor Filters. The Gabor filter is a sinusoidally modulated Gaussian function. Mathematically, a 2D Gabor function g, is the product of a 2D Gaussian and a complex exponential function. The general expression is given by: g θ , λ,σ1,σ2 (x, y) = exp -1/2{x y}M {x y} T exp{jπ (x cosθ + y sinθ) [START_REF]Cancer Facts and Figures[END_REF] λ

_ where M = diag(,σ 1 -2
,σ 2 -2 ). The parameter θ represents the orientation, λ is the wavelength, and σ 1 and σ2 represent scale at orthogonal directions. When the Gaussian part is symmetric, we obtain the isotropic Gabor function: gθ,λ,σ(x, y) = exp{-x2 + y2/2σ 2 ) exp { jπ (x cosθ + y sinθ)

(2) λ A bank of Gabor filters was applied to the images. The angle of orientation in the bank of filters was varied from θ ϵ [-π/2;π/2], the frequency of sinusoidal function was set to f = 0.022, and the variance S was set to 0.05. These parameters were chosen experimentally. A total of 19 Gabor filters were applied to the images.

AD Map Generation

Two maps were created from the output of the Gabor filtration: a maximum amplitude map A(x, y) and a dominant orientation map O(x,y). The maximum amplitude for each point (x,y) in the responses to the Gabor filters was extracted to form the amplitude map A(x,y). The angle of orientation of the maximum amplitude was considered to be the dominant angle and was used in creating the orientation map O(x,y). Both maps are used to create the final Architectural Distortion map AD(x,y).

To create the AD map, histogram analysis was performed on O(x,y). Both the amplitude and the orientation map were divided into blocks of 16 by 16 pixels. The maximum amplitude for each block was obtained from A(x,y). Histogram analysis was performed on each block of O(x,y). It is expected that blocks containing ADs should show several dominant peaks (corresponding to dominant orientation angles), while blocks of normal tissue should only have one dominant peak. Therefore, in the histogram analysis an experimentally determined threshold T 1 was set to determine if a peak is dominant or not. If a peak was within a certain percentage of the maximum orientation angle, it was classified as a dominant orientation angle.

The AD map was created by multiplying the maximum grayscale amplitude in block with the number of dominant angles and then dividing by a characteristic number N. The characteristic number was established experimentally. Mathematically,

AD i = Dn i * A max (3) N

Where i is the index number of the block, Dn is the number of dominant angles in block i, A max is the maximum grayscale amplitude in block i, and N is the characteristic number. AD is a value that quantifies the probability of the presence of AD in block i [START_REF] Jasionowska | A two-step method for detection of architectural distortions in mammograms[END_REF].

AD Map Analysis

In the final stage, the AD maps were analyzed. The blocks with the highest AD values were selected as regions of interest (ROI). The image containing ROIs was then compared with the ground truth mammogram. A true positive (TP) was recorded when a ROI was within the region identified as AD by the ground truth. A false positive (FP) was recorded when a ROI was not within the region identified as AD by the ground truth.

RESULT

Our data consisted of the 19 images identified with AD in the MIAS database. Table 1 shows the number of true positives (TP) and false positives (FP) identified in each image. Images in which ADs were unidentified have 0 TPs. In one case (Image 18) the AD has been marked twice (2 TPs). We achieved a sensitivity of 78.9% (15 cases of AD were detected) with a false positive per image (FPI) of 18.4 (see Table 2). 10 out of 19 ADs were classified as malignant, while 9 was classified as benign AD in the ground truth data. Of the 10 malignant ADs, the AD Map analysis correctly identified 80% (8 out of 10) malignant ADs, and 78% (7 out of 9) of benign ADs. The total number of ROIs obtained was 291, of which 16 ROIs positively identified AD. The results of the AD Map analysis algorithm are compared with two commercial CAD systems in Table 2. The sensitivity for the AD Map analysis algorithm is much higher than that of the commercial CAD systems. However, the FPI is also relatively high. Most AD detection algorithms employ two stages: Identification of ROIs in the first stage and reduction of the number of FPs in the second stage [START_REF] Jasionowska | A two-step method for detection of architectural distortions in mammograms[END_REF]14]. The AD map analysis algorithm is currently a one stage detection algorithm. A second stage would need to be included to reduce the FPI to a more reasonable number. shown in Figure 3b.The maximum amplitude map A(x,y) obtained after Gabor filtering is shown in Figure 3c. The histogram of the orientation map O(x,y) is shown in Figure 3d. The generated AD map is shown in Figure 3e. The selected ROI from the AD map analysis is shown in Figure 3f. Finally, the selected ROIs overlaid over the original image is shown in Figure 3f. The ground truth is also labeled in Figure 3f. In this particular case, one of the ROIs is located within the region identified as AD by the ground truth.

DISCUSSION

The sensitivity of the Gabor filtered probability map algorithm is far greater than that of the two commercial CAD systems shown in Table 2. The high number of FPI can be attributed to several factors. First as was duly noted, no segmentation to separate the breast region was performed on the images in the preprocessing step. Nonidentification of the breast boundary (skin-air) caused the algorithm to flag ROIs outside of the breast region. This influenced the number of FPI in each image. The boundary between breast and film showed up consistently as a bright vertical line in the AD maps (see Figure 3e). This had to be excluded in the AD map analysis. Secondly, the pectoral muscle in the MLO view shows up as a bright white region. Since the algorithm relies on the fact that AD sites possess maximum grayscale values, the pectoral muscle region is often flagged as a ROI even when AD is not present. Thirdly, images in the MIAS database contain labels identifying the image as either an MLO view or CC view (see Figure 1a). These labels were imprinted on the films at different locations and so the algorithm did not remove them in the preprocessing step. These also contributed to the high FPI.

There are several things that can be done to improve the AD probability map algorithm. Most importantly, a second stage needs to be incorporated into the algorithm to reduce the number of FPI. This can be accomplished by filtering and thresholding in the Radon Domain, feature extraction and classification of ROIs [START_REF] Suckling | The Mammographic Image Analysis Society Digital Mammogram Database[END_REF][START_REF] Jasionowska | A two-step method for detection of architectural distortions in mammograms[END_REF]. The number of filters in the Gabor banks can be increased to include more orientation angles. A breast segmentation algorithm can also be implemented in the preprocessing stage to extract the breast region from the surrounding to reduce the FPI due to the labels and the breast-film boundary. Finally, it should be noted that the algorithm has been applied to a single dataset from the MIAS database. The results could vary if applied to a different set of images.

CONCLUSION

We have developed an algorithm to detect ADs in mammograms. The Gabor filtered probability map algorithm has a sensitivity of 79% and an FPI of 18, more than two times that of commercial CAD systems. Future work involves the development of a second stage in the algorithm to reduce the FPI value and application of the algorithm to a different set of database images.

Fig. 1 .

 1 Fig. 1. (a) Original MIAS image (b) Ground truth: Original image with a boundary containing AD overlaid. (b) Enlarged ground truth image

Fig. 3 .

 3 Fig. 3. AD Detection Algorithm. (a) Original MIAS image (b) Preprocessed Image (c) Maximum Amplitude Map A(x,y) (d) Orientation Map O(x,y) (e) AD map (f) ROIs with highest AD values (g) AD detection results overlaid over original image

Table 1 .

 1 True Positives (TP) and False Positives (FP) identified in Images

	Image	No. of True Positives	No. of False positives
	No.	(TP)	(FP)
	1	1	41
	2	1	5
	3	1	5
	4	1	10
	5	1	17
	6	0	>> 41
	7	1	10
	8	1	17
	9	1	40
	10	1	40
	11	1	33
	12	1	13
	13	1	19
	14	0	>> 41
	15	0	>> 41
	16	1	11
	17	0	>>41
	18	2	7
	19	1	7

Table 2 .

 2 Comparison of AD Map Analysis Algorithm to Commercial CAD systems.

	AD Map Analy-	Image Checker	CADx	Second
	sis		R 2	Look
	Sensitivity (%)	79	38	21
	No of Images	19	45	45
	FPI	18	0.7	1.27
	% of malignant ADs	80%	-	-
	detected
	% of benign ADs de-	70%	-	-
	tected
	Total No of ROIs	291	-	-
	Total No of ROIs	16	-	-
	with AD
	(a)	(b)		(c)
	(d)	(e)	(f)	(g)