Di Ge

Lu Zhang

Christine Cavaro-Ménard

Patrick Le Callet

Numerical Stability issues on Channelized Hotelling Observer under different background assumptions

This paper addresses the numerical stability issue on the Channelized Hotelling Observer (CHO). The CHO is a well-known approach in the medical image quality assessment domain. Many researchers have found that the detection performance of the CHO does not increase with the number of channels, contrary to expectation. And to our knowledge, nobody in this domain has found the reason. We illustrated that this is due to the ill-posed problem of the scatter matrix and proposed a solution based on Tikhonov regularization. Although Tikhonov regularization has been used in many other domains, we show in this paper another important application of Tikhonov regularization. This is very important for researchers to continue the CHO (and other channelized model observer) investigation with a reliable detection performance calculation.

I. INTRODUCTION

The necessity and the importance of objective image quality measurement has long been acknowledged in the medical imaging community [START_REF] Barrett | Objective assessment of image quality: Effects of quantum noise and object variability[END_REF][START_REF] Barrett | Objective assessment of image quality. II. Fisher information, Fourier crosstalk, and figures of merit for task performance[END_REF][START_REF] Barrett | Objective assessment of image quality. III. ROC metrics, ideal observers, and likelihood-generating functions[END_REF][START_REF] Rolland | Effect of random background inhomogeneity on observer detection performance[END_REF]. To this end, numerical observers are used [START_REF] Narayan | Prediction of human observer performance by numerical observers: an experimental study[END_REF][START_REF] Abbey | Human-and modelobserver performance in ramp-spectrum noise: effects of regularization and object variability[END_REF][START_REF] Bruyant | Human and numerical observer studies of lesion detection in ga-67 images obtained with map-em reconstructions and anatomical priors[END_REF] to predict the outcome of a human observer performing certain diagnostic tasks (e.g., detection, localization). Their usefulness relies on the approximation to the performances achieved by human observers and in particular by clinical experts. Studies of numerical observers are largely motivated by their being less expensive and cumbersome than human observers.

The detection task is a binary classification with two underlying hypotheses (presence or absence of a signal in the image). The ideal observer (IO) was firstly proposed for the optimization of the binary classification problem, by choosing the hypothesis resulting in the greatest likelihood of the image data. The IO achieves the maximum AUC (area under the Receiver Operating Characteristic curve) among the observers, but the computation of likelihoods are tractable only for simplified cases. The Hotelling Observer (HO, aka Optimum Linear Discriminant) on the other hand maximizes the SNR among all linear observers while using only the first and secondorder statistics of the image data, i.e., the mean and the covariance matrix of the vectorized image data. The inversion of the (large) covariance matrix and its accurate estimate are the main difficulties in practical applications. To avoid the inversion problem of big matrix, data reduction method based on Channelized Hotelling Observers (CHO) was proposed to approach the SNR performance of the HO under sufficient statistics conditions [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF]. The channel design is also a mathematical model (based on filter banks [START_REF] Daly | The visible differences predictor: an algorithm for the assessment of image fidelity[END_REF][START_REF] Watson | Model of visual contrast gain control and pattern masking[END_REF] ) to approach the human perceptual and cognitive behaviors. Despite its simplicity in numerical implementation and flexibility in channel design to incorporate human per-ception models, the CHO yields counter intuitive results when it is implemented blindly. Increasing the number of channels could worsen the detection performance, making it difficult to choose the optimal channel number with respect to human experts' diagnoses. since the stable rela- tionship is a prerequisite in making such decisions. In the present study, we investigate the cause of the numerical instability and find that it is due to the direct inversion method applied to ill-conditioned scatter matrices. We propose a regularization method well established in the world of linear algebra for the CHO implementation and its performances are shown to be more stable while increasing channel numbers. To illustrate the numerical instability and the efficiency of our proposed regularization method, simulation tests are performed on signal known exactly (SKE) detection task in which all the signal attributes (e.g. intensity amplitude, size, shape, orientation and location) are known and do not vary throughout the experiment.

II. PROBLEM FORMULATION

We describe briefly here the CHO in the SKE detection task and the derivation of its scalar test statistic. Specific attentions are given to the problem of performance instability w.r.t. increasing channel numbers. The detection problem is considered as the validation of one of the two exclusive hypotheses : H 0 (signal absent) and H 1 (signal present). The observed image data g is given by:

H h : g = hx + b, h = 0, 1 (1)
where the known signal is denoted by x, the background by b, and the absence or the presence of the signal is controlled by the binary variable h. In the following, a 2D image (test image data, signal or background) is represented by a column vector using vertical concatenation. Let the observed image have √ M × √ M pixels, its vectorized version g is then a M × 1 vector.

A. Signal model and CHO

The current work considers the basic detection task of a rotationally symmetric nodule signal profile in image quality assessment and the image generation process follows [START_REF] Abbey | Human-and modelobserver performance in ramp-spectrum noise: effects of regularization and object variability[END_REF]. The background correlation structure and the imaging system are both rotationally symmetric.

We use the 2D Gaussian signal model proposed in [START_REF] Zhang | The effects of anatomical information and observer expertise on abnormality detection task[END_REF][START_REF] Goossens | The use of steerable channels for detecting asymmetrical signals with random orientations[END_REF] :

[x] p = a exp - 1 2 (p -q) t D -1 (p -q) (2)
where both p and q are 2-dimensional column vectors containing the sampling grid index.

[x] p denotes the signal intensity at p. The signal intensity attenuation is hereby modeled by a Gaussian function of peak amplitude a, centered at q. The diagonal matrix D such that diag(D) = [bσ 2 , σ 2 ] t specifies the ellipse's scale σ and shape b (b = 1 in the spheric case). Note that when all signal parameters are known exactly in advance by observers, the detection task is said to be a SKE task. Four artificial background models are mostly used in the literature: white Gaussian background, correlated Gaussian background, lumpy background [START_REF] Park | Performance of a channelized-ideal observer using laguerre-gauss channels for detecting a gaussian signal at a known location in different lumpy backgrounds[END_REF], and clustered lumpy background [START_REF] Castella | Mammographic texture synthesis: second-generationclustered lumpy backgrounds using agenetic algorithm[END_REF]. Note that all of these backgrounds are rotationally symmetric in the statistical sense. The superposition of the signal and the background yields simulated test images, presented in Fig. 1. The CHO pre-processes an image by channelization (dimension reduction by matrix multiplication): In practice, the SKE CHO can be implemented by a training stage and a test stage. In the training stage, two sets of image are considered : signal absent and signal present. The output template writes :

g ′ = U t g (
w CHO = S ′ -1
x ′ (4) using channelized signal and scatter matrix estimate ( x ′ and S ′ respectively) :

x ′ = ḡ1 ′ -ḡ0 ′ , (5)
S ′ = 1 2 (Var(g ′ 1 ) + Var(g ′ 0 )) (6)
The empirical means and covariances are noted as • and Var(•) respectively while the subscript • 0 and • 1 denote the signal class of the input signal. An inversion of the p × p matrix S ′ in Eq. ( 4) is needed.

In the test stage, the above-trained template w CHO is used to calculate a test statistic λ CHO linearly w.r.t. g ′ :

λ CHO = w t CHO g ′ (7)
A figure of merit can be calculated in order to characterize the CHO detection performance, such as :

SNR = ( λ0 -λ1 ) 2 (Var(λ 0 ) + Var(λ 1 ))/2 (8)
where test statistics are obtained for both classes (λ 0 and λ 1 respectively). Another widely used figure of merit for binary classifier (AUC) can be obtained by summing the area under the ROC (receiver operating characteristic) curve. In the following, we choose to show graphically the performance vs number of channels using the SNR measure since it is more capable of quantifying the separation of the two classes if the measure of the AUC is sufficiently close to unity, as argued in [START_REF] Park | Incorporating Human Contrast Sensitivity in Model Observers for Detection Tasks[END_REF]. Note the two measures are equivalent if the test statistics follow the independent and identically distributed (i.i.d.) Gaussian distributions [START_REF] Barrett | Objective assessment of image quality. III. ROC metrics, ideal observers, and likelihood-generating functions[END_REF]. Gallas et al. [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF] has shown that the CHO maximizes the SNR among all linear observers with reduced dimensions under the sufficient statistics conditions, i.e., channels should be chosen to effectively represent the HO template.

B. CHO simulation results

Gallas et al. [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF] has proposed the Laguerre-Gauss (LG) channels to adapt to the rotational symmetric properties of the signal and background model. The LG functions are the product of Laguerre polynomials and Gaussian (exponential) weightings. By combining the Gaussian decay, the LG functions in polar coordinates writes:

u n (r | a U ) = √ 2 a U exp - πr 2 a 2 U L n 2πr 2 a 2 U ( 9
)
where r ∈ Ê + is the radial distance variable and a u the spread parameter of the LG channel. These functions are then sampled and vectorized to form the column vectors U n of U in Eq. ( 3). It is argued in [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF] that an optimal choice of the parameter a U is related to the spread of the signal and the correlation length of the lumpy background. 2-D LG functions in Eq. ( 9) are calculated with different values of a U for a 128×128 image and are shown in Fig. 2. L-G channels having smaller a U tend to have more local effects. Performances of CHO detection using direct inverse of the scatter matrix (cf Eq. ( 4)) are shown in Fig. 3. L-G channel numbers vary from 2 to 40 while ranges of the spread parameter a u are adjusted to include the optimal value in each case. Two Gaussian signal profiles (with σ = 2 and σ = 5 respectively) are tested over the four backgrounds in Fig. 1. In all cases, 2000 training images and another 2000 testing images are generated in simulation. forming a sufficient database to accurately estimate the inverse covariance matrix considering the learning curve by Chan Based upon this simple simulation dataset, we show the difficulties in explaining test results of classical CHO implementations if the number of channels increases.

While the white Gaussian background allows relative performance stability for all L-G channels, the detection performances do not necessarily improve by increasing the channel numbers for the other background sets. On the contrary, important downgrades are frequently observed in Fig. 3 (2)-( 4) and this is a general result since no particular channel numbers do worse than others. However it is expected that by increasing the channel numbers 9)) and the y-axis ranges are adjusted for each case. The Channel numbers all range from 2 to 40.

an asymptote in the figure of merit should correspond to the performance of the ideal linear observer, i.e., without any dimension reduction [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF]. The problem is that the inverse of the covariance matrix S ′ -1 in Eq. ( 6) is not accurately estimated since more channels require more cases, as was demonstrated by Chan et al. [START_REF] Chan | Classifier design for computer-aided diagnosis: Effects of finite sample size on the mean performance of classical and neural network classifiers[END_REF] and in some cases it may not even be invertible (see discussion on a u in the next section). While L-G channels with bigger a u tend to yield more stable performances in general, in some cases such as Fig. 3 (2) smaller a u could have interesting potentials if it were not for the stability issue. Note the zoom-in effect of the y-axis for each scenario in Fig. 3.

This result is particularly frustrating since the L-G channels are chosen because continuous L-G functions form a basis on the space of rotationally symmetric square-integrable functions in Ê 2 : any rotational sym- metric function y(r) can be exactly expanded by :

y(r) = ∞ n=0 α n u n (r | a U ),
with α n the inner product defined by :

α n = a U √ 2 y(r)u n (r | a U )dr.
It is expected that the more channels are applied, the better detection performance should be achieved. We investigate in the following the numerical cause and propose a method accordingly. Recall that the CHO as well as any model observers are developed to be comparable to human observers and not to surpass them in performances. It is worth noting that the motivation of this study is the maximization of the detection performance of the L-G CHO, and by doing so achieve a stabilized relationship w.r.t. the channel numbers so that the choice of parameters like {p, a U } based on the comparison with human observers' performances is unbiased.

III. MATRIX CONDITION AND INVERSION

To our knowledge, no previous study has been made to achieve reasonable detection performances w.r.t. increasing channel numbers even though it is necessary priorly to the determination of optimal channel parameter. Indeed, the inversion of the covariance matrix ( S ′ ) becomes more sensitive to the spread parameter and finite sample size if more channels are added (curse of dimensionality). Chan et al. [START_REF] Chan | Classifier design for computer-aided diagnosis: Effects of finite sample size on the mean performance of classical and neural network classifiers[END_REF] studied the effect of finite sample size on classifier designs for computer-aided diagnosis (CAD) algorithms and showed empirically the linear dependences between the AUC and the reciprocal of the number of training samples. Generally speaking, the sample size available for classifier design is limited and thus introduces variance and bias into the detection performance. In the current study, we show that the "dimension curse" is arriving more rapidly than expected with small spread parameters (cf Fig. 4), such that increasing the number of training images to achieve better estimate of the inverse covariance matrix is no longer an option : either unrealistic in the case of medical imaging or too much time consuming in simulation tests. Apart from the sample number issue, the real covariance matrix itself could not be invertible (ill-posed). The L-G channels are not a basis on a sampled grid. What is worse, these channels become almost linearly dependent when the a u parameter is pushed to smaller values due to the sampling grid effect (see the aliasing for the case a u = 5, p = 17 in Fig. 2). This is typically what's happening in Fig. 4 with a u ≤ 1 : the condition numbers of the covariance matrix (the ratio between the biggest eigenvalue in absolute value and that of the smallest) start with high levels and are extremely fast growing afterwards. We also observe that :

1. the condition numbers grow exponentially by increasing p of L-G channels for all spread parameters;

2. they run quickly out of the dynamic range (10 20 ) given the float point precision ǫ = 10 -16 of the simulation settings, making them extremely illconditioned matrices (cond(•) ≫ ǫ -1 );

3. larger a U yields slower increase but remains exponential. However, smaller a U might have better performances in some cases (cf Fig. 3) if stabilized, which makes the regularized solution inevitable if we aim to optimize the CHO performances;

4. the other backgrounds yield similar results as that of the lumpy background case in Fig. 4. Higher order L-G channels are thus penalized by the fast growing condition number of the covariance matrix estimate in CHO applications. It is also interesting to note that without further adaptations, it is reasonable to privilege L-G channels with higher values for the spread parameter a U because of their numerical stability of the inverse covariance matrix S ′ . However, as shown in Fig. 3 and argued in [START_REF] Gallas | Validating the use of channels to estimate the ideal linear observer[END_REF] its optimal value is dependent on the spread of the known signal s and background correlation properties. This is also the reason why this numerical instability issue should be addressed so that the choice of the optimal a U and the performances vs channel number curves are unbiased. In the following, we propose a numerical fix to achieve a more meaningful curve for the figure of merit, the latter being the prerequisite to finding the optimal CHO model w.r.t. the human perceptual and cognitive tasks [START_REF] Daly | The visible differences predictor: an algorithm for the assessment of image fidelity[END_REF].

IV. REGULARIZED INVERSION

The empirical estimate matrix S ′ is by construction (cf Eq. ( 6)) non-negative definite and allows the following SVD decomposition :

S ′ = VDV t , (10)
where V is an orthonormal matrix containing normalized eigenvectors (VV t = I), and D = diag(d 1 , . . . , d p ) is a diagonal matrix whose main diagonal values are eigenvalues in the descending order d 1 ≥ . . . ≥ d p > 0. We recall that the condition number is the ratio

d 1 /d p . The direct inverse writes S ′ -1 = VD -1 V t .
The Tikhonov regularization is particularly adapted to the evaluation of moderate dimension and ill-conditioned problems such as Eq. ( 4) and provides a direct closedform solution. Instead of calculating Eq. ( 4) by direct inverse as the solution of a linear problem :

S ′ • w = x ′ , (11)
a bias is introduced deliberately by considering the following problem :

min 1 2 S ′ • w -x ′ 2 + 1 2 η w 2 , η > 0 (12)
for which the positive coefficient η controls the regularization level of the solution. From the SVD point of view, the regularization method can be applied to replace S ′ -1 in Eq. ( 4) with w † CHO = S † x ′ where

S † = S ′ t S ′ + ηI -1 S ′ t = V     d1 η+d 2 1 . . . dp η+d 2 p     V t and η is chosen such that d 2 1 ≫ η ≫ d 2 p
for the extremely ill-conditioned cases. We note that for S ′ with relatively small condition numbers, the regularization has negligible effects by imposing d 2 1 ≫ η : S † is close to the direct inverse of S ′ and thus the detection performances do not change significantly. On the contrary, for extremely illconditioned cases the presence of η avoids the division by zero. Indeed, the regularization technique attempts to achieve numerical stability by penalizing the template norm w CHO 2 so that the latter would not grow out of control as in the case of the direct inverse by Eq. ( 4). Evidently, bigger η values could further regularize the template solution norm, at the cost of a more important bias.

We present in Fig. 5 the figure of merit (SNR) for the detection performances on the simulated data used in Fig. 3. Detection performances are more smoothly increasing if not strictly monotonic. Most importantly, the regularized inversion method stabilizes L-G channels with relatively small spread parameters that in some cases are the optimal L-G CHO (cf the lower panel of Fig. 5(2) and the upper panel of Fig. 5(3)).

V. CONCLUSION

We presented in the current study the idea of using Tikhonov regularization method to address the extremely ill-conditioned scatter matrix inversion problem, the main cause of the numerical instability observed in the SKE detection tasks. Rotationally symmetric backgrounds and signals are chosen together with adapted L-G channels to illustrate the principle. It is indeed a generic problem : other channels such as the DoG channels and the Gabor channels that have shown capable of closely tracking human performances are equally subject to the phenomenon. Our goal is to achieve stable performance results with respect to the channel numbers for a given dataset such that the choice of optimal parameters should not be influenced. In that sense, the Tikhonov regularization method achieves that goal by obtaining a stable performance asymptote that characterizes the upper limit of linear observers.

FIG. 1 .

 1 FIG. 1. Simulated test images (128 × 128) with gaussian lesion profile in the center of different backgrounds: (a) White Gaussian, (b) Correlated Gaussian, (c) Lumpy, (d) Clustered lumpy.

  3) for which U M×p contains p channel (p ≪ M ) responses in each column. All calculations afterwards are per-formed on the channelized data, denoted by (•) ′ in the following.

FIG. 2 . 2 a

 22 FIG. 2.LG functions for aU = 5, 10, 20, 40 (column) and the order of channel = 0, 3, 9, 17 (line). The grayscale is normalized by the factor of √ 2a U for each column.

( 1 )( 3 )( 4 )

 134 FIG.3. L-G CHO detection performances (SNR) using different L-G channel numbers with direct inverse (by Gaussian elimination) of the scatter covariance matrix S ′ . Four backgrounds are tested, each with 2000 training image inputs and 2000 testing image. L-G spread parameters au (c.f. Eq. (9)) and the y-axis ranges are adjusted for each case. The Channel numbers all range from 2 to 40.

FIG. 4 .

 4 FIG. 4. Example of the condition number of the matrix S ′ with increasing dimension p for the five classes of L-G channels (au = [1 5 10 20 40]). The test is realized on a lumpy background identical to the lower panel (signal profile σ = 5) of Fig. 3(3)

( 1 )( 3 )( 4 )

 134 FIG.[START_REF] Narayan | Prediction of human observer performance by numerical observers: an experimental study[END_REF]. SNR with proposed Tikhonov regularization for all backgrounds with simulation data identical to Fig.3. The y-axis ranges are adjusted for each case given specific backgrounds and signal sizes.


Hafiz Ahmed

Rosane Ushirobira

Denis Efimov

Damien Tran

Jean-Charles Massabuau

Dynamical model identification of population of oysters for water quality monitoring

The measurements of valve activity in a population of bivalves under natural environmental conditions (16 oysters in the Bay of Arcachon, France) are used for a physiological model identification. A nonlinear auto-regressive exogenous (NARX) model is designed and tested. The model takes into account the influence of environmental conditions using measurements of the sunlight intensity, the moonlight and tide levels. A possible influence of the internal circadian/circatidal clocks is also analyzed. Through this application, it is demonstrated that the developed dynamical model can be used for estimation of the normal physiological rhythms of permanently immersed oysters and considered for detection of perturbations of these rhythms due to changes in the water quality, i.e. for ecological monitoring.

I. INTRODUCTION

Our world's environmental quality has been changing rather quickly and water quality everywhere is changing as well. Nowadays, national and international legislation has strict recommendations on protection of aquatic environment against the discharge of dangerous substances 1 . The realization of a protection program for aquatic systems needs a largely distributed network of ecological monitoring for water quality. Bio-indicators are increasingly used and show high efficiency, for instance in their ability to reveal the presence of low concentrations of contaminants through bioaccumulation in tissues 2 . However, if the goal is to sample at a high frequency, 24/7, with avoidance of missing short term events, most existing systems for ecological monitoring of aquatic environment are limited because they become rapidly too expensive and unrealistic. Indeed, frequently an intensive exploitation of human resources for collection of samples, chemical analysis and measurement of toxicity is involved [START_REF] Telfer | Review of environmental impact assessment and monitoring in aquaculture in europe and north america[END_REF]. A solution is to develop unmanned systems, able to work at high frequency by remote control, 24/7, either as bio-sensors by themselves (integrating information) or early warning detectors to trigger a sampling campaign. Today, networks of such online sensors, operating at large scale do not exist and are still a matter of research.

To fulfill that objective, an installation of numerous on-line remote sensors is required, working at high frequency for instant collection of information on a daily basis in marine environment [START_REF] Kröger | Sensing the sea[END_REF]. Behavioral and physiological responses of wildlife to pollution are very sensitive and can be estimated for an indirect ecological monitoring. However, a limiting 1 http://europa.eu/legislation_summaries/environment/water_protection_ management/l28002b_en.htm; http://water.epa.gov/scitech/swguidance/standards/criteria/aqlife/index.cfm; 2 For example, see the US Mussel Watch, http://ccma.nos.noaa.gov/about/coast/nsandt/welcome.html.

factor today is that an on-line and continual physiological monitoring of wildlife requires transmission of high volumes of data, whose analysis needs sufficiently accurate models of animal behavior in natural conditions [START_REF] Galang | Real-time detection of water pollution using biosensors and live animal behaviour models[END_REF]. Other difficulties lie in the fact that animals may be heavily influenced by environment, group interactions and internal rhythms (e.g. feeding, breathing, spawning).

Observation of the opening and closing activities of bivalves is a possible way to evaluate their physiological behavior in reaction to environment. The deviations from a normal state can be used for detection of a contaminant in surrounding water. The pioneer work that analyzes bivalve's activities through the record of their valve movements (e.g., valvometry) was realized by Marceau [START_REF] Marceau | Contraction of molluscan muscle[END_REF] with smoked glazed paper. Today, valvometers are commercially available and are mainly based on the principle of electromagnetic induction, like the Mossel Monitor [START_REF] Kramer | The valve movement response of mussels: A tool in biological monitoring[END_REF] or the Dreissena Monitor [START_REF] Borcherding | The Zebra Mussel Dreissena Polymorpha: Ecology, Biological Monitoring and First Applications in the Water Quality Management[END_REF]. In recent years, the interest for modeling and estimation of behavior of marine animals directly in real marine conditions has intensively increased [START_REF] Riisgard | Valvegape response times in mussels (mytilus edulis): Effects of laboratory preceding feeding conditions and in situ tidally induced variation in phytoplankton biomass[END_REF], [START_REF] Robson | Mussels flexing their muscles: A new method for quantifying bivalve behavior[END_REF], [START_REF] Garcia-March | Shell gaping behavior of pinna nobilis l., 1758: Circadian and circalunar rhythms revealed by in situ monitoring[END_REF], [START_REF] Galang | Real-time detection of water pollution using biosensors and live animal behaviour models[END_REF].

A remarkable monitoring solution has been realized in the EPOC CNRS UMR laboratory in Arcachon, France [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: Application to cadmium[END_REF], [START_REF] Sow | Water quality assessment by means of hfni valvometry and high-frequency data modeling[END_REF], [START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF], where a new framework for noninvasive valvometry has been developed and implemented since 2006. The designed method is strongly based on bivalve's respiratory physiology and ethology. The obtained measurements have been deeply analyzed using statistical methods handling very large amount of data, which is the key to analyze the bivalve's behavior at a large scale. The obtained data is available on-line on the site of the project3 . This system allows the bivalves to be studied in their natural environment with minimal experimental constraints. The obtained arrays of data of opening and closing activities of bivalves were used for analysis from different points of view using statistical approaches (estimation of probability density functions for time and amplitude of opening) [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: Application to cadmium[END_REF], [START_REF] Sow | Water quality assessment by means of hfni valvometry and high-frequency data modeling[END_REF], [START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF].

The goal of the present study is the identification of a physiological dynamical model for bivalves using EPOC's data. A population of 16 oysters living in the Bay of Arcachon (France) is studied, the influence of external forces (like the sunlight and the moonlight, the tide level) are taken into account. Since the behavior of oysters is also guided by internal circadian/circatidal rhythms, their influence is also incorporated in the developed NARX model. The type and the origin of internal circadian/circatidal rhythms of bivalves is rather uncertain, that is why a side result of the present study consists in verifying different hypotheses for modeling of circadian/circatidal rhythms for oysters. Van der Pol oscillator model has been used in this regard.

The outline of the paper is as follows. A brief description of the measurement scheme and experiments is given in Section 2. The model identification procedure is presented in Section 3. The verification of the model and its utility for an ecological monitoring of water quality are discussed in Section 4.

II. MEASUREMENT SYSTEM DESCRIPTION

The monitoring site is located in the Bay of Arcachon, France, at the Eyrac pier (Latitude: 44 • 40 N, Longitude: 1 • 10 W). Sixteen Pacific oysters, Crassostrea gigas, measuring from 8 cm to 10 cm in length were permanently installed on this site. These oysters were all from the same age group (1.5 years old) and came from the same local supplier. They also all grew in the Bay of Arcachon. They were immersed on the sea bottom (at 3 m to 7 m deep in the water, depending on the tide activity).

The electronic equipment has been first described in [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: Application to cadmium[END_REF], [START_REF] Chambon | Influence of the parasite worm polydora sp. on the behaviour of the oyster crassostrea gigas: A study of the respiratory impact and associated oxidative stress[END_REF]. A considerable advantage of the developed monitoring system is that it is completely autonomous without in situ human intervention for one full year. Each animal is equipped with two light coils (sensors), ≈ 53 mg each, fixed on the edge of each valve. These coils measure 2.5 × 2.5 × 2 mm and were coated with a resin sealing before fixation on the valves. One of the coils emits a high-frequency sinusoidal signal, that is received by another coil. Measurements are performed every 0.1 sec successively (with the frequency 10 Hz) for each of the sixteen animals. This means that the behavior of a particular oyster is measured every 1.6 sec. Every day, 54000 triplets (1 distance, 1 stamped time value, 1 animal number) are collected for each oyster. The strength of the electric field produced between the two coils is proportional to the inverse of distance between the point of measurement and the center of the transmitting coil, which leads to estimation of the distance between coils. After some scaling, the obtained distance measurements characterize a relative opening of valves for an animal, i.e. a value between 0 and 100 is received, where 0 and 100 correspond to complete closing and opening in percentage respectively.

On the shore, a second electronic unit manages the data acquisition and transmission. This unit is equipped with a GSM/GPRS modem and uses Linux operating system for driving the first control unit immersed in the water, managing the data storage with a time stamp, accessing the Internet, and transferring the data. At the end of a working day (or any programmed experimental duration), the collected data is transmitted to a central workstation server located in the Marine Station at Arcachon, France. The valve activity data is stored in a central database. A schematic description of the monitoring system is presented in Figure 1.

The collected tidal data includes a measurement of the height of the water column every hour, while the times and levels of low and high tides are provided by the hydrographic 4 . To take into account the external diurnal/nocturnal rhythms of the bivalve, the sun and the moon positions at the sampling site have been used. It is worth to note that data losses in transmission were sometimes observed due to impairment of data transfer.

III. MODEL IDENTIFICATION

Summarizing the description given in the previous section, after a specified duration of the experiment, the electronic system provides a matrix of values of relative opening of valves for each animal P i,j ∈ [0, 100] for 1 ≤ i ≤ N , 1 ≤ j ≤ n, where n = 16 is the number of oysters in the experiments and N is the number of measured points, N = 54000 × N days with N days > 0 the number of days used for identification of the model. In our case, the data of 2007 collected at the Eyrac pier in Arcachon has been used, so N days = 365. In addition, the values S i , M i , W i are provided for 1 ≤ i ≤ N and characterize the sun and the moon positions with respect to the horizon in degrees and tide levels in meters respectively. The signals S i and M i take negative values for the corresponding positions below the horizon line.

Since N ≈ 2 × 10 7 , to simplify computations and the presentation in this study, the measurements have been averaged on hourly basis, i.e. define δN = 3600 P q,j ,

s k = δN -1 kδN q=(k-1)δN +1 S q , m k = δN -1 kδN q=(k-1)δN +1 M q , w k = δN -1 kδN q=(k-1)δN +1 W q ,
correspond to hourly average opening of valves, the sun and the moon positions, the water level respectively. Examples of obtained signals p k,j for three animals (j = 1, 2, 3), s k , m k and w k are shown in figures 2 and 3 with zoomed versions.

As we can conclude from these plots, different oysters have largely different behavioral activity and the correlation between them is rather weak. Define the averaged opening

h k = n -1 n j=1 p k,j ,
the derived signal is presented in Figure 4. From this plot, we can deduce that the averaged on hourly basis openings of valves for population h k have a more regular and smoother behavior. This is a reason to consider the signal h k more suitable for the forthcoming analysis and model identification. Thus an averaged physiological population dynamics has to be modelled for oysters. A motivation of this choice is standard in modeling biological systems: each animal is largely perturbed by (stochastic) exogenous and endogenous forces, which are hard to predict or measure, while averaged signal h k is less sensitive to these perturbations and represents mainly the regular physiological behavior of the animals. Naturally, this is especially true if the population size n is big, but this is not the case in the considered application due to various experimental complexities forcing to limit the population size to 16. It is well known [START_REF] Tran | Field chronobiology in molluscan bivalves: How moon and sun cycles interactions deeply drive oyster activity rhythms[END_REF] that the opening activity of oysters h k is externally governed by the sunlight, moon oscillations and the tide activity, denoted here by s k , m k and w k respectively. Other sources of rhythmicity are internal (feeding, breathing etc.) and supervised by internal clocks. The reference signal generated by circadian rhythm genetics, which is directly influencing h k , will be denoted by c k . The generic structure of a dynamical physiological model for population of oysters can be presented as follows:

h k = F (h k-1 , . . . , h k-ν ; s k , . . . s k-µ ; (1) m k , . . . m k-µ ; w k , . . . w k-µ ; c k , . . . c k-µ ; θ) + ǫ k , c k = D(ξ k ), (2)
ξ k = G(ξ k-1 ; s k , . . . s k-µ ; m k , . . . m k-µ ; w k , . . . w k-µ ),
where ν ≥ 0 represents the number of past events taken into account by the animals in order to determine the opening distance at the next time instant k (i.e. in this model k -1 is the current time instant and h k-1 is the current valve state); µ ≥ 0 is the number of previous values for the positions of sun/moon and tide levels that are used in the decision on the value of h k (the model [START_REF] Telfer | Review of environmental impact assessment and monitoring in aquaculture in europe and north america[END_REF] assumes that the instant values of s k , m k and w k can be used for this decision); ǫ k ∈ R is the disturbance representing additional uncertain influencing forces and approximation errors (it is assumed that this term is sufficiently small if the model has been well identified); θ ∈ R q , q > 0 is the vector of constant parameters of the model ( 1), and F : R ν+4µ+4+q + → R + is a function defining the physiological model structure; ξ k ∈ R p is the state of the circadian oscillator (2); the functions G : R p+3µ+3 → R p and D : R p → R define the structure of the circadian rhythm model [START_REF] Kröger | Sensing the sea[END_REF]. The model ( 2) assumes that circadian oscillations are entrained by the external cues s k , m k and w k (or some of them depending on the type of the rhythm under consideration, see below). The structure scheme of the model is given in Figure 5.

The identification procedure consists in selecting an appropriate structure for the physiological model F (•) and the circadian model G(•), D(•), with posterior calculation of the corresponding vector of parameters θ.

Further in this section two problems are considered. First, the issue of circadian rhythms modeling is discussed and a solution is proposed. Second, an ARMAX structure is selected for F (•) and the corresponding model is designed.

A. Models of circadian clocks

A circadian rhythm is a biological process that displays an endogenous and self-sustained oscillation with a period of about 24 hours. These rhythms are driven by a circadian clock.Circadian rhythms have been observed in almost all forms of living organisms. There exist many mathematical models of circadian clocks [START_REF] Klerman | On mathematical modeling of circadian rhythms, performance, and alertness[END_REF], [START_REF] Smolen | A reduced model clarifies the role of feedback loops and time delays in the drosophila circadian oscillator[END_REF], [START_REF] Goldbeter | Computational approaches to cellular rhythms[END_REF]. Different hypothesis on the nature of clocks presented in oysters [START_REF] Tran | Field chronobiology in molluscan bivalves: How moon and sun cycles interactions deeply drive oyster activity rhythms[END_REF] can also be issued, but not yet supported by a mathematical model. First, like many other animals, oysters may have oscillations with the period of 24 hours externally governed by sunlight (i.e. circadian clock). Second, they may also have a second exogenous rhythm with the period 12.4 hours and entrained by tides (circatidal clock). Third, the circadian clocks may be regulated by the moon with the period 24.8 hours (circalunidian clock). The diversity of hypothesis follows the variety of habitation areas of bivalves, e.g. arctic zones, ocean conditions with strong tides or oceans where tides are almost negligible. Finally, all these factors may entrain the animals simultaneously, and the circadian clock may have a varying period from 12.4 to 24.8 hours. In this work, for brevity, only the first hypothesis is considered, and the corresponding mathematical model is developed below.

Conventional models of circadian oscillators [START_REF] Klerman | On mathematical modeling of circadian rhythms, performance, and alertness[END_REF], [START_REF] Smolen | A reduced model clarifies the role of feedback loops and time delays in the drosophila circadian oscillator[END_REF], [START_REF] Goldbeter | Computational approaches to cellular rhythms[END_REF] have different complexity. There exist also methods to design generic oscillating systems [START_REF] Efimov | Yakubovich's oscillatority of circadian oscillations models[END_REF] of any complexity. However, as it was observed in [START_REF] Forger | Reconciling mathematical models of biological clocks by averaging on approximate manifolds[END_REF] and [START_REF] Efimov | Controlling the phase of an oscillator: a phase response curve approach[END_REF], any oscillating system on the limit cycle can be fully characterized by its current phase and by the form of the limit cycle (i.e. amplitude of oscillations). If the exact genetic or biochemical nature of circadian oscillations is not clear, then to preserve the model simplicity, the circadian rhythm can be represented by any oscillator, and a particular attention has to be paid to a relation of the phase of this oscillator and environmental rhythms. In [START_REF] Antle | Gates and oscillators: a network model of the brain clock[END_REF] and [START_REF] Forger | Reconciling mathematical models of biological clocks by averaging on approximate manifolds[END_REF], variants of Van der Pol oscillator [START_REF] Kanamaru | Van der pol oscillator[END_REF] have been used for modelization of circadian rhythms for different animals. Van der Pol oscillator has a simple planar mathematical model with only two parameters that allows the form of the limit cycle and the period of oscillations to be tuned:

ξ1 = α[ξ 2 + µ(ξ 1 - 4 3 ξ 3 1 ) + u], ξ2 = -αξ 1 ,
where ξ 1 ∈ R and ξ 2 ∈ R are the model states, α > 0 determines the model oscillation velocity (adjustment of this parameter regulates the period of oscillation for the model); µ > 0 is a parameter to determine the shape of oscillations (µ ≈ 0 for almost harmonic oscillations, µ ≥ 2 leads to a nonlinear profile); u ∈ R is the model input that can be used for the entrainment. In our case the parameter α has to be selected in order to ensure the period of oscillations of 24 hours for the first hypothesis, for instance. The form of oscillations on the limit cycle, i.e. the value of parameter µ, is less important since in the model ( 1), (2) the shape of c k can be taken into account later in the function F (•) and by tuning the parameters θ (the value µ = 2 is used in this work for all computations). Therefore, the functions G(•), D(•) in (2) can be defined based on Van der Pol equations discretized using the Euler method:

ξ 1,k = ξ 1,k-1 + T α[ξ 2,k-1 + µ(ξ 1,k-1 - 4 3 ξ 3 1,k-1 ) + u k ], ξ 2,k = ξ 2,k-1 -T αξ 1,k-1 , c k = max{0, ξ 1,k },
where T = 1 hour is the sampling time, c k is the positive part of ξ 1 and u k is the resetting input, which has to be specified.

The entrainment regulates the phase of circadian oscillations (synchronizes the exogenous cues and endogenous rhythms), the mechanism of entrainment and phase resetting can be well analyzed using the Phase Response Curve approach [START_REF] Efimov | Controlling the phase of an oscillator: a phase response curve approach[END_REF], [START_REF] Efimov | Phase resetting control based on direct phase response curve[END_REF]. For example, for the first hypothesis we select u k = ρ max{0, s k }/ max 1≤k≤L |s k |, where ρ > 0 is the scaling parameters that is selected in order to harmonize the amplitude of the input u k and the dimension of the limit cycle. The selection max{0, s k } is applied since for the circadian rhythm the daylight is the principal zeitgeber, and the light is emitted when the sun is above the horizon only.

To generate the signal c k using the proposed model it is necessary to properly assign the initial phase for the circadian clock model (initial position on the limit cycle), which has to be coordinated with the current sun activity. A possible solution is to apply to the model the input s s pattern extracted from the first month of the experiment (or the last one in our case, since the experiment duration is 1 year and the sun position has annual periodicity) several times in order to entrain the oscillator.

B. ARMAX model

The ARMAX (auto-regressive-moving-average with exogenous inputs) model is one of the most popular structures used for identification in various fields of science [START_REF] Nelles | Nonlinear System Identification[END_REF]. Its advantages include linearity with respect to all signals and parameters, different effective methods for calculation of the parameters, simplicity of stability analysis of the obtained model and robustness with respect to perturbations. Therefore, the function F (•) in ( 1) is selected in the following form:

To evaluate the obtained accuracy of the designed model two performance costs J ǫ and J p have been calculated. The cost J ǫ estimates the average amplitude of ǫ k (the approximation errors in the model): The cost J p evaluates the prediction quality of the model on an infinite time interval:

ǫ j = h j+ν -ĥj+ν , ĥj+ν = Φ T j θ, 1 ≤ j ≤ L -ν, J ǫ = L-ν j=1 ǫ 2 j L -ν .
e j = h j+ν -hj+ν , 1 ≤ j ≤ L -ν, J p = L-ν j=1 e 2 j L -ν , where hk = h k , 1 ≤ k ≤ ν, hν+j = ν i=1 θ h i hν+j-i + µ z=0 θ s i sν+j-z + µ z=0 θ m i mν+j-z + µ z=0 θ w i w ν+j-z + µ z=0 θ c i c ν+j-z , 1 ≤ j ≤ L -ν
are the estimates of the valve positions h k generated independently by the designed model in the presence of the same inputs. Examples of the obtained estimates ĥk and hk (for ν = µ = 72, i.e. the oysters have 3 days of memory) are shown in Figure 6, J ǫ = 12.06 and J p = 23.18. As we can conclude from this example, the model demonstrates a sufficiently good accuracy of representation of the physiological behavior of oysters.

B. Application to ecological monitoring

Various problems related to an ocean ecological monitoring system are already mentioned in the introduction. Valve activity model of the oyster population can be used as a tool to overcome those problems. The area of population of bivalves is rather large, and the developed physiological model can serve for water quality monitoring, if we would compare the real measurements h k obtained by an embedded electronic unit on a bay and the estimates provided by the model hk . For example, the following daily monitoring residual can be computed:

r k = 1 24 23
z=0
|h k-z -hk-z | for all 24 ≤ k ≤ L, where an averaging on 24 hours is used in order to decrease sensitivity of the residual with respect to measurement noises. For the dataset collected at Since the output of a storm sewer is located at about 10-15 meters from the oysters, we suggest that some unknown contaminant, washed by the rain, could have reach them. This simple case study justifies the utility of the developed approach for ecological monitoring.

V. CONCLUSION

The paper presents the first development of a dynamical physiological model for oysters, which takes into account an influence of external rhythms (daylight, moonlight and tides) and internal circadian clocks as well. A generic "black-box" modeling approach has been used, and NARX structure of the model has been selected. The obtained model has been successfully applied for monitoring of water quality in the area of installation of the population of oysters, the obtained results are very encouraging.

In future works more sophisticated model structures have to be analyzed, and more hypothesis on the nature of circadian rhythms in bivalves have to be verified. The vector of inputs in the model can be extended by including the water temperature and the current cloudiness.

Fig. 1 .

 1 Fig. 1. [12] Synoptic representation of the system, from field to laboratory: (1) Oyster equipped with two electrodes; (2) connecting cables from oysters to first level electronic card; (3) first level electronic card in a waterproof case (immersed); (4) electrical connection between the first and the second electronic cards; (5) second electronic card (emerged); (6) GPRS antennae; (7) GPRS and Internet connection; (8) Marine Station of Arcachon; (9) daily update on Internet

Fig. 2 .Fig. 3 .

 23 Fig. 2. Hourly averaged valve openings for three oysters

Fig. 4 .

 4 Fig. 4. Averaged hourly valve opening for the population

Fig. 5 .

 5 Fig. 5. Structure scheme of the generic model (1),[START_REF] Kröger | Sensing the sea[END_REF]

Fig. 6 .

 6 Fig. 6. The results of the model numerical verification

Fig. 7 .

 7 Fig. 7. Residual for water quality monitoring

http://molluscan-eye.epoc.u-bordeaux1.fr/

http://www.shom.fr/