Ezequiel Geremia

Olivier Clatz

Bjoern H Menze

Ender Konukoglu

Antonio Criminisi

Nicholas Ayache

Spatial Decision Forests for MS Lesion Segmentation in Multi-Channel Magnetic Resonance Images

Keywords: Multi-sequence MRI, Segmentation, Multiple Sclerosis, Random Forests

A new algorithm is presented for the automatic segmentation of Multiple Sclerosis (MS) lesions in 3D Magnetic Resonance (MR) images. It builds on a discriminative random decision forest framework to provide a voxel-wise probabilistic classification of the volume. The method uses multi-channel MR intensities (T1, T2, FLAIR), knowledge on tissue classes and long-range spatial context to discriminate lesions from background. A symmetry feature is introduced accounting for the fact that some MS lesions tend to develop in an asymmetric way. Quantitative evaluation of the proposed methods is carried out on publicly available labeled cases from the MICCAI MS Lesion Segmentation Challenge 2008 dataset. When tested on the same data, the presented method compares favorably to all earlier methods. In an a posteriori analysis, we show how selected features during classification can be ranked according to their discriminative power and reveal the most important ones.

Introduction

Multiple Sclerosis (MS) is a chronic, inflammatory and demyelinating disease that primarily affects the white matter of the central nervous system.

Automatic detection and segmentation of MS lesions can help diagnosis and patient follow-up. It offers an attractive alternative to manual segmentation which remains a time-consuming task and suffers from intra-and inter-expert variability. MS lesions, however, show a high variability in appearance and shape which makes automatic segmentation a challenging task. MS lesions lack common intensity and texture characteristics, their shapes are variable and their location within the white matter varies across patients.

A variety of methods have been proposed for the automatic segmentation of MS lesions. For instance, in [START_REF] Anbeek | Probabilistic Segmentation of White Matter Lesions in MR Imaging[END_REF] and [START_REF] Admiraal-Behloul | Fully automatic segmentation of white matter hyperintensities in mr images of the elderly[END_REF], the authors propose to segment white matter signal abnormalities by using an intensity-based k-nearest neighbors method with spatial prior and a fuzzy inference system, respectively. A similar classifier combined with a template-driven segmentation was proposed in [START_REF] Wu | Automated segmentation of multiple sclerosis lesion subtypes with multichannel mri[END_REF] to segment MS lesions into three different subtypes (enhancing lesions, T1 black holes, T2 hyperintense lesions). A false positive reduction based on a rule-based method, a level set method and a support vector machine classifier is presented in [START_REF] Yamamoto | Computer-aided detection of multiple sclerosis lesions in brain magnetic resonance images: False positive reduction scheme consisted of rule-based, level set method, and support vector machine[END_REF] along with a multiple-gray level thresholding technique. Many general purpose brain tissue and brain tumor segmentation approaches can be modified easily for MS lesion segmentation.

In (Bricq et al., 2008b), for example, the authors present an unsupervised algorithm based on hidden Markov chains for brain tissue segmentation in 2 MR sequences. The method provides an estimation of the proportion of white matter (WM), grey matter (GM) and cerebro-spinal fluid (CSF) in each voxel. It can be extended for MS lesions segmentation by adding an outlier detector (Bricq et al., 2008a).

Generative methods were proposed consisting in a tissue classification by means of an expectation maximization (EM) algorithm. For instance, the method presented in [START_REF] Datta | Segmentation and quantification of black holes in multiple sclerosis[END_REF] aims at segmenting and quantifying black holes among MS lesions. The EM algorithm can be modified to be robust against lesion affected regions, its outcome is then parsed in order to detect outliers which, in this case, coincide with MS lesions [START_REF] Van Leemput | Automated segmentation of multiple sclerosis lesions by model outlier detection[END_REF]. Another approach consists in adding to the EM a partial volume model between tissue classes and combining it with a Mahalanobis thresholding which highlights the lesions [START_REF] Dugas-Phocion | Hierarchical segmentation of multiple sclerosis lesions in multi-sequence MRI[END_REF]. Morphological postprocessing on resulting regions of interest was shown to improve the classification performance [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. In [START_REF] Freifeld | Multiple sclerosis lesion detection using constrained GMM and curve evolution[END_REF], a constrained Gaussian mixture model is proposed, with no spatial prior, to capture the tissue spatial layout. MS lesions are detected as outliers and then grouped in an additional tissue class. Final delineation is performed using probability-based curve evolution. Multi-scale segmentation can be combined with discriminative classification to take into account regional properties [START_REF] Akselrod-Ballin | An integrated segmentation and classification approach applied to multiple sclerosis analysis[END_REF]. Beyond the information introduced via the spatial prior atlases, these methods are limited in their ability to take advantage of long-range spatial context in the classification task.

To overcome this shortcoming, we propose the use of an ensemble of dis-criminative classifiers. Our algorithm builds on the random decision forest framework which has multiple applications in bioinformatics [START_REF] Menze | A comparison of random forest and its Gini importance with standard chemometric methods for the feature selection and classification of spectral data[END_REF], and, for example, more recently also in the image processing community [START_REF] Andres | Segmentation of SBFSEM volume data of neural tissue by hierarchical classification[END_REF][START_REF] Yi | Discriminative, semantic[END_REF][START_REF] Criminisi | Regression forests for efficient anatomy detection and localization in CT studies[END_REF]. Adding spatial and multi-channel features to this classifier proved effective in object recognition [START_REF] Shotton | Textonboost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[END_REF], brain tissue segmentation in MR images [START_REF] Yi | Discriminative, semantic[END_REF], myocardium delineation in 3D echocardiography [START_REF] Lempitsky | Random forest classification for automatic delineation of myocardium in real-time 3D echocardiography[END_REF] and organ localization in CT volumes [START_REF] Criminisi | Regression forests for efficient anatomy detection and localization in CT studies[END_REF].

Applying multi-channel and context-rich random forest classification to the MS lesion segmentation problem is novel, to our knowledge. The presented classifier also exploits a specific discriminative symmetry feature which stems from the assumption that the healthy brain is approximately symmetric with respect to the mid-sagittal plane and that MS lesions tend to develop in asymmetric ways. We then show how the forest combines the most discriminative channels for the task of MS lesion segmentation.

Materials

This section describes the data, algorithms and notations which are referred to in the rest of the article.

MICCAI Grand Challenge 2008 dataset

The results in this article rely on a strong evaluation effort. This section presents the MICCAI 1 Grand Challenge 2008 datasets, which is the largest 1 MICCAI is the annual international conference on Medical Image Computing and Computer Assisted Intervention.

dataset publicly available, and explains the way our method is compared against the winner of the challenge [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. In the rest of the article, the MICCAI Grand Challenge 2008 on MS Lesions Segmentations will be referred as MSGC.

Presentation

The MSGC (Styner et al., 2008a) aims at evaluating and comparing algorithms in an independent and standardized way for the task of MS lesion segmentation. The organizers make publicly available two datasets through their website. A dataset of labeled MR images which can be used to train a segmentation algorithm, and an unlabeled dataset on which the algorithm should be tested. The website offers to quantitatively evaluate the segmentation results on the unlabeled dataset using the associated private ground truth database, and to publish the resulting scores. This project is an original initiative to provide an unbiased comparison between MS lesions segmentation algorithms. In the rest of the article, the dataset for which labels are publicly available will be referred to as public dataset, whereas the dataset for which data is not available will be referred to as private dataset.

Data

The public dataset contains 20 cases, 10 from the Children's Hospital in Boston (CHB) and 10 from the University of North Carolina (UNC), which are labeled by a CHB expert rater. The private dataset contains 25 cases, 15 from CHB and 10 from UNC. The private dataset was annotated by a single expert rater at CHB and jointly by 2 expert raters UNC. For each case, the centers provided 3 MR volumes: a T1-weighted image, a T2-weighted image and a FLAIR image. These were co-registered and sampled to fit the isotropic 0.5 × 0.5 × 0.5 mm 3 resolution.

Both private and public datasets gather anatomical images from two different centers, CHB and UNC, and shows high variability in intensity contrast (cf. Section 5.2), image noise and bias field. Both public and private datasets contain highly heterogeneous cases and could thus be considered as realistic test cases.

Evaluation

Quantitative evaluation is carried out on the private dataset using a set of known metrics defined in (Styner et al., 2008a) and summed up in Table 1.

The two full sets of expert segmentations were used as reference for method comparison.

Top-ranked methods

The challenge results highlight four top-ranked methods each reflecting a different approach to the task of MS lesion segmentation. A k-nearest neighbor classification of brain tissue relying on spatial location and intensity value was proposed in [START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF]. This method provides a voxel-wise probabilistic classification of MS lesions. In (Bricq et al., 2008a), the authors present an unsupervised segmentation algorithm based on a hidden Markov chain model. The method takes into account neighborhood information, MR sequences and probabilistic priors in order to delineate MS lesions. Alternatively, the iterative method proposed in [START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF][START_REF] Shiee | A topology-preserving approach to the segmentation of brain images with multiple sclerosis lesions[END_REF] jointly performs brain tissue classification and MS lesion segmentation by combining statistical and topological atlases. Finally, in [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF], the

Name

Definition

Unit Best Worse T N R T N F P + T N % 100 0

T P R T P T P + F N % 100 0

F P R F P F P + T N % 0 100 authors show that a global threshold on the FLAIR MR sequence, infered using an EM brain tissue classification, suffices to detect most MS lesions.

P P V T P T P + F P % 100 0 V O V ol(Seg ∩ GT ) V ol(Seg ∪ GT ) % 0 100 V D V ol(Seg) -V ol(GT ) V ol(GT ) % 0 < ∞ SD v∈∂(Seg) min u∈∂(GT ) d(u, v) + v∈∂(GT ) min u∈∂(Seg) d(u, v) card(Seg ∪ GT ) mm 0 < ∞ Table 1:
The final segmentation is then constrained to appear in the white matter by applying morphological operations.

The method proposed in [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF] won the MICCAI MS Segmentation Challenge 2008. For this specific method, the segmentation results on public and private datasets were made available by the authors and will be used as reference. In the rest of the article, methods will be identified by their reference.

Data preprocessing

We sub-sample and crop the images so that they all have the same size, 159 × 207 × 79 voxels, and the same resolution, 1 × 1 × 2 mm 3 . Sub-sampling and cropping intend to reduce the time spent learning the classifier. The preprocessing procedure corrects for RF acquisition field inhomogeneities [START_REF] Prima | Maximum likelihood estimation of the bias field in MR brain images: Investigating different modelings of the imaging process[END_REF] and performs inter-subject intensity calibration [START_REF] Rey | Détection et quantification de processus évolutifs dans des images médicales tridimensionnelles : application à la sclérose en plaques[END_REF]. Spatial normalization is also performed by aligning the mid-sagittal plane with the center of the images [START_REF] Prima | Computation of the mid-sagittal plane in 3d brain images[END_REF].

Spatial prior is added by registering the MNI atlas [START_REF] Evans | 3D statistical neuroanatomical models from 305 MRI volumes[END_REF] to the anatomical images, each voxel of the atlas providing the probability of belonging to WM, GM and CSF (cf. Figure 1). The Image Fusion module of MedINRIA (INRIA, 2010) is used to perform the affine registration of the MNI atlas onto every single case.

Notations

The multi-channel aspect of the method presented in this article requires to carefully define and name each channel. MR images from the MSGC dataset will be noted I s where the index s ∈ {T 1, T 2, F LAIR} stands for an MR sequence. Registered spatial priors will be noted P t where the index t ∈ {W M, GM, CSF } stands for a brain tissue class.

Although having different semantics, anatomical images and spatial priors will be treated under the unified term signal channel and denoted C ∈ {I T 1 , I T 2 , I F LAIR , P W M , P GM , P CSF }. The data consists of a collection of voxel samples described by their spatial position x = (x, y, z). Voxels can be evaluated in all available signal channels.

The value of the voxel x in channel C is denoted C(x).

Methods

This section describes our adaptation of the random decision forests to the segmentation of MS lesions and illustrates the visual features employed.

Context-rich decision forest

Our detection and segmentation problem can be formalized as a binary classification of voxel samples into either background or lesions. This classification problem is addressed by a supervised method: discriminative random decision forest, an ensemble classifier using decision trees as base classifiers.

Decision trees are discriminative classifiers which are known to suffer from over-fitting [START_REF] Breiman | Classification and Regression Trees[END_REF]. A random decision forest (Amit and Geman, 1997) achieves better generalization by growing an ensemble of many independent decision trees on a random subset of the training data and by randomizing the features made available to each node during training [START_REF] Breiman | Random forests[END_REF].

Forest training

The training data consists of a set of labeled voxels

T = {x k , Y (x k )}
where the label Y (x k ) is given by an expert. When asked to classify a new image, the classifier aims to assign every voxel x in the volume a label y(x).

In our case, y(x) ∈ {0, 1}, 1 for lesion and 0 for background.

The forest has T components with t indexing each tree. During training, all observations (voxels) x k are pushed through each of the trees. Each internal node applies a binary test [START_REF] Shotton | Textonboost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[END_REF][START_REF] Yi | Discriminative, semantic[END_REF][START_REF] Lempitsky | Random forest classification for automatic delineation of myocardium in real-time 3D echocardiography[END_REF][START_REF] Criminisi | Regression forests for efficient anatomy detection and localization in CT studies[END_REF] as follows:

t τ low ,τup,θ (x k ) =    true, if τ low ≤ θ(x k ) < τ up f alse, otherwise (1)
where θ is a function identifying the visual feature extracted at position During training, each node p is optimized using the partition of the training data T p it receives as input. At the end of the training process, each node p is assigned to the optimal binary test t λ * p , where λ * p = (τ * low , τ * up , θ * ) p . The optimality criterion is the information gain, denoted IG, as defined in [START_REF] Quinlan | C4.5: Programs for Machine Learning[END_REF])

IG(λ, T p ) = H(T p ) -H(T p |(t λ (x k )) k ) (2)
where T p ⊂ T and where H denotes the entropy. More precisely, the term

H(T p |(t λ (x k )) k
) measures the error made when approximating the expert labeling Y by the binary test t λ . The optimal parameter λ * p maximizes the information gain

λ * p = arg max λ IG(λ, T p ) (3)
for node p. As a result, the optimal binary test is the test discriminating lesion from background voxels such as maximizing the information gain.

Only a randomly sampled subset Θ of the feature space is available at each node for optimization, while the threshold space is uniformly discretized.

The optimal λ * = (τ * low , τ * up , θ * ) is found by exhaustive search jointly over the feature and threshold space. Random sampling of the features leads to increased inter-node and inter-tree variability which improves generalization [START_REF] Breiman | Random forests[END_REF].

Trees are grown to a maximum depth D. At the node level, a leaf node is generated when the information gain is below a minimal value IG min .

As a result of the training process, each leaf node l t of every tree t receives a partition T lt of the training data. The partition T lt can be divided into two sets respectively containing background and lesion voxels and defined as

T b lt = {(x, Y (x)) ∈ T lt |Y (x) = b} (4)
where b ∈ {0, 1} stands for the background and lesion class, respectively.

Subsequently, the following empirical posterior probability is defined

p lt (Y (x) = b) = |T b lt | |T lt | (5)
and stored at the leaf node.

Figure 2 illustrates how the decision trees partition the data in the feature space and how resulting probabilities are stored in leaf nodes.

Prediction

When applied to a new test volume

T test = {x k },
p(y(x) = b) = 1 T T t=1 p lt (Y (x) = b) (6)
which is a mean over all the trees in the forest. This probability may be thresholded at a fixed value τ posterior if a binary segmentation is required.

A posterior map P b is obtained by applying the same prediction procedure to all voxels. Thus for every voxel x, P b (x) = p(y(x) = b) is the posterior Larger opacities indicate larger probability of a pixel belonging to a class while uncertain regions are indicated by less saturated colours. The white line plots the locus of points for which p(y(x) = green) = p(y(x) = red). We observe that 1) forest posteriors mimic the maximum-margin behavior, 2) uncertainty increases when moving away from the training data.

Advantages

The probabilistic random forest framework presented in this section shows considerable advantages over other classifiers, e.g. support vector machines (SVMs). Indeed, it combines efficient probabilistic classification and transparent feature selection as detailed below.

When a new case is presented to the classifier, every voxel goes through a sequence of decisions on different channels. As a result, the posterior probability affected to each voxel measures the confidence of this voxel being an MS lesions in the multi-channel space infered from training.

Trees from the same forest are all independent from one another. The posterior map can thus be computed in parallel: each tree computes its own posterior map, they are then combined to form the final result (cf.

Section 3.1.2). The training of the random forest can also be parallelized in a trivial manner by learning each tree independently from the others.

Moreover, unlike more "black-box" supervised methods such as SVMs or neural networks, the random forest framework enables us to enter the learned trees and identify the most discriminative features. In Section 5.4, we take advantage of this property to draw a detailed analysis of the most discriminative visual features for the task of MS segmentation.

The following properties motivate the use of random forests framework for the task of MS lesions segmentation: 1) when applying the random forest the binary trees can be evaluated extremely fast in prediction; 2) thanks to parallelism training on MR volumes is fast; 3) as a result of the training process, an optimal sequence of decisions, including most informative channels and visual features, is effortlessly available (cf. Section 5.4.2).

In many applications, random forests have been found to generalize better than SVM or boosting (Yin et al., 2010). The generalization power increases monotonically with increasing the forest size. Unlike SVM and boosting, random forests also estimate the confidence of the prediction as a by-product of the training process. Figure 3 shows the feature space of an exemplary segmentation problem. It illustrates that the random forest also generalizes well in regions of the feature space with sparse data support, which is beneficial.

The uncertainty increases in those areas in feature space which have little data support and on the boundary between classes. This is an important expected behavior, similar to that of Gaussian Processes [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Another interesting property of random forests lies in the way they separate the feature space. In Figure 3, probability maps show the highest uncertainty values on voxels equidistant of the two classes. This is a feature the random forest classifier shares with maximum-margin classifier such as SVMs (Bernhard Schlkopf and Smola, 1999).

Visual features

This section aims at presenting the visual features and at arguing the underlying motivation.

Two kinds of visual features are computed:

1. local features:

θ loc C (x) = C(x) (7)
where C is an intensity or a prior channel, and C(x) is the value of channel C at position x;

2. context-rich features comparing the voxel of interest with distant regions.

The first context-rich feature compares the local voxel value in channel C 1 with the mean value in channel C 2 over two 3D boxes R 1 and R 2 within an extended neighborhood:

θ cont C 1 ,C 2 ,R 1 ,R 2 (x) = C 1 (x) - 1 V ol(R 1 ) x ′ ∈R 1 C 2 (x ′ ) - 1 V ol(R 2 ) x ′ ∈R 2 C 2 (x ′ ) (8)
where C 1 and C 2 are both intensity or prior channels. The regions R 1 and R 2 are sampled randomly in a large neighborhood of the voxel x (cf. Figure 4).

The sum over these regions is efficiently computed using integral volume processing [START_REF] Shotton | Textonboost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[END_REF]. The random sampling of the features is part of the random forest framework and were used in previous work [START_REF] Criminisi | Decision forests with longrange spatial context for organ localization in CT volumes[END_REF][START_REF] Criminisi | Regression forests for efficient anatomy detection and localization in CT studies[END_REF][START_REF] Yi | Discriminative, semantic[END_REF].

The second context-rich feature compares the voxel of interest at x with its symmetric counterpart with respect to the mid-sagittal plane, noted S(x):

θ sym C (x) = C(x) -C(S(x)) (9)
where C is an intensity channel.

A new version of the symmetry feature loosens the hard symmetric constrain as the size of the neighborhood increases, in order to take into account the fact that the brain is not perfectly symmetric. Instead of comparing with the exact symmetric S(x) of the voxel, the minimal difference value within a close neighborhood of S(x), denoted S, is chosen. Three different sizes are considered for this neighborhood: 6, 26 and 32 neighbors respectively (cf.

Figure 4). We obtain a softer version of the symmetric feature which reads:

θ sym C,S (x) = min x ′ ∈S {C(x) -C(x ′ )} (10)
To summarize, three visual features are introduced: local, neighborhood and symmetry features. They can be thought as meta-features to be used on top of standard image filters, such as local moments, gradients or textures, rather than replacing them. In our case, visual features are evaluated on raw images and spatial priors, but any other channel could be added, e.g.

intensity gradient as described in [START_REF] Yi | Discriminative, semantic[END_REF].

The presented features can be applied individually to every voxel, unlike e.g. geometric moments. This is an essential property which enables voxel-wise image classification. The presented method not only provides a voxel-wise classification, but also integrates neighborhood information in the classification process thanks to context-rich features. The use of 3D boxes to integrate neighborhood information is motivated by the fact that they are extremely efficient through integral volume processing [START_REF] Shotton | Textonboost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[END_REF].

Experiments and results

Results presented in this section aim at evaluating the segmentation results and comparing the context-rich random forest approach to methods presented during the challenge [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF][START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF]Bricq et al., 2008a;[START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF]. Experiments described here are discussed in Section 5.

Exhaustive segmentation results are available for both public and private datasets under the following url:

Results on the public MSGC dataset

For quantitative evaluation, the 20 available cases from the public dataset are classified and compared to other methods [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF][START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF]Bricq et al., 2008a;[START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF]. A three-fold crossvalidation is carried out on this dataset: the forest is trained on 2 3 of the cases and tested on the other 1 3 , this operation is repeated 3 times in order to collect test errors for each case.

The binary classification is evaluated using two measures, true positive rate (T P R) and positive predictive value (P P V ), both equal 1 for perfect segmentation (cf. Table 1). Tables in supplemental material report extensive results allowing comparison on every case of the MSGC public dataset. It shows that the learned context-rich random forest achieves better T P R in all cases (cf. top bar plot), and better P P V in 70% of the cases (cf. center bar plot). Computed p-values for the pair-sample t-test show that these improvements are significative for both T P R (p = 1.3 • 10 -7 ) and P P V (p = 0.0041) scores. [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF] on the public dataset. Relative improvement over [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF], defined as RI = (score RF -score other )/score other , are significant for both T P R (p = 1.3 • 10 -7 ) and P P V (p = 0.0041) scores. Significant improvements over [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF] are highlighted in bold. The MSCG website carried out a complementary and independent evaluation of our algorithm on the previously unseen private dataset. The results, reported in Table 3, confirm a significant improvement over [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. The presented spatial random forest achieves, on average, slightly larger true positive (T P R), which is beneficial (cf. Table 1), and comparable false positive (F P R) rates but lower volume difference (V D) and surface distance (SD) values (cf. Table 3). Pair-sample p-values were computed for the t-test on the private dataset. Results show significant improvement over the method presented in [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF] on SD (p = 4.2 • 10 -6 ) for the CHB rater, and on SD (p = 6.1 • 10 -3 ) for the UNC rater.

Results on the private MSGC dataset

Discussion

Interpreting segmentation results

Quantitative evaluation of segmentation results, for both public (cf. Tables 2) and private (cf. Tables 3,4, 5 and 6) datasets, show that the presented random forest framework compares favorably to top-ranked methods. More specifically, they show significant improvement on the algorithm presented in [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. Exhaustive results are available in the supplemental material allow case-by-case comparison with other methods.

The MSGC website [START_REF] Styner | MS lesion segmentation challenge[END_REF] Independent quantitative evaluation confirms improvement on the algorithm presented in [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. Boldface highlights significant improvements. The spatial random forest achieves, on average, slightly larger true positive (T P R), which is beneficial and comparable false positive (F P R) rates but lower volume difference (V D) and surface distance (SD) values.

rater. The method presented in (Bricq et al., 2008a) and our context-rich random forest approach are first and second respectively. They show very close scores, 82.1354 and 82.0755 respectively. A score of 82 places the accuracy of the method just below that of a human expert. In addition, we know that the reliability of automatic methods is generally higher than that of human experts. Approaching the performance of a human rater is thus even more interesting. [START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF]. The relative mean improvement on the algorithm presented [START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF] on the private dataset is defined as follows RI = (score RF -score Anbeek )/score Anbeek . The RI and the p-values a and is reported on top for each metric of the MSGC, associated p-values are reported below.

Independent quantitative evaluation confirms improvement on the algorithm presented in [START_REF] Anbeek | Automated MS-lesion segmentation by K-Nearest neighbor classification[END_REF]. Boldface highlights significant improvements. The spatial random forest achieves, on average, better results on UNC than on CHB labels: higher true positive (T P R), which is beneficial and lower false positive (F P R) rates, lower volume difference (V D) and surface distance (SD) values.

Although segmentation results include most MS lesions delineated by the expert (cf. Figures 5 and6), we observe that some MS lesions are missing.

Missed MS lesions are located in specific locations which are not represented in the training data, e.g. in the corpus callosum (cf. Figure 5, slice 38).

This is a limitation of the supervised approach. In this very case, however, the posterior map highlights the missed lesion in the corpus callosum as belonging to the lesion class with high uncertainty. Low confidence (or high Table 5: Average results computed by the MSGC on the private dataset and compared to the method presented in (Bricq et al., 2008a). The relative mean improvement on the algorithm presented (Bricq et al., 2008a) on the private dataset is defined as follows RI = (score RF -score Bricq )/score Bricq . The RI and the p-values a and is reported on top for each metric of the MSGC, associated p-values are reported below.

Independent quantitative evaluation confirms improvement on the algorithm presented in (Bricq et al., 2008a). Boldface highlights significant improvements. The spatial random forest achieves, on average, slightly larger true positive (T P R), which is beneficial, but also slightly larger false positive (F P R) rates, and lower volume difference (V D) and surface distance (SD) values.

uncertainty) reflects the incorrect spatial prior inferred from an incomplete training set. Indeed, in the training set, there is no example of MS lesions appearing in the corpus callosum.

On the contrary, the random forest is able to detect suspicious regions [START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF]. The relative mean improvement on the algorithm presented [START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF] on the private dataset is defined as follows RI = (score RF -score Shiee )/score Shiee . The RI and the p-values a and is reported on top for each metric of the MSGC, associated p-values are reported below.

with
Independent quantitative evaluation confirms improvement on the algorithm presented in [START_REF] Shiee | Multiple sclerosis lesion segmentation using statistical and topological atlases[END_REF]. Boldface highlights significant improvements. The spatial random forest achieves, on average, slightly larger true positive (T P R), which is beneficial but slightly larger false positive (F P R) rates, and lower volume difference (V D) and surface distance (SD) values.

appearance model and spatial prior implicitly learned from the training data points out that hyper-intense regions in the FLAIR MR sequence which lay in the white matter (cf. Section 5.4) can be considered as MS lesions with high confidence.

Recent histopathological studies have shown that grey matter regions are also heavily affected by the MS disease [START_REF] Geurts | Grey matter pathology in multiple sclerosis[END_REF]. In our case, the public dataset does not show any MS lesion in the grey matter of the brain. Subsequently, the decision forest learns that MS lesions preferably appear in the white matter. Adding new cases showing grey matter MS lesions in the training set would allow the forest to automatically adapt the segmentation to include this kind of lesions. This observation stresses the necessity of gathering large and heterogeneous datasets for training purposes.

When focusing on quantitative measures, we observe that cases UNC01 and UNC06 from the public dataset show surprisingly low scores (cf. Table 2). The labels by the CHB expert for these two cases are abnormal: the ground truth is mirrored with respect to the anatomical images. This may be considered as a label error and explains the low scores for these two specific cases. The MSGC website confirmed this observation and subsequently corrected the online database. We also observe that learning on the whole public dataset and testing on the private dataset (cf. Table 3) produces better average results than the three-fold cross-validation carried out on the public dataset (cf. Table 2). Again this illustrates the benefit of learning the classifier on large enough datasets capturing better the variability of the data.

Influence of preprocessing

Data normalization is critical in order to ensure that the feature are evaluated in a coherent way in all the images presented to the forest. The evaluation of context-rich features, θ cont , is sensitive to rotation: spatial normalization is performed using rigid registration [START_REF] Prima | Maximum likelihood estimation of the bias field in MR brain images: Investigating different modelings of the imaging process[END_REF]. In the same way, the evaluation of intensity-based features requires inter-case intensity calibration.

Classification results for cases from the CHB (cf. Figures 5 and6) and the UNC (cf. Figures 7 and8) centers are obtained with the same forest. It is mandatory to use the same preprocessing as during training (cf. Section 2.2).

By doing so, the cases from different datasets, e.g. T1-weighted and FLAIR images in Figures 5 and6, show very similar intensity values for a specific brain tissue and a given MR sequence. However, we observe that the contrast in the T1-weighted and FLAIR images is more marked in case CHB05 (cf. Figure 5) than in case UNC02 (cf. Figure 6). Despite contrast changes, classification results are coherent. This illustrates the stability of our method, the random forest framework together with its preprocessing step to slight inter-image contrast variations.

The trees are generated in parallel on 30 nodes and gathered to form the forest. Cropping and sub-sampling the training images aims at reducing, by a factor larger than 10, the execution time spent to learn a single tree. On IBM e325 dual-Opterons 246 at a maximum frequency of 2Ghz, learning a tree on 20 sub-sampled images and with parameters fixed in Section 4.2 on a single CPU takes, on average, 8 hours.

Influence of forest parameters

The number of the trees and their depth, respectively denoted T and D, characterize the generalization power and the complexity of the nonparametric model learned by the forest. This section aims at understanding the contribution of each of these meta-parameters.

A 3-fold cross-validation on the public dataset is carried out for each parameter combination. Segmentation results are evaluated for each combination using two different metrics: the area under the receiver operating characteristic (ROC) curve and the area under the precision-recall curve.

The ROC curve plots T P R vs. F P R scores computed on the test data for every value of τ posterior ∈ [0, 1]. The precision-recall curve plots P P V vs.

T P R scores computed on the test data for every value of τ posterior ∈ [0, 1].

Results are reported in Figure 7.

We observe that 1) for a fixed depth, increasing the number of trees leads to better generalization; 2) for a fixed number of trees, low depth values lead to underfitting while high depth values lead to overfitting; 3) overfitting is reduced by increasing the number of trees.

This analysis was carried out a posteriori. Forest parameters were indeed selected in a safety-area with respect to under-and overfitting. The safety-area corresponds to a sufficiently flat region in the evolution of the areas under the ROC and the precision-recall curve. As shown in Figure 8, increasing the number of trees tends to benefit the generalization power of the classifier. We also observe that the performance of the classifier stabilizes for large enough forests.

Analysis of feature relevance

During training, features considered for node optimization form a large and heterogeneous set (cf. Section 3.2). Unlike other classifiers, random forests provide an elegant way of ranking these features according to their the ROC curve. We observe that, for a fixed depth D = 14, increasing the number of trees improves generalization as stated in [START_REF] Breiman | Random forests[END_REF]. The increase in performance stabilizes around the value T = 30.

discriminative power. In this section, we aim at better understanding which are the most discriminative channels and visual cues (local, context-rich or symmetric) used in the classification process.

Most discriminative visual features

The first approach consists in counting the nodes in which a given feature was selected. We observe that local features were selected in 24% of the nodes, context-rich features were selected in 71% of the nodes whereas symmetry features were selected in 5% of the nodes (cf. Figure 9). In this case, no distinction is made as for the depth at which a given feature was selected.

Context-rich features exhibit high variability (900 of them are randomly sampled at every node). This variability combined with their ability to highlight regions which differ from their neighborhood explains why they were chosen. Together with local features, context-rich features learn a multichannel appearance model conditioned by tissue spatial priors. Symmetry features are under-represented in the forest and thus prove to be the least discriminative ones. This is due to the fact that a large proportion of periventricular MS lesions tend to develop in a symmetric way. Nevertheless, symmetric features appear in top levels of the tree (up to third level) which indicates that they provide an alternative to local and context-rich features when these two fail.

A finer estimation of the feature importance consists in weighting the counting process. For a given feature, instead of only counting the nodes in which it appears, we also take into account the proportion of lesion voxels it helps discriminating: "the larger the proportion of lesion voxels it helps to discriminate, the larger the weight of the feature". This leads us to define for a fixed depth value d, the importance of a given feature type, denoted (IF T ), as:

IF T (α) = 1 |T 1 | p |T 1 p | • χ α (θ p ) α ∈ {loc, cont, sym} (11)
where α denotes a feature type, p indices the nodes in layer d, T 1 is the training set of lesion voxels which partition T 1 p reached node p, and χ is the indicator function such that:

χ α (θ p ) =    1, if θ p is of type α 0, otherwise (12)
The feature importance evaluates to 21.1 % for local features, 76.6 % for context-rich features and 2.3 % for symmetry features. Results are comparable to those obtained only by counting the features in the forest, but the real advantage of this measure is to allow us to draw depth-by-depth feature importance analysis in a normalized way.

The feature importance as a function of the depth of the tree is reported in Figure 10. Presented results are averaged values over a forest containing T = 30 trees. Again, we observe that context-rich features are predominantly selected as the most discriminative, which confirms the trend reported in Figure 9. However, as shown in Figure 10, the preponderance of context-rich features is not uniform throughout the tree. Indeed, local features are the most discriminative in layers 0 and 2. A careful analysis of selected channels helps understanding why local features are selected in the top layers of the tree (cf. Section 5.4.2).

The selected context-rich features show high variability. More specifically, the long-range regions are distributed all over the neighborhood. Depthby-depth analysis does not show any specific pattern in the position of the regions with respect to the origin voxel. In addition, the volume of the regions also show high variability. The observed heterogeneity of selected contextrich features aims at coping with the variability of MS lesions (shape, location and size).

The symmetry feature is under-represented in the forest. Its discriminative power is thus very low compared to local and context-rich features.

This observation induces two complementary interpretations to explain why symmetry features are the least significant: 1) most of MS lesions appear in peri-ventricular regions and in a symmetric way, 2) most of MS lesions can be clearly identified by their signature across MR sequences and their relative position in the white matter of the brain. However, in deeper layers of the trees, the symmetry feature is more significant and tends to classify ambiguous asymmetrical regions. When looking into the selected features, we also notice that the hard symmetric constraint is preferred over the loose symmetric constrain (cf. Section 3.2). Indeed, the feature importance evaluates to 1.6% for the hard symmetric feature, and to 0.7% for the loose symmetric feature. Moreover, in the rare cases where the loose constrain is selected, the 6-neighbors version predominates (cf. Section 3.2). This observation supports the idea that considering brain hemispheres as symmetric is an accurate approximation in our specific setting (cf. Sections 2 and 3).

Most discriminative channels

The second approach focuses on the depth at which a given feature was selected. For every tree in the forest, the root node always applies a test on the FLAIR sequence (θ loc F LAIR ). It means that out of all available features, containing local, context-rich and symmetry multi-channel features, θ loc F LAIR was found to be the most discriminative. At the second level of the tree, a context-rich feature on spatial priors (θ cont W M,GM ) appears to be the most discriminative over all trees in the forest. It aims at discarding all voxels which do not belong to the white matter.

The optimal decision sequence found while training the context-rich forest can thus be thought as a threshold on the FLAIR MR sequence followed by an intersection with the white matter mask (cf. Figure 11). Interestingly, this sequence matches the first and second steps of the pipeline proposed by the challenge winner method [START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]. Note that in our case, it is automatically generated during the training process. Deeper layers in the trees, then, refine the segmentation of MS lesions by applying more accurate decisions. The feature importance (cf. Equation 11) can be extended in a straightforward way to be parametrized not only by the type of feature (local, context-rich, symmetric) but also by the channel. When globally looking at the selected channels (cf. Figure 12), we notice that their importance varies throughout the tree: first layers, as mentioned before, favor detection of bright spots in the white matter by successively testing the FLAIR MR sequence, spatial priors on WM and GM tissues and finally testing on the T2 MR sequence; deeper layers take into account other modalities to adjust the segmentation.

Conclusion

We demonstrated the power of the RF formulation applied to the difficult task of MS lesion segmentation in multi-channel MR images. We presented the first layer performs a local test on the FLAIR MR sequence in order to detect bright spots, the second one discards all voxels which do not belong to the white matter by using context-rich information over the WM and GM channels. Note that a large spectrum of available channels is tested throughout the tree. three kinds of 3D features based on multi-channel intensity, prior and contextrich information. Those features are part of a context-rich random decision forest classifier which demonstrated improved results on one of the state of the art algorithms on the public MS challenge dataset. In addition, the random decision forest framework provided a means to automatically select the most discriminative features to achieve the best possible segmentation.

Future work could include the use of more sophisticated features to reduce even further the preprocessing requirements. The context-rich random forest framework presented in this article is generic which is an additional strength of the method. It can be applied as is to any other segmentation task, e.g. brain tumors segmentation in multi-sequence MR images of the brain.

Finally, one could investigate an extension of the proposed approach to larger multi-class problems in order to try to simultaneously segment brain tissues (WM, GM, CSF) along with MS lesions.

  The evaluation metrics true negative rate (T N R), true positive rate (T P R), false positive rate (F P R) and positive predictive value (P P V ) are defined using the following notations: true positives (T P ), true negatives (T N ), false positives (F P ) and false negatives (F N ). The volume overlap (V O) and the relative absolute volume difference (V D) evaluates the differences between the segmentation (Seg) and the ground truth (GT ) by computing their volume (V ol). The average symmetric surface distance (SD) measures how close the segmentation and the ground truth are from each other using the Euclidean distance d on the set of boundary voxels noted ∂. The best, respectively worse, column contains the metric score of the perfect segmentation, respectively of a completely-off segmentation.

Figure 1 :

 1 Figure 1: Case CHB07 from the public MSGC dataset. From top to bottom: three axial slices of the same patient. From left to right: preprocessed T1-weighted (I T 1 ), T2weighted (I T 2 ) and FLAIR MR images (I F LAIR ), the associated ground truth GT and the registered white matter atlas (P W M ).

x

  k . There are several ways of defining θ, either as a local intensity-based average, local spatial prior or context-rich cue. These are investigated in more detail in the next section. The value of the extracted visual feature is thresholded by τ low and τ up . The voxel x k is then sent to one of the two child nodes based on the outcome of this test.

  each voxel x k is propagated through all the trees by successive application of the relevant binary tests. When reaching the leaf node l t in all trees t ∈ [1..T ], posteriors p lt (Y (x) = b) are gathered in order to compute the final posterior probability defined as follows:

Figure 2 :

 2 Figure 2: Decision trees encode feature space partitions. (a) A decision tree of depth D = 2 is considered in this example. Decision node 1 and leaf nodes 5 and 6 are colored to track the partitions of the training data in the feature space. The black cross stands for an unseen sample (voxel) which is classified while propagated down the tree. (b) A zoom on node 2 shows that its binary test, denoted by t τ low,2 ,τup,2,θ2 , is optimized over a partition of the training data, denoted by T 2 = {v k,2 , Y (v k,2 )}. The leaf node 6 encloses the class distribution of the set of voxels reaching it during training. The classes are background in blue and lesion in red. (c) The dots stand for the training voxels and are colored according to their class. The black cross denotes a voxel from an unseen volume considered for prediction. Every decision node in the forest applies an axis-aligned feature test.Here we focus on decision nodes 0 and 2 using features θ 0 and θ 2 , respectively.

Figure 3 :

 3 Figure 3: Posterior maps learned from two distinct synthetic training sets. In both cases, the training data consists of two classes, green and red, and is used to learn a large forest, here T = 350. The posterior map is obtained by classifying a dense grid in the feature space and is then overlayed with the associated training data (shown as points).

Figure 4 :

 4 Figure 4: 2D view of context-rich features. (a) A context-rich feature depicting two regions R 1 and R 2 with constant offset relatively to x. (b-d) Three examples of randomly sampled features in an extended neighborhood. (e) The symmetric feature with respect to the mid-sagittal plane. (f) The hard symmetric constraint. (g-i) The soft symmetry feature considering neighboring voxels in a sphere of increasing radius. See text for details.

Forest

  parameters are fixed to the following values: number of random regions |Θ| ≃ 950, number of trees T = 30, tree depth D = 20, lower bound for the information gain IG min = 10 -5 , posterior threshold τ posterior = 0.5. Parameters T and D are set here to maximum values, Section 5.3 explains how these parameters can be optimized in order to improve segmentation results.

A

  context-rich random forest was learned on the whole public dataset from the MS Lesion Challenge, i.e. 20 labeled cases. Forest parameters are fixed to the following values: number of random regions |Θ| ≃ 950, number of trees T = 30, tree depth D = 20, lower bound for the information gain IG min = 10 -5 , posterior threshold τ posterior = 0.5. Considerations that lead to these parameter values are detailed in Section 5.3.

Figure 5 :Figure 6 :

 56 Figure 5: Segmenting Case CHB05 from the public MSGC dataset. From left to right: preprocessed T1-weighted (I T 1 ), T2-weighted (I T 2 ) and FLAIR MR images (I F LAIR ) overlayed with the associated ground truth GT , the posterior map P osterior = (P lesion (v k )) k displayed using an inverted grey scale and the FLAIR sequence overlayed with the segmentation (Seg = (P osterior τ posterior ) with τ posterior = 0.5). Segmentation results show that most of lesions are detected. Although some lesions are not detected, e.g. peri-ventricular lesion in slice 38, they appear enhanced in the posterior map. Moreover the segmentations of slices 38 and 42 show peri-ventricular regions, visually very similar to MS lesions, but not delineated in the ground truth.

Figure 7 :

 7 Figure 7: Influence of forest parameters on segmentation results. Both curves were plotted using mean results from a 3-fold cross validation on the public dataset. Left: the figure shows the influence of forest parameters on the area under the precision-recall curve. Right: the figure shows the influence of forest parameters on the area under the ROC curve.The ideal classifier would ensure area under the curve to be equal to 1 for both curves. We observe that 1) for a fixed depth, increasing the number of trees leads to better generalization; 2) for a fixed number of trees, low depth values lead to underfitting while high values lead to overfitting; 3) overfitting vanishes by increasing the number of trees.

Figure 8 :

 8 Figure 8: Influence of the number of trees on segmentation results. Both curves were plotted using mean results from a 3-fold cross validation on the public dataset. Top: the figure shows the influence of the number of trees on the area under the precision-recall curve. Bottom: the figure shows the influence of the number of trees on the area under

Figure 9 :

 9 Figure 9: Ranking features according to the proportion of nodes in which they appear. Context-rich features are selected in 71% of the nodes, local features are selected in 24% of the nodes whereas symmetry features are selected in 5% of the nodes.

Figure 10 :

 10 Figure 10: Type of feature selected by layer of the tree. For a fixed depth, the red circle stands for the importance of the context-rich feature (θ cont ), while the green circle stands for the importance of the local feature (θ loc ). For clarity, symmetry features (θ sym ) are omitted as they are under-represented in the forest. The blue line monitors the proportion of training samples of the lesion class which do not reside in leaf nodes, foreach layer of the tree. We observe that context-rich features are predominantly selected as the most discriminative ones except in layers 0 and 2.

Figure 11 :

 11 Figure 11: Combination of features and channels learned by the forest to discriminate MS lesions. The first layer of all trees in the forest performs a threshold on the FLAIR MR sequence. The second one discards all voxels which do not belong to the white matter. The posterior map is obtained by using a forest with trees of depth 2 and thus highlights hyper-intense FLAIR voxels which lie in peri-ventricular regions.

Figure 12 :

 12 Figure 12: Channel importance as a function of the depth of the tree. Plots draw the channel importance drawn as a function of the depth of the tree for both local (top) and long-range features (bottom). For a fixed depth, only the most discriminative channel is depicted. Note how successive layers of the tree test complementary channels:

Table 2 :

 2 Metric [%] Souplet et al. Context-rich RF RI [%] p-value T P R 19.21 ± 13.68 39.39 ± 18.40 105 1.3 • 10 -7 Comparison of context-rich random forests with method presented in

	P P V	29.55 ± 16.26 39.78 ± 20.19	35	0.0041

Table 3 :

 3 Average results computed by the MSGC on the private dataset and

	gathers the results of the meth-

compared to the method presented in

[START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]

. The relative mean improvement on the algorithm presented in

[START_REF] Souplet | An automatic segmentation of T2-FLAIR multiple sclerosis lesions[END_REF]

on the private dataset is defined as follows RI = (score RF -score Souplet )/score Souplet . The RI and the p-values a and is reported on top for each metric of the MSGC, associated p-values are reported below.

Table 4 :

 4 Average results computed by the MSGC on the private dataset and compared to the method presented in

	Rater Metric [%] Anbeek et al. Context-rich RF RI [%] p-value
		V D	46.93 ± 50.41 52.94 ± 28.63 +12.8	0.62
	CHB	SD T P R 59.14 ± 21.79 58.08 ± 20.03 -1.80 7.85 ± 11.00 5.27 ± 9.54 -32.8 5.9 • 10 -3 0.84
		F P R 78.51 ± 20.56 70.01 ± 16.32 -10.8 8.6 • 10 -2
		V D	100.7 ± 132.3 50.56 ± 41.41 -49.8 9.1 • 10 -2
	UNC	SD T P R 48.86 ± 21.64 51.35 ± 19.98 +5.10 9.77 ± 9.00 5.6 ± 6.67 -42.7 3.6 • 10 -4 0.38
		F P R 83.19 ± 18.59 76.81 ± 11.70 -7.68	0.11

  Rater Metric [%] Bricq et al. Context-rich RF RI [%] p-value

		V D	73.03 ± 78.80 52.94 ± 28.63 -27.5	0.19
	CHB	SD T P R 46.70 ± 19.94 58.08 ± 20.03 +24.4 1.0 • 10 -3 6.65 ± 6.55 5.27 ± 9.54 -20.8 0.20
		F P R	51.06 ± 25.23 70.01 ± 16.32 +37.1 8.3 • 10 -5
		V D	51.33 ± 27.00 50.56 ± 41.41 -1.50	0.92
	UNC	SD T P R 39.50 ± 16.06 51.35 ± 19.98 +30.0 1.2 • 10 -3 6.61 ± 5.23 5.6 ± 6.67 -15.2 0.17
		F P R	60.80 ± 22.75 76.81 ± 11.70 +26.3 8.6 • 10 -4

Table 6 :

 6 high certainty. Suspicious regions are visually very similar to MS lesions and widely represented in the training data, but they are not delineated by the expert, e.g. the left frontal lobe lesion again in Figure 5, slice 38. The Rater Metric [%] Shiee et al. Context-rich RF RI [%] p-value Average results computed by the MSGC on the private dataset and compared to the method presented in

		V D	84.17 ± 120.8 52.94 ± 28.63 -37.1	0.22
	CHB	SD T P R 55.40 ± 23.60 58.08 ± 20.03 +4.83 7.95 ± 16.65 5.27 ± 9.54 -33.5 9.0 • 10 -3 0.64
		F P R 68.85 ± 23.75 70.01 ± 16.32 +1.69	0.76
		V D	69.63 ± 115.2 50.56 ± 41.41 -27.4	0.42
	UNC	SD T P R 49.79 ± 24.54 51.35 ± 19.98 +3.14 7.104 ± 8.93 5.6 ± 6.67 -21.1	0.15 0.76
		F P R 74.28 ± 20.08 76.81 ± 11.70 +3.41	0.50

  Tuning the meta-parameters of the forest on the training data is not a valid practice. Using out-of-bag samples for forest parametrization is indeed preferable. Due to the fact that little training data is available for the MS lesion class, available labeled data was exclusively used to train the forest. From this perspective, the forest parameters were set to arbitrary but high enough values to avoid under-and overfitting: T = 30 and D = 20.


Edna Graneli
email: edna.graneli@lnu.se

Johannes A Hagström

Maria O P Moreira

Edna Granéli

Clarisse Odebrecht

Domoic acid production and elemental composition of two Pseudo-nitzschia multiseries 1 strains, from the NW and SW Atlantic Ocean, growing in phosphorus-or nitrogen-2 limited chemostat cultures. 3 4

Keywords: Pseudo-nitzschia, domoic acid, chemical composition, nutrient limitation, toxin Pseudo-nitzschia, Domoic acid, Chemical composition, Toxin, Nutrient limitation 18 19

come

Domoic acid production and elemental composition of two Pseudo-nitzschia multiseries strains, from the NW and SW Atlantic Ocean, growing in phosphorus-or nitrogen-limited chemostat cultures.

and Brazil (BR). The algae were grown as chemostat cultures at 0.2, 0.3 and 0.4 d -1 under 5 nitrogen (N) and phosphorus (P) deficient conditions. The level of deficiency significantly 6 affected the atomic C:N, C:P, C:Si and N:P ratios in both strains. In both strains P per cell 7 was 2-4 X higher in the N-than in the P-deficient cultures. The opposite was not found for N 8 in the P-deficient cultures, as shown by the N:P ratios and C:N ratios. The C:N and C:P ratios 9

were significantly lower in the CA strain, and this did not change due to level of deficiency. 10

The concentration and production of DA per cell per day, were significantly higher for both 11 strains under P-deficiency as expected since the toxin is rich in N. However, DA was also 12 produced by both strains during continuous cell division under N-deficiency. High or low 13 bacterial densities associated with P. multiseries did not increase or decrease DA production. 14

Our data imply that more attention needs to be given to the N:P ratios and concentrations in 15 the waters where these algae occur, as both N-and P-deficiency affect DA production and 16 cellular DA concentrations. 17 Mass mortalities of marine birds and mammals occur yearly in many coastal areas of 3 the world during and after blooms of the diatom Pseudo-nitzschia spp., especially during 4 strong upwelling. Several Pseudo-nitzschia spp. produce the responsible neurotoxin domoic 5 acid (DA), which causes Amnesic Shellfish Poisoning (ASP) in humans (cf. Bates and 6 Trainer, 2006). The most severe incident to date, from a human point of view, occurred in 7 1987 off Prince Edward Island (PEI), SE Canada. DA from P. multiseries Hasle 8 (Bacillariophyceae) bio-accumulated in blue mussels (Mytilus edulis Linnaeus). More than 9 100 humans consuming the mussels suffered from e.g. nausea and amnesia, and 3 people died 10 (cf. Landsberg, 2002). In the Galician Rias (NW Spain), DA originating from blooms of 11 Pseudo-nitzschia spp. was first detected in 1994 [START_REF] Míguez | First detection of domoic acid in Galicia 17 (NW Spain)[END_REF]. Pseudo-nitzschia spp. 12 and Nitzschia navis-varingica from brackish waters in Vietnam, Japan and the Philippines, 13 and from marine or brackish waters in Ireland, Denmark, Spain and Sweden are also known 14 to produce DA [START_REF] Lundholm | Studies on the marine 11 planktonic diatom Pseudo-nitzschia. 2. Autecology of P. pseudodelicatissima based on 12 isolates from Danish coastal waters[END_REF][START_REF] Cusack | Confirmation 8 of domoic acid production by Pseudo-nitzschia australis (bacillariophyceae) isolated 9 from Irish waters[END_REF][START_REF] Kotaki | Wide distribution of Nitzschia navis-varingica, 27 a new domoic acid-producing benthic diatom found in Vietnam[END_REF]. This shows 15 that a possible ASP outbreak with devastating effects on both ecosystems and economies is a 16 worldwide problem. 17

Published research has shown that a Pseudo-nitzschia species may be toxic in one part 18 of the world but not in another [START_REF] Bates | Bloom dynamics and physiology of 26 domoic-acid-producing Pseudo-nitzschia species[END_REF]. Studies on the influence of nutrient 19 availability on DA production by P. multiseries have so far been focused on phosphorus (P) [START_REF] Fehling | Growth dynamics of non-toxic 21 Pseudonitzschia delicatissima and toxic P. seriata (Bacillariophyceae) under simulated 22 spring and summer photoperiods[END_REF] and silicon (Si) deficiencies [START_REF] Bates | Bloom dynamics and physiology of 26 domoic-acid-producing Pseudo-nitzschia species[END_REF]. These studies have most of the time been 21

performed in batch cultures. This was for instance done when investigating the effect of 22 various nutrient deficiencies, including nitrogen (N), and when comparing toxin production in 23 different strains of the same Pseudo-nitzschia species [START_REF] Bates | Controls on domoic acid production by the diatom Nitzschia pungens 23 f. multiseries in culture: nutrients and irradiance[END_REF]Fehling et al., 24 2004). Continuous cultures on the other hand allow a dynamic equilibrium between nutrient 25 input and growth, thereby simulating natural conditions in the lab better than most other 26 methods. Continuous culturing has been successfully used to study e.g. the effect of P-27 deficiency on P. multiseries (Pan et al., 1996b) Since DA is an amino acid, i.e. N is an 28 essential element, it has been argued that some strains stop producing the toxin during N 29 deficiency [START_REF] Bates | Ecophysiology and metabolism of ASP toxin production[END_REF] different rates in a chemostat (0.2 -1.2 d -1 ). However, to our knowledge, there is a lack of 1 knowledge regarding the possible differences between strains within the same species isolated 2 from different locations and grown under identical continuous culture conditions with 3 different N:P ratios. 4

The questions we focused on during this experiment were therefore: (1) Are there any 5 physiological differences between two strains of P. multiseries from different origins (SW 6 and NW Atlantic Ocean) concerning DA production when growing in N-and P-deficient 7 conditions? (2) What is the effect of growth rate (i.e. level of nutrient deficiency) on cell 8 physiology and DA production for the two strains? The experiment was carried out using two non-axenic strains of the toxic diatom 15

Pseudo-nitzschia multiseries Hasle. One strain originating from the Provasoli- Guillard 16 National Center for Culture of Marine Phytoplankton (CCMP 1660), USA and had been 17 isolated from waters of eastern Prince Edward Island (Canadian SE coast) in December 1993. 18

The other strain was isolated in June 2001, from the southern part of the Atlantic Ocean 19 (Brazil) and was obtained from the Phytoplankton Laboratory of the Fundação Universidade 20

Federal do Rio Grande (FURG), Brazil. Both strains had been identified as P. multiseries 21 using electron microscopy. In addition, an 18S rRNA gene sequence analysis of both strains 22 was done prior to the experiment in order to verify that the two strains indeed are the same 23 species (Fig. 1). 24

The partial 18S rRNA gene was amplified from Pseudo-nitzschia spp. cultures by 25 adding 5 µL of cells resuspended in TE buffer (10 mM Tris, 0.5 mM EDTA, pH 8.0) to the 26 polymerase-chain-reaction (PCR) mix. The PCR mix contained 200 mM deoxy-nucleoside 27 triphosphate (dNTP), 5% dimethylsulfoxide, 5 pmol each of primers EUK1 and 518r, Taq 28 DNA polymerase reaction buffer with 1.5 mM MgCl2 and 1 U of a mixture of Taq DNA 29 polymerase (Roche) in a total reaction volume of 50 µL [START_REF] Janson | Evidence for multiple 19 species within the endoparasitic dinoflagellate Amoebophrya ceratii as based on 18S 20 rRNA gene-sequence analysis[END_REF]. The 30 amplification was performed with a GeneAmp PCR System 2400 (Perkin Elmer, USA) with 31 an initial incubation of 5 min at 94°C. The 30 subsequent cycles were 30 seconds at 94°C, 30 32 seconds at 50°C and 1 min at 72°C. A final incubation of 5 min at 72°C finished the PCR 33 amplification. The gel electrophoresis of the PCR products was performed according to The PCR products were later sequenced on both strands after purification with a 2

QIAquickTM PCR kit (Qiagen) using the primers EUK1 and 518r. The sequencing reaction 3 was performed using dye terminator chemicals (BigDye: Perkin Elmer) according to the 4 instructions provided by the manufacturer. Sequences were read using a DNA auto sequencer 5

Model ABI 377-XL (Perkin Elmer). The sequences reported here, and those retrieved from 6 GenBank, were aligned using CLUSTAL W 1.8 [START_REF] Thompson | CLUSTAL W -improving the 19 sensitivity of progressive multiple sequence alignment through sequence weighting, 20 position-specific gap penalties and weight matrix choice[END_REF], with the gap and 7 extension penalty set to 6.6. Phylogenetic analyses were performed on the PHYLO_WIN 8 software [START_REF] Galtier | SEAVIEW and PHYLO_WIN: Two graphic 32 tools for sequence alignment and molecular phylogeny[END_REF]. The tree was edited using TREEVIEW 1.5 for PC (Page, 9 1996). The Jukes and Cantor method was used to estimate the distances (Jukes and Cantor, 10 1969). 11 12

Culture maintenance 13 14

The two strains were maintained in f/2 medium, with the trace metal solution modified 15 according to medium L1 [START_REF] Guillard | Culture methods[END_REF]. However, the nitrogen concentration was reduced 16 so that the 16N:1P Redfield ratio was attained (see below). The cultures were grown under a 17 16:8 hours light:dark cycle. Photosynthetically Active Radiation (PAR) was adjusted to 200 18 µmol photons m -2 s -1 (Philips TLD 36W/830 cool white fluorescent tubes) and the 19 temperature was kept at 18 ± 0.5°C. This was followed by sub-culturing for one month in f/10 20 medium, with the same light, temperature and adjustments of nutrient ratio as above. These 21 sub-cultures were later used as inocula for the chemostat cultures. All media were prepared in 22 10 L Nalgene polycarbonate carboys using filtered (Munktel glass fibre filter, 1.2 µm pore 23 size) aged Baltic Sea surface water (6 psu) adjusted to 31 ± 1 psu. Proportional quantities of 24 the major salts (NaCl, MgCl 2 6H 2 O, CaCl 2 2H 2 O, Na 2 SO 4 , NaHCO 3 and K 2 SO 4 ) were used 25 to achieve the desired salinity. The salinity was checked using a calibrated refractometer. To 26 avoid precipitation in the final media only nitrogen (NaNO 3 ) was added before the water was 27 autoclaved (90 min., 121°C). Stock solutions of phosphorus (KH 2 PO 4 ), iron 28 Both strains of P. multiseries were grown in duplicate under N deficiency (36 µM 3 NaNO 3 , 36 µM K 2 HPO 4 ; atomic N:P 1:1) and P deficiency (580 µM NaNO 3 , 2 µM K 2 HPO 4 ; 4 atomic N:P 290:1). All eight 1-L chemostat culture flasks, containing 920 -996 mL culture, 5

were aerated with sterile air (0.2 µm pore size Sartorius Minisart® single use syringe filter 6 unit) to keep algal cells in suspension, and to ensure mixing as new medium was added drop-7

wise to the surface with peristaltic pumps (see below). Measurements of temperature, pH, 8 chlorophyll a (chl a), and microscopic observation of cells in all 8 cultures were made daily. 9

To avoid differences in media composition between culture flasks of one nutrient 10 treatment they were supplied with media from the same bottle. All cultures were kept under a 11 16:8 hours light: dark cycle at 200 µmol photons m -2 s -1 . The dilution rates (i.e. cell division 12 rate, related to the level of nutrient deficiency) of the cultures were checked daily by 13 measuring the out-flow of each culture. The first dilution rate was set to 0.2 d -1 (most 14 deficient). As soon as steady state conditions had been achieved in all cultures, a sample 15 ( 50% of volume) from each culture was taken (see below for details about sub-sampling) 16

and the dilution was thereafter adjusted to 0.3 d -1 . Sampling was repeated and the dilution rate 17 was finally adjusted to 0.4 d -1 (least deficient). The final sampling was conducted after a total 18 of 119 days in chemostat culture (Fig. 2). We considered the cultures to be at steady state 19 when the coefficient of variation of daily measurements of chl a of five consecutive days was 20 not greater than 8 % (usually < 5%). 21

The sample collected from each culture at steady state was sub-sampled in order to 22 determine abundance, volume, and morphology of the cells. Concentrations of chl a, 23

dissolved inorganic nutrients (N-NO 2 -+ N-NO 3 - , N-NH 4 + , P-PO 4 3-
, and reactive silicate) and 24 bacterial abundance in the cultures were also determined. Finally, total particulate carbon 25 (TPC), nitrogen (TPN), phosphorus (TPP), biogenic particulate silicon (BPSi) and domoic 26 acid (DA; particulate and dissolved) concentrations were analyzed (see below). 27 28

Analytical procedures and measurements 29 30

Samples for estimation of cell densities and determination of cell volume were 31 preserved with acid Lugol's solution. Cells were counted according to [START_REF] Utermöhl | Zur vervollkommnung der quantitativen phytoplankton methodik[END_REF], using 32 a Nikon Diaphot -300 inverted microscope. The samples were homogenized using a mini 33 shaker (MS2, IKA Tamro Med-Lab), which also disrupted any chains, before making dilution 34

V = (A × B × C × D) + [(1 -D) × A × (B + b)/2 × C]
Eq. 1 4

where V is the cell volume, A the length of the valve (apical axis), B is the valve width 5 (transapical axis in valve view), C is the valve depth (prevalvar axis in girdle view), D is the 6 ratio of the linear part of the cell length to the total cell length, and b is the cell breadth on its 7 extreme, corresponding to the minor base of the trapezoid. The BR strain had an average b = 8 0.5 µm, and b = 1.0 µm in the CA strain. Measurements were made using a light microscope 9

with phase contrast at 1000X magnification (ZEISS, model AXIOVERT). The microscope 10 was equipped with a digital camera and an image processing and measurement program 11 Bacterial abundance in the cultures was determined using a flow cytometer 23 (FACSCalibur, Becton Dickinson). Culture sub-samples of 1 mL (duplicates) were fixed with 24 pre-filtered (0.2 µm pore size filter) formaldehyde (2% final concentration) and stained for 30 25 min at room temperature in the dark with SYTO 13, a nucleic acid stain (Molecular Probes) at 26 a final concentration of 5 µM, following the method by [START_REF] Troussellier | Flow cytometric 23 discrimination of bacterial populations in seawater based on SYTO 13 staining of 24 nucleic acids[END_REF]. [START_REF] Valderrama | Methods of nutrient analysis[END_REF] 28

(

Nutrient analyses 29 30

Analyses of dissolved inorganic nutrients in the culture filtrates and inflow media 31 followed the methods modified by [START_REF] Valderrama | Methods of nutrient analysis[END_REF]. A Whatman GF/C glass fiber filter 32 was used to obtain filtrate for N and P analyses. An all-plastic filtration apparatus and 33 samples for analyses of the particulate fraction were filtered onto pre-combusted (2 hours, 1 450°C) Whatman GF/C glass fiber filters for total carbon, nitrogen and phosphorus (TPC, 2 TPN and TPP, respectively). The filters were then dried (24 hours, 60°C) and analyzed using 3 a CHN elemental auto analyzer (Fisons Instruments model NA 1500) for TPC and TPN. 4

Determination of TPP followed the method of [START_REF] Solórzano | Determination of total dissolved phosphorus and 8 particulate phosphorus in natural waters[END_REF]. Analysis of BPSi 5 followed the method originally described by [START_REF] Paasche | Silicon content of five marine plankton diatom species measured with a 20 rapid filter method[END_REF]. The filters, with algal cells, were 6 stored in plastic petri dishes at -20°C, and dried (2 hours, 80°C) between thawing and 7

analyses. 8 Domoic acid in the cells was determined by filtering duplicate sub-samples onto 9

Whatman GF/F glass fiber filters. The filters were freeze-dried to avoid loss of DA during the 10 transport to the Department of Food Chemistry, Friedrich-Schiller University, Jena, Germany 11 for analysis. The filtrates of the above samples were also collected for dissolved DA 12 measurements. Sodium azide (NaN 3 , final concentration 10 mg L -1 ), was used to preserve the 13 samples, which were stored at + 4°C prior to shipping. The DA analyses (particulate and 14 Cell volume, calculated as mean of all nutrient treatments and dilution rates for each 30 strain (n = 12), was higher (835.4 ± 108.3 µm 3 ) for the BR strain (apical axis around 100 µm) 31 than the CA strain (262.2 ± 33.6 µm 3 , apical axis around 32 µm). There was no trend in final 32 cell abundance between treatments (Fig. 3). The smaller Canadian cells were expected to give 33 a higher cell yield, and the lack of difference in cell abundance between strains is probably 34 due to the fact that the CA strain had been in culture much longer than the BR strain (see 1 below). We observed a general increase of chl a concentration as the cultures became less N-2 deficient, ranging from 24 -113 µg L -1 . Both extremes were observed in the BR strain, with 3 the lowest concentration during N-deficiency at 0.2 d -1 and highest during P-deficiency at 0.4 4 d -1 (Fig. 2). There was a significant increase in chl a C -1 (2-way ANOVA, P < 0.001), as the 5 cultures of both strains became less nutrient deficient (Fig. 3). There was on the other hand a 6 decreasing trend in bacterial abundance with reduced level of nutrient deficiency in all flasks. 7

The highest abundance, 12.2 × 10 6 bacteria mL -1 , was found in the N-deficient BR culture at 8 0.2 d -1 (i.e. high deficiency) whereas the P-deficient CA culture had the lowest recorded 9 bacterial abundance, 1.5 × 10 6 mL -1 , at the intermediate level of deficiency (Fig. 3). 10 11

Nutrients 12 13

The amount of N cell -1 in the BR strain was higher in the P-deficient than in the N-14 deficient treatment (Fig. 4). This was not observed in the CA strain. The CA strain had 15

highest N content at the intermediate level of deficiency in both nutrient deficient treatments. 16 No significant difference in the amount of P per cell between the strains within the same 17 nutrient treatment was found, but a significant difference was observed between treatments 18

(2-way ANOVA, P < 0.001). The N-deficient cultures of both strains had 2 -4 times more P 19 per cell than the corresponding P-deficient cultures (Fig. 4). Cellular content of C and Si was 20 also 2 -4 times higher in the BR than in the CA strain, in particular during P-deficiency, and 21 this difference was significant (2-way ANOVA, P < 0.05) (Fig. 4). Both C and Si per cell 22 decreased as level of deficiency was reduced in the BR strain. 23

The atomic C:N and C:P ratios were significantly higher in the BR strain, with the 24 exception of the C:P ratio under N-deficiency, and this did not change due to level of 25 deficiency (2-way ANOVA, P < 0.05) (Fig. 5). The C:N ratios were similar for the same 26 strain under both N-and P-deficiency. This shows that the cells contained low amounts of N 27 even in the P-deficient treatment when N was supplied in excess to their need. The C:P ratios 28 however, show a "luxury consumption" of P under N-deficiency as the C:P ratios are below 29 the Redfield atomic ratio of 106. The C:Si ratios decreased in the N-deficient BR culture as 30 the culture became less nutrient deficient (Fig. 5). There was a significant difference in N:P 31 ratios between N-deficient and P-deficient cultures (2-way ANOVA, P < 0.01). N-deficient 32 cultures did not show a change in atomic N:P ratios due to level of deficiency, while there 33 was an increase in the two P-deficient cultures with highest ratio in the central point of the 34 5). However, the N:P ratios of all cultures growing under N-1 deficiency were below, and all but one culture growing under P-deficiency had N:P ratio 2 above, the atomic Redfield ratio. 3 4 Particulate, dissolved and total domoic acid, and domoic acid production rate 5 6 Domoic acid was detected in both strains and nutrient treatments (Fig. 6). The 7 particulate (PDA) concentration per cell and the total DA (TDA = dissolved + particulate 8 fractions) production cell -1 day -1 was higher in P-deficient than in N-deficient cultures for 9 both strains. The range in PDA concentrations at cellular N:P ratios ≤7.8 were 0.08 -0.26 pg 10 cell -1 with values ranging between 0.93 and 3.45 pg cell -1 for cellular N:P ratios ≥14.7 (Fig. 11 6). PDA in the BR strain cells ranged between 0.11 and 3.45 pg cell -1 , and PDA 12 concentrations in the CA strain ranged between 0.08 and 2.99 pg cell -1 (Fig. 6). Total DA 13 production cell -1 day -1 followed the same trend as PDA cell -1 . That is, there was a clear 14 difference in TDA production depending on the N:P ratios (0.04 to 0.07 and 0.40 to 0.97 pg 15 TDA cell -1 day -1 for cellular N:P ratios ≤7.8 and ≥14.7, respectively). Bacterial abundance in 16 the culture flasks did not seem to have an important effect on neither PDA cell -1 nor TDA 17 production cell -1 day -1 in this experiment since no trends could be observed (Fig. 6). Level of 18 nutrient deficiency (d -1 ) only had an effect on PDA cell -1 and TDA production in the P-19 deficient BR strain with a decrease in toxin as dilution increased. No such trend could be 20 observed in the CA strain or the N-deficient BR strain (Fig. 6). 21

22

DISCUSSION 23 24

In this paper we present results from chemostat culture studies comparing the cell 25 physiology and domoic acid production for two strains of Pseudo-nitzschia multiseries. The 26 strains were isolated from marine waters at two opposite latitudes: Canada (CA) and Brazil 27 (BR). The strains were grown in parallel under N-or P-deficiency at three different levels of 28 nutrient deficiency (0.2, 0.3 and 0.4 d -1 ). The CA strain had been cultured in laboratory for ca. 29 10 years, whereas the BR strain was isolated 6 months prior to the experiment. Diatoms are 30 known to decrease in cell size while in culture [START_REF] Davidovich | Sexual reproduction in the pennate diatoms 11 Pseudonitzschia multiseries and P. pseudodelicatissima (Bacillariophyceae)[END_REF]. Reduced cell 31 sizes, cell volumes (about 31 % of the BR strain), and number of cells per chain in the CA 32 cultures were therefore expected. The CA strain also had more morphological deformations 33 such as "lobed" or "sickle-shaped" cells, and some with one or more "undulations", than the A deficiency in N or Si in phytoplankton leads to an increase in cellular carbon content, 3 resulting from carbon accumulation when cell division ceases (Flynn and Martin-Jezequel, 4 2000), as well as under P-deficiency (Granéli and Johansson, 2003; Granéli and Salomon, 5 2010; Pan et al., 1996a;[START_REF] Sakshaug | Chemical composition and 5 alkaline phosphatase activity of nutrient-saturated and P-deficient cells of four marine 6 dinoflagellates[END_REF]. [START_REF] Harrison | Marine diatoms grown 12 in chemostats under silicate or ammonium limitation. III. Cellular chemical composition 13 and morphology of Chaetoceros debilis, Skeletonema costatum, and Thalassiosira 14 gravida[END_REF] observed that chl a:C 6 ratios were lower under N-deficiency than for NP-sufficient diatom cells (Skeletonema 7 costatum, Chatoceros debilis, Thalassiosira gravida) when growing in chemostats. This was 8 due to a decrease in chl a cell -1 (i.e. not increase in C) the more N-limited the cells were. In 9 line with the findings by Pan et al. (1996a) we also observed an increase in cellular chl a:C in 10 both strains and nutrient treatments as the cells became less deficient. The amount of C cell -1 11 in this experiment, on the other hand, decreased as the cells were subject to reduced 12 deficiency in the BR strain, something not observed for CA cells. Cellular N, P, C and Si in 13 the experiment reported here were generally higher in the BR strain; and, especially for the 14 BR strain, decreased as the level of nutrient deficiency decreased. The higher N, P, Si, and C 15 cellular concentrations in the BR strain can be explained by the larger size of this strain in 16 comparison with the CA strain. Furthermore, the elemental ratios measured for both strains 17 are all within the range of what has been reported for P. multiseries by others (e.g. Pan et al., 18 1996a,b). production, even when the strains had been isolated from the same location. Our results show 1 that, not only when submitted to P-deficiency but also even under N-deficiency and when 2 there is a continuous cell division (supply of N even in small amounts), both BR and CA 3 strains produced DA. However, DA production was several-fold greater at cellular N:P ratios 4 above 14.7 (P-deficiency and N-sufficiency) than below 7.8 (N-deficiency and P-sufficiency). 5

Further, investigations of environmental conditions that may promote Pseudo-nitzschia spp. 6 blooms, affect the physiology, and/or trigger/enhance DA production are therefore necessary. depletion, as reported by Pan et al. (1996b). The results reported here resemble those reported 20 for two clones of Pseudo-nitzschia sp. cf. pseudodelicatissima isolated from the northern Gulf 21 of Mexico by [START_REF] Pan | Pseudo-nitzschia sp. cf. pseudodelicatissima -a confirmed producer of 25 domoic acid from the northern Gulf of Mexico[END_REF]. Their data suggest that a rapidly growing population, i.e. not a 22 stressed population due to nutrient deficiency, yields maximum net DA production rates and 23 DA cell quotas. Production of DA at steady state was also found to be positively correlated to 24 the N:P ratios, suggesting that P-deficiency promoted DA production (Pan et al., 1996a). We 25 found that N was indeed a requirement for high DA production for both the Canadian and 26

Brazilian strain during chemostat culturing since the amount of toxin was low at cellular N:P 27 ratios ≤7.8, i.e. N-deficiency. 28

A literature review reveals a great variation in P. multiseries cellular DA. Pan et al. 29 (1996a) report a range of 0.3 -4.5 pg DA cell -1 during P deficiency. Furthermore, they argue 30 that the optimum N:P ratio for P. multiseries is not the Redfield atomic ratio of 16:1 but that it 31 is closer to 10:1, and that DA production was only halted at N:P ratios less than about 7. This 32 is in line with the findings in our experiment reported here, i.e. 0.08 -3.45 pg DA cell -1 with 33 highest concentrations at high N:P ratios. There was no general trend in DA concentration as dilution rate changed, i.e. N:P ratio seem to have greater impact on DA production than 1 growth rate. [START_REF] Bates | Domoic acid toxicity of large new cells of 30 Pseudo-nitzschia multiseries resulting from sexual reproduction[END_REF], investigating DA toxicity of large new cells of P. multiseries 2 after induced sexual reproduction in cultures, observed that the parental clones had almost lost 3 their capacity to produce DA (< 0.1 pg DA cell -1 ) while their offspring cells (large new cells) 4

were able to produce up to 33 to 67 (5 months old clones) and 65 pg DA cell -1 (one 12 months 5 old clone). According to our results, culture age was not the most important factor since the 6 10-year older CA isolate produced similar levels of DA as the 6 months old isolated BR strain 7 (0.08 -2.99 and 0.11 -3.45 pg DA cell -1 , respectively). DA production in this experiment 8 (0.04 -0.97 pg DA cell -1 day -1 ) was also found to be similar to what has been reported 9 elsewhere. Lewis et al. (1993) found 0.01 pg DA cell -1 day -1 , and Subba Rao et al. (1990) 10 report a maximum of 2.0 pg DA cell -1 day -1 . This variation in DA cellular concentration is 11 within the range found for other Pseudo-nitzschia species [START_REF] Fehling | Growth and domoic acid 24 production by Pseudo-nitzschia seriata (Bacillariophyceae) under phosphate and 25 silicate limitation[END_REF]. 12 Other authors have found a positive correlation between associated bacterial abundance 13 and DA production. That is, non-axenic cultures had a higher DA production rate than axenic 14 cultures [START_REF] Douglas | Production of domoic acid, a neurotoxin amino acid, by 14 an axenic culture of the marine diatom Nitzschia pungens f. multiseries Hasle[END_REF][START_REF] Douglas | Domoic acid production by 17 axenic and non-axenic cultures of Nitzschia pungens f. multiseries[END_REF]. [START_REF] Stewart | Bacterial involvement in determining domoic acid levels in 10 Pseudonitzschia multiseries cultures[END_REF], studying bacterial 15 contribution to DA levels in culture, concluded that the total amount of DA at any one given 16 time is dependent on the Pseudo-nitzschia spp. production rate, nutrient status, bacterial DA 17 production stimulation, and bacterial degradation of DA. In our experiment, however, the 18 bacteria associated with P. multiseries did not seem to affect DA production, neither 19 positively nor negatively. 20

21

SUMMARY 22 23

Our results show that the level of deficiency did have a significant effect on atomic 24 N:C, P:C, Si:C and N:P ratios in the P. multiseries cells from Brazil or Canada. Furthermore, 25 DA concentrations per cell and production rates by the CA and BR strains are in the range of 26 those observed previously by other authors, across all three levels of deficiency and N:P ratios 27 used in our experiments. However, domoic acid concentration per cell, and DA production 28 per cell per day, were significantly different depending on level of nitrogen or phosphorus 29 deficiency, being higher under P-deficiency as expected due to the high content of N in the 30 toxin. However, DA was always produced in all treatments by P. multiseries in this chemostat 31 culture experiment. Hence, devastating DA contamination of, for example, shellfish can 32 occur at a wide range of growth rates and, in particular, nutrient ratios during and after a 33 prolonged bloom. The elemental value for diatoms (Brzezinski, 1985) of each nutrient ratio is indicated with a horizontal line. Other explanations as in Fig. 2.

-

  [START_REF] Thessen | Intra-and interspecies differences in 16 growth and toxicity of Pseudo-nitzschia while using different nitrogen sources[END_REF] studied the effect of NO 3 of three Pseudo-nitzschia species in batch culture.Kudela et al. (2004) 31 studied the effect of Si-limitation on different clones of P. australis and P. multiseries 32 growing in batch-, semicontinuous-, and continuous cultures. They also studied the effect of 33

( 32 M

 32 FeCl 3 6H 2 O/EDTA), and trace metals were autoclaved separately. Silicon (Na 2 SiO 3 5H 2 O, 29 always 140 µM final concentration to avoid limitation), vitamins, and the salt NaHCO 3 were 30 filtered through a 0.2 µm pore size non-pyrogenic and hydrophilic sterile single use syringe 31 filter unit (Sartorius Minisart®). Finally, media pH was adjusted to 8.1 ± 0.1 with sterile 1.2 HCl and/or 2 M NaOH, all done aseptically. 33

  and counting. Cell volume estimation of P. multiseries was based on frustule 1 measurements, assuming a rectangular box shape with double trapezoidal prisms in their 2 extremes[START_REF] Maldonado | The effect of Fe and Cu 14 on growth and domoic acid production by Pseudo-nitzschia multiseries and 15 Pseudonitzschia australis[END_REF]. The cell volume 3

  SPOT Insight, Diagnostic Instruments Inc., version 3.5 for Windows). Cell volume 12 estimation was conducted on the samples preserved in Lugol's solution at the Phytoplankton 13 Laboratory of the Fundação Universidade Federal de Rio Grande (FURG), Brazil. 14 Chlorophyll a was determined by filtering sub-samples through a Gelman Sciences type 15 A/E glass fibre filter. Extraction was done in 6 mL 95% ethanol according to Jespersen and 16 Christoffersen (1987) with the following method modifications. The filter in ethanol was 17 sonicated in an ice-cooled water bath in the dark for 15 minutes. Extraction was allowed to 18 complete under quiescent conditions for at least 2 hours in the dark. Prior to modification of 19 the original method, tests for comparison with 6 hour extraction in room temperature were 20 done and no significant difference was observed (data not shown). Chlorophyll a was 21 measured using a Turner Designs Fluorometer model 10 AU. 22

19 20 21 2 - 25 RESULTS 26 27

 192022526 dissolved) followed[START_REF] Furey | 29 Determination of domoic acid in shellfish by liquid chromatography with electrospray 30 ionisation and multiple tandem mass spectrometry[END_REF] using an HPLC coupled to a mass-spectrometer (API 15 165, Applied Biosystems/MDS Sciex). Total DA was calculated as the sum of particulate and 16 dissolved DA per cell. Domoic acid production rate (pg TDA cell -1 day -1 ) was calculated by 17 multiplying total DA concentration (pg TDA cell -1 ) with dilution rate per day. 18Statistics way ANOVA, using GraphPad Prism 5 for Mac OS X, was used to compare results 22 between level of deficiency and algal strain. Results were considered significantly different at 23 P < 0.05. 24 Algal cell volume, numbers, chlorophyll a concentrations and bacterial abundance28 29

  deficiency (Fig.

11 BR

 11 strain. Thus we considered differences in cell morphology observed between the two 1 strains a consequence of their culture age before the start of the experiment. 2

19 20Nutrient ratios and Domoic acid 21 22 A

 1922 number of studies have investigated the effect of various sources of N (e.g. urea, the physiology of Pseudo-nitzschia spp. and strains in batch cultures(Bates, 24 1998;[START_REF] Bates | Ecophysiology and metabolism of ASP toxin production[END_REF] Trainer, 2006, Howard et al., 2007;[START_REF] Thessen | Intra-and interspecies differences in 16 growth and toxicity of Pseudo-nitzschia while using different nitrogen sources[END_REF]. Laflamme (1993) 25 studied the effect of Si-and P-limitation on DA production of P. multiseries in chemostats. 26[START_REF] Kudela | Laboratory analyses of nutrient stress and 30 toxin accumulation in Pseudo-nitzschia species from Monterrey Bay, California[END_REF] used chemostat cultures to determine DA production of P. multiseries 27 strains, from one location, during Si-and N-deficiency. However, until now the effect of N-28 deficiency on general physiology and DA production in two P. multiseries strains isolated 29 from the NW and SW Atlantic Ocean, and growing in parallel chemostat cultures, has not yet 30 been investigated. Since DA is an amino acid, and N is an essential element required for its 31 synthesis, it has been argued that DA production would cease during N-deficiency(Bates et 32 al., 1991). However,[START_REF] Kudela | Laboratory analyses of nutrient stress and 30 toxin accumulation in Pseudo-nitzschia species from Monterrey Bay, California[END_REF] detected low levels of DA in N-limited chemostats. 33

7

 7 Among these are, light conditions[START_REF] Bates | The effects of iron on domoic 33 acid production by Pseudo-nitzschia multiseries[END_REF][START_REF] Fehling | Growth dynamics of non-toxic 21 Pseudonitzschia delicatissima and toxic P. seriata (Bacillariophyceae) under simulated 22 spring and summer photoperiods[END_REF], iron sufficiency or 8 deficiency[START_REF] Bates | The effects of iron on domoic 33 acid production by Pseudo-nitzschia multiseries[END_REF][START_REF] Rue | Domoic acid binds iron and copper: a possible role for the 3 toxin produced by the marine diatom Pseudo-nitzschia[END_REF][START_REF] Wells | Domoic 31 acid: The synergy of iron, copper, and the toxicity of diatoms[END_REF], several other trace 9 metals[START_REF] Wells | Domoic 31 acid: The synergy of iron, copper, and the toxicity of diatoms[END_REF], and the studies must be conducted under identical conditions on 10 several strains. Continuous cultures as an instrument to understand algal ecophysiology have 11 an advantage over batch cultures by allowing a dynamic equilibrium between nutrient input, 12 and cell death and growth; similar in principle to the conditions in natural waters. As in 13 nature, nutrient input rate during continuous culturing, rather then their concentrations at a 14 given time, determines growth limitation[START_REF] Rhee | Continuous culture in phytoplankton ecology[END_REF]. 15There are several examples from the literature where N-and/or P-deficiency cause an 16 increase in cellular toxin concentration in a number of algal species from different taxa, 17 including P. multiseries (Johansson 2003; Granéli and Flynn, 2006; Granéli and Salomon, 18 2010). It has been argued that DA production in Pseudo-nitzschia spp. is associated with Si 19

2 3 2 diatoméeFig. 1 .

 221 Fig. 1. Phylogenetic analysis of 18S rRNA gene sequences from Pseudo-nitzschia spp. and other diatoms. The Jukes and Cantor method was used to estimate the distances. The numbers at the branches indicate the percentage of 500 bootstraps; resamplings that support the branch values below 50% are not shown. The Pseudo-nitzschia sp. sequences are a sister group to Nitzschia apiculata. Note that the three Brazilian strains of Pseudo-nitzschia pungens (PNP1), Pseudo-nizschia pungens (PNP2) and Pseudo-nitzschia pungens (PNP3) sequences form one tight cluster and the Brazilian Pseudo-nitzschia multiseries (PNM1) and Canadian Pseudonitzschia multiseries (CCMP 1660) form a second tight cluster with zero branch lengths. The scale indicates the observed sequence divergence of 10%. Thus the divergence between Brazilian and Canadian P. multiseries is less than 1%.

Fig. 2 .

 2 Fig. 2. Chlorophyll a concentrations (µg L -1 ) of Pseudo-nitzschia multiseries (a Canadian and a Brazilian strain) growing under nitrogen (N) or phosphorus (P) deficiency in chemostat cultures during steady state (mean ± SD, n=2). Error bars are occasionally smaller than the size of the symbol, and symbols are superimposed. The dilutions used were 0.2, 0.3 and 0.4 per day in order to obtain high (0.2 d -1 ), intermediate (0.3 d -1 ), and low (0.4 d -1 ) N or P deficiency. Only the days when steady-state was obtained for each treatment is plotted.

Fig. 3 .

 3 Fig. 3. Observed Pseudo-nitzschia multiseries cell abundance (upper panel), chlorophyll a per unit cellular C (middle panel), and bacterial abundance (lower panel) in relation to level of N or P deficiency. Other explanations as in Fig. 2.

Fig. 4 .

 4 Fig. 4. Intracellular concentrations of N, P, C, and Si in relation to level of nutrient deficiency.Other explanations as in Fig.2.

Fig. 5 .

 5 Fig. 5. Elemental ratios of C:N, C:P, C:Si, and N:P in relation to level of nutrient deficiency.

Fig. 6 .

 6 Fig. 6. Concentration of Domoic acid (PDA) cell -1 , and total domoic acid cell -1 day -1 (TDA = particulate + dissolved DA in the culture bottles) in relation to cellular (atomic) N:P ratios (left hand panels), bacterial abundance (middle panels), and dilutions (right hand panels) (mean, n=2).

http://mc.manuscriptcentral.com/jplankt Journal of Plankton Research

Verlag, Berlin and Heidelberg, Germany, pp. 81-93.

Brzezinski, M.A.(1985) The Si:C:N ratio of marine diatoms: Interspecific variability and the

effect of some environmental variables. J. Phycol., 21, 347-357.

Springer-Verlag, Berlin and Heidelberg, Germany, pp. 229-241.

Prymnesium parvum cells grown under N-or P-deficient conditions. Harmful Algae, 2,

135-145.

Granéli, E. and Salomon, P., (2010) Factors influencing allelopathy and toxicity in

Oceanography, and Intergovernmental Oceanographic Commission of UNESCO, Paris,

Landsberg, J.H.(2002) The effects of harmful algal blooms on aquatic organisms. Rev. Fish.

Sci., 10, 113-390.