Cédric Gouy-Pailler
email: cedric.gouypailler@gipsa-lab.inpg.fr

Reza Sameni

Marco Congedo

Christian Jutten

Iterative Subspace Decomposition for Ocular Artifact Removal from EEG Recordings

In this study, we present a method to remove ocular artifacts from electroencephalographic (EEG) recordings. This method is based on the detection of the EOG activation periods from a reference EOG channel, definition of covariance matrices containing the nonstationary information of the EOG, and applying generalized eigenvalue decomposition (GEVD) onto these matrices to rank the components in order of resemblance with the EOG. An iterative procedure is further proposed to remove the EOG components in a deflation fashion.

Introduction

Electroencephalography (EEG) is a widely used technique for analyzing and interpreting human cerebral activity. EEG signals are usually interpreted by means of spectral and topographical measures that reflect global activity of the brain network. However, EEG measures are always contaminated by non-cerebral signals, which may disturb the interpretation of the brain activity. This issue has become a recurrent problem, for example in Brain-Computer Interfaces (BCI), where it has been proved to decrease classification error rates [START_REF] Thulasidas | Effect of ocular artifact removal in brain computer interface accuracy[END_REF].

Ocular artifacts generally occur during blinking or saccades of the eye, and are featured by high amplitude transient artifacts that defect the EEG. They are best recorded by an Electrooculogram (EOG) electrode or a pair of electrodes located close to the eyes. The high amplitude peaks are not seen on all channels, but mainly (and almost exclusively) on the fronto-paretal channels in combination with the occipital electrodes. These peaks are considered as one of the most considerable artifacts in EEG studies [START_REF] Koles | The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG[END_REF].

A common way of removing these artifacts is to apply independent component analysis (ICA) on multichannel EEG recordings and to remove the components which show a maximal correlation with a reference EOG channel [START_REF] Jung | Removing electroencephalographic artifacts by blind source separation[END_REF][START_REF] Ille | Artifact correction of the ongoing EEG using spatial filters based on artifact and brain signal topographies[END_REF]. However, it is not always possible to associate the components extracted by ICA to the EOG in an automatic and unsupervised manner. Moreover, EEG recordings can be rather noisy, and since ICA is based on a measure of independence (and not a measure of signal "cleanness"), the noise in the input channels can be amplified in the output, which again makes the detection of the true EOG component rather difficult. Lastly, most ICA methods are blind to Gaussian noise and as a consequence spread the Gaussian noise among the components.

Another common way of removing EOG artifacts from the EEG is to use a subtraction-based approach [START_REF] Romero | A comparative study of automatic techniques for ocular artifact reduction in spontaneous EEG signals based on clinical target variables: A simulation case[END_REF]. Here, the idea is to use the clean EOG recordings to remove ocular signals from EEG by a simple subtraction of a scaled EOG. Yet there is no evidence that EOG recordings are free of EEG. Thus by subtracting EOG we can also remove EEG signals of interest.

In recent works we have shown the applicability of Generalized Eigenvalue Decomposition (GEVD) for separating pseudo-periodic maternal ECG from fetal ECG signals recorded from the maternal abdomen [START_REF] Sameni | Multichannel electrocardiogram decomposition using periodic component analysis[END_REF]. In that work, one of the advantages of GEVD over other source separation techniques was the ability of ranking the extracted components in order of periodicity, which provided a means of automatic and unsupervised ECG decomposition and filtering. In this work, we develop a similar idea based on GEVD for the automatic detection and removal of EOG artifacts from multichannel EEG recordings.

The remainder of this paper is organized as follows: in section 2, we present a general nonlinear framework to decompose signals into independent subspaces using GEVD; the results of this method are presented in section 3 over simulated and real signals. The last section is devoted to conclusion and perspectives.

Method

Linear Transform

Suppose that we have an array of N EEG channels x(t), and a reference EOG channel denoted by EOG(t). Due to the spiky nature of the EOG, it is possible to detect the onset and offset times of the EOG artifact from the reference EOG channel. To do so, we define E(t), the averaged power of the EOG signal (or, alternatively the variance) within a sliding window of length w around t, as

E(t) . = 1 w w/2 τ =-w/2 EOG(t -τ ) 2 . (1)
Using this definition, an EOG is detected whenever E(t) passes some predefined threshold th. The active periods of the EOG may therefore be defined as

t a . = {t|E(t) > th} . (2)
As we will note later, the procedure of finding the offsets and onsets of the EOG does not need to be perfect, and the results can be further improved in a recursive procedure.

We now seek linear transforms of the multichannel recordings x(t), that maximally resemble the EOG (in the sense that the power of the extracted signals are concentrated during the active time periods). Therefore, if we denote the linear mixture as y(t) = b T x(t), we can try to maximize the cost function

ζ(b) = E ta y 2 (t a ) E t {y 2 (t)} , (3)
where E t {•} represents averaging over t. Here, the idea is to find linear mixtures of the input signals, with a maximal energy during the EOG activation time t a , while minimizing the global component energy. Equation ( 3) may be rewritten as

ζ(b) = b T E ta x(t a )x(t a ) T b b T E t {x(t)x(t) T } b , (4)
which is in the form of the Rayleigh Quotient [START_REF] Strang | Linear Algebra and Its Applications[END_REF][START_REF] Golub | Matrix Computations[END_REF], and may be solved by the joint diagonalization of two covariance matrices: the covariance matrix of the EEG channels over the whole dataset, and the covariance matrix of the data during the active periods of the EOG, respectively defined as

C x . = E t {x(t)x(t) T } , (5)
A x . = E ta {x(t a )x(t a ) T } . (6)
The intuition behind this method is to achieve decorrelated components that are at the same time globally and locally decorrelated. We already have a sense about the global decorrelation, which is achieved by sphering C x . In addition, the diagonalization of A x assures that the achieved components are also locally decorrelated over the active EOG epochs, too. This assures that the later extracted components have no redundancy up to second order statistics. The matrix that jointly diagonalizes the matrix pair (A x , C x ) is in fact the solution to the GEVD problem written as

U A x U T = Λ U C x U T = I , ( 7
)
where Λ is a diagonal matrix containing the generalized eigenvalues on its diagonal in descending order, and U is the matrix containing the generalized eigenvectors on its columns 3 . The decomposed signals may now be found by

y(t) = U T x(t) , (8)
where the elements of y(t) correspond to a linear transformation of the original data x(t) that are ranked according to their resemblance with the EOG activation epochs. This means that y 1 (t) contains the most information regarding the EOG while y N (t) is the least resembling the EOG. Another interpretation for the method is that y 1 (t) is the most dominant component of the EOG while having the least contribution in the overall energy of the signals (corresponding to the smallest eigenvalue of C x ). On the other hand y N (t) is the main non-EOG component, with the least EOG contamination (corresponding to the smallest eigenvalue of A x , i.e., having the least contribution of the active EOG epochs).

Signal/Noise Separation

Up to now, the components have been ranked according to their resemblance with the EOG through a linear transformation. The next step is to remove the ocular artifacts from the most contaminated components using a linear or nonlinear transformation. The result of this noise removal, denoted by z .

= (z 1 , • • • , z N ) T , can be expressed as ∀i ∈ [1, ..., N ] z i (t) = y i (t) -f i [y i (t)] , (9)
The denoising transform must be carefully chosen to remove ocular artifacts from the most contaminated signals, while preserving the non-EOG components. A simple but efficient choice is

f : f i (u) = u , i ∈ [1, ..., M ] f i (u) = 0 , i ∈ [M + 1, ..., N ] . (10)
This transform removes the first M (M ≪ N ) components contaminated by ocular artifacts and keeps the remaining N -M components unchanged. Therefore, in order to find the EOG free components, we can eliminate the first few components of y(t) and transform back the rest of the components using the inverse of the matrix U . The number of eliminated components M depends on the number of expected dimensions of the EOG subspace. Note that due to the elimination of the first M components, the rank of the multichannel signals are reduced to N -M .

Iterative Improvements

So far, we have removed the most dominant EOG components through a combination of a linear transformation, denoising, and back-projection. The results may be further improved by repeating the upper mentioned method in a recursive procedure. To do so, we can use y 1 (t), the component which most resembles the EOG, to re-estimate the onset and offsets of the EOG and its activation epochs using a smaller energy threshold th, and to recalculate C x , A x , and the other steps of the algorithm in each iteration. By repeating this procedure in several iterations, a better estimate of the EOG will be achieved. This iterative extension of the algorithm is of special interest for the cases in which a good EOG reference is not available. Therefore, we start with a coarse EOG onset and offset estimate; but we improve this estimate in the next iterations.

Experiments

Simulated Data. In order to test the method, we generate artificial EEG contaminated by ocular artifacts. The signals are generated as follows:

x(t) = EEG(t) + β (k 1 • • • k N ) T EOG(t) , (11)
where ( 1) EEG(t) is generated using a multivariate autoregressive filter (Yule-Walker order p = 8, cf. [START_REF] Anderson | Multivariate autoregressive models for classification of spontaneous electroencephalographic signals during mental tasks[END_REF] for details), trained using N = 4 channels of artifactfree EEG data, (2) EOG(t) results from the convolution of a typical blink segment, extracted from real EOG data, and an ensemble of Dirac impulses having a Poisson distribution with parameter λ = 0.2 Hz, (3) β is a scale parameter to adjust the signal-to-noise ratio (SNR) and ( 4)

(k 1 , • • • , k N )
T is a vector of gain under a linear hypothesis between EOG and EEG, for all i, k i ∼ N (1, 0.3) . A multivariate autoregressive filter has been preferred over a multi-univariate autoregressive filter, because the multivariate filter naturally takes into account spatial correlations of the data and learns those correlations from the real data.

Our approach is compared with a classical ICA approach for EOG removal, based on FastICA [START_REF] Hyvärinen | Independent Component Analysis[END_REF]. The input of the algorithm consists of the N EEG channels x(t) and the EOG channel EOG(t). FastICA seeks the most independent N + 1 components. We then remove the component which is the most correlated to the EOG, by setting it to zero. Components are then back-projected onto the sensor space.

The method proposed in this paper and FastICA, abbreviated GEVD and ICA, respectively, are evaluated by comparing the first N components of the back-projected signals and the EEG generated by the multivariate autoregressive filter. The idea behind this evaluation procedure is that a perfect denoising would lead to perfect recovery of the multivariate autoregressive process EEG(t), free of any EOG contamination. We therefore define the error signals as a function of the denoised processes resulting from ICA and GEVD (x ICA and x GEVD , respectively) as

ǫ ICA (t) = x ICA -EEG(t) (12) ǫ GEVD (t) = x GEVD -EEG(t) , (13)
from which we can define the following performance index (in decibels):

Q = 10 log 1 N N i=1 V(ǫ GEVDi (t)) V(ǫ ICAi (t))
where V(•) denotes the variance operator. This index is built such that a positive value indicates that ICA outperforms GEVD, while a negative index indicates that GEVD outperforms ICA. Figure 1 shows the results obtained for simulated data. Two signal lengths are considered (50000 and 100000 points). Two hundred simulations are done for each box-plot. The SNR tuned by β, is varying between -10dB and 30dB. These values are chosen by considering that a blinking artifact can have an amplitude of up to 100 times stronger than the EEG. Such a condition would yield an SNR of between -10dB to 0dB, depending on the blinking rate.

In the GEVD method, only the first component given by GEVD is removed. This choice is due to the rather simple method used to generate the EOG and EEG mixture. Moreover, in this case the iterative EOG improvement method is not used, since the EOG onsets and offsets are already known. We highlight the fact that in our simulation scenario, the simulation condition for the ICA method is better than it is for real signals, since we provide it with the clean EOG channel, directly.

Our method based on GEVD is shown to clearly outperform a classical ICA approach for a large set of SNRs (t-test highly significant). This effect is much more significant for low SNRs, showing that our method is of particular interest for bad conditionings of the signals. Real Data. We also evaluate our algorithm using a BCI experiment dataset4 [START_REF] Brunner | BCI competion 2008 -Graz dataset A[END_REF]. The data consists of 22 EEG channels and three EOG channels. The interelectrode distance is about 3.5 cm. The signals are sampled with 250 Hz and bandpass filtered between 0.5 and 100 Hz. An additional 50 Hz notch filter is used to suppress power-line noise.

The back-projected signals after removing the first EOG component are depicted in Fig. 2. In this example, due to the quality of one of the EOG channel, no improvements are achieved by iterating the method. This figure only shows two of the 22 denoised EEG channels. The first one, Fz, is known to be highly contaminated by EOG artifacts because it is close to the eyes. Fig. 2 shows that the EOG component is removed from this channel and also that the EEG signal outside blink contamination is kept perfectly unchanged. Channel Cz, known to be less contaminated than Fz is also drawn (Fig. 2 middle). While not as clear as in the first presented channel, we can see that the most important part of the EOG contamination is removed from the channel.

To illustrate the interest of the iterative method, we then show that due to the high amplitude features of the EOG contamination, we can even use some EEG channels to evaluate the active periods of EOG. This is of particular interest when no EOG channel is provided. Such a situation is illustrated in Fig. 3 where we have started with a noisy EEG reference and improved this reference after two iterations of GEVD.

Conclusion

We presented an automatic method to remove ocular artifacts from EEG measurements. Our algorithm outperformed a classical approach based on independent component analysis of EEG data. This method is a special case of a more general framework designed to generally decompose multivariate signals into independent subspaces. This flexible framework allows to choose different criteria for determining the optimal linear transform. In this paper, the nonstationarity of the EOG was used; but in other applications other prior information such as the periodicity or spectral contrast may be used.

In future works, further improvements can be achieved by replacing our simple component nulling function by some wavelet-based denoising method applied on the first few channels. In that case, it is expected that we separate the EOG from the EEG, without reducing the rank of the EEG signals. 3. Iterative improvements of the estimation of EOG. We start from a noisy EEG channel (up). The active periods is then determined according to the energy of the noisy EEG channel. We then select the first component to iterate the method: the active periods of the (i + 1) th iteration are determined using the energy of the first component extracted in the i th iteration.

Fig. 1 .

 1 Fig.1. Method comparison: for different values of initial signal to noise ratios between EEG and EOG, we performed 200 simulations. The performance index is positive when ICA performs better than GEVD whereas it is negative otherwise.

Fig. 2 .

 2 Fig. 2. Denoising of two specific channels, Fz and Cz electrodes (left and right, respectively). From top to bottom, the signals are (1) initial EEG, (2) denoised EEG, (3) residuals, and (4) EOG recordings.

Fig.

  Fig.3. Iterative improvements of the estimation of EOG. We start from a noisy EEG channel (up). The active periods is then determined according to the energy of the noisy EEG channel. We then select the first component to iterate the method: the active periods of the (i + 1) th iteration are determined using the energy of the first component extracted in the i th iteration.

Note that, contrary to the eigenvalues of symmetric matrices that are mutually orthogonal, generalized eigenvectors, i.e., the columns of U , are not generally orthogonal to each other; but following[START_REF] Koles | The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG[END_REF] they are "Cx-orthogonal"[START_REF] Strang | Linear Algebra and Its Applications[END_REF], p. 344. It means that each column of U is orthogonal to every column of Cx.

BCI competition IV, dataset 2a, subject A05[START_REF] Brunner | BCI competion 2008 -Graz dataset A[END_REF]

⋆ This work was supported by Égide, the French Ministry of Defense (Délégation Générale pour l'Armement -DGA), the Open-ViBE project, and the European program of Cooperation in the Field of Scientific and Technical Research.


Yandika Sirgabsou
email: yandika.sirgabsou@laas.fr

Claude Baron

Laurent Pahun

Phillipe Esteban

Software fault propagation patterns for model-based safety assessment in autonomous cars

Keywords: MBSA, MBSE, Safety analysis, Embedded software, Fault propagation pattern

The development of driver assistance and autonomous driving systems for vehicles has started to revolutionize the transportation sector, promising comfort, and safety. While significant technological progress has already been made in this area, many challenges remain. Among these challenges, ensuring safety has become even more critical due to the increasing use of complex, communicating, and reconfigurable embedded software. Current solutions to address safety include the use of model-based approaches for safety analyses instead of the traditional document-based safety analysis that is both informal and inefficient when faced with complexity. To this end, and in the context of automotive embedded software, we propose to rely on the use of fault patterns to improve the construction of software models used to conduct safety analyses. This paper makes a methodological proposal that improves current practices in terms of facilitated model construction and reusability, and that has been validated on the study of an automotive software component.

I. Introduction

The rapid development of embedded systems has led to numerous innovations in various systems in our modern society such as autonomous vehicles and highly computerized systems in airplanes. The technological challenges related with complexity and societal needs for guaranteeing safety induced by this trend opened new avenues for research in systems engineering but also exacerbated existing problems as they relate to the use of critical software and its contribution to systems safety.

To cope with these issues, industrials developing safety critical systems are looking for new methods and tools for designing and sharing design ideas more efficiently while ensuring system safety as required by standards and regulations. In the past decades, most of systems and software engineering development processes relied on document-based methods that relied on informal design documents to convey design ideas and artefacts from one development stage to the other. The informal aspect of these practices (such as manual analysis based on informal documents that are subject to the interpretation of the safety analyst) makes them prone to errors and less efficient in regard of the complexity of today's systems architectures. Although they are still widely used, these methods are now being challenged and model-based are more and more favored. In this context, embedded systems manufacturers are turning towards model driven engineering as part of both their systems and software development as well as systems safety assessment. In Systems Engineering (SE), this had led to the adoption of MBSE (Model-Based Systems Engineering), a systems engineering practice aimed at describing both a problem (need) and its solution through models, concepts and languages [START_REF] Baron | Towards a continuous certification of safety-critical avionics software[END_REF]. Its adoption can now be considered a success story as we witness that more industrials developing safety critical systems are turning towards the MBSE approach in Systems and Software Engineering. Examples include Dassault with its integrated 3DS MBSE solution or SIEMENS that integrates MBSE within its Product Lifecycle management (PLM ) solution. In Systems safety, a similar trend has led to the development of Model based Safety Assessment (MBSA), a practice that enables the capture, through specific formalisms and languages, of a systems safety related model (that describes the failure behavior and unifies all the safety property of a system in a single model), on the basis of which different safety analyses can be made.

However, despite the discipline being an early pioneer in the use of models, the wide adoption of model-based approaches for safety assessment has remained embryonic. In automotive, the ISO 26262 standard [2], titled "Road vehicles -Functional safety", requires performing safety analyses not only at system level but also at software level, to ensure the safe behavior of the embedded software. Moreover, in the context of autonomous driving, embedded software assumes various critical safety functionalities. Unfortunately, today, the currents practices in safety analyses do not focus enough on embedded software even though the software implements the logic of some of the critical safety mechanisms. As a result, in the software context, safety analyses are either not performed or if performed, only done through traditional document-oriented approaches. Therefore, there is a need for more focused and rigorous methods for safety analyses at software level.

This paper aims to propose a methodology, specially aimed at improving the practice of automotive embedded software safety analysis thought the use of fault patterns within the MBSA approach. In the next section, a state of the art of current MBSA approaches is presented. In section 3, a methodological proposal based on the use of fault pattern is made and applied to a case study in section 4. The results are discussed in section 5 and a conclusion in made in section 6.

II. State of the art

In systems safety, safety assessment has been dominated by document-centric methods and processes since the 60s. Classical methods such as Fault Tree Analysis (FTA) and Failure Modes and Effects Analysis (FMEA) have been complementarily used by experts in the industry for systems safety analyses. In these classical approaches, safety analyses are manually performed based on paper-style artefacts such as systems drawings or spreadsheets that are found in design documents. While these practices use a sort of well-defined semantics (such as Boolean operators and FTA symbols), their representation is often very far from the systems they describe.

With the introduction of Model-Based Safety Analysis (MBSA) starting around early 90s, through the earliest model-oriented safety analysis techniques such as FPTN [START_REF] Fenelon | An integrated tool set for software safety analysis[END_REF], Figaro [START_REF] Bouissou | Knowledge Modelling and Reliability Processing: Presentation of the Figaro Language and Associated Tools[END_REF] or the AltaRica [START_REF] Arnold | The AltaRica Formalism for Describing Concurrent Systems[END_REF] language, the focus has shifted to model-based approaches that no more base safety analyses on paper-style documents, but on a formal model of the system under design. This move has led, nowadays to an academic trend that seeks to address the interdisciplinarity and the consistency of the use of broad models through MBSE and MBSA throughout Systems Engineering and Systems Safety. In this trend one possibility is to directly associate the MBSA safety models with the MBSE system models. This allows conducting safety analysis directly on an extended version of the MBSE model. Examples of methods that are based on this extended approach include the approach included in the xSAP/Nu-SVM [START_REF] Bozzano | Improving System Reliability via Model Checking: The FSAP/NuSMV-SA Safety Analysis Platform[END_REF] [7] safety analysis tool. Another example is the Hierarchically Performed Hazard Origin Studies [START_REF] Papadopoulos | Hierarchically Performed Hazard Origin and Propagation Studies[END_REF] based on extended SIMULINK models [START_REF] Jiang | Model-based safety analyses of embedded system using stateflow[END_REF]. Although a clear advantage of the extended approach resides in the consistency (between system design and safety) it enables through the use of a shared system and safety model, it is argued that basing safety analysis on an extended model can lead to false assumptions and thus leading to hidden safety flaws (despite this claim not being shared by all systems engineering and safety practitioners). Furthermore, the safety analysis resulting from such model could be difficult to exploit because of the complexity of the extended model. In practice, safety assessment models are primarily used for mathematical or probabilistic calculations such as minimal cuts or monté carlo simulations. The more complex the model, the more computing power is required to effectively perform these calculations. Furthermore, the result of the safety analysis can be difficult to exploit due to the source models being blurred and overloaded (making them lose their ability to support a seamless communication) while the generated formal models are incomplete and uselessly complex [START_REF] Rauzy | Foundations for model-based systems engineering and model-based safety assessment[END_REF].

As an alternative to the described first approach, a more appealing trend in academic research is to build a separate MBSA model that needs to be kept consistent with the MBSE model through additional measures such as model synchronization [START_REF] Legendre | Toward Model Synchronization Between Safety Analysis and System Architecture Design in Industrial Contexts[END_REF]. This last approach is often based on dedicated safety modeling languages such as AltaRica [START_REF] Point | AltaRica: Constraint automata as a description language[END_REF], Figaro [START_REF] Bouissou | Knowledge Modelling and Reliability Processing: Presentation of the Figaro Language and Associated Tools[END_REF] or SAML [START_REF] Gudemann | A Framework for Qualitative and Quantitative Formal Model-Based Safety Analysis[END_REF]. These safety focused languages allow unambiguous representation of systems for the needs of safety analysis using well defined syntax and semantics. However, early experiments feedback suggests that, in the case of the dedicated model approach, the MBSA model construction can be challenging for safety analysts especially for complex systems. In such case, it is imperative to find the right level of details in modeling for fear of having a model too complex that may well be at the limit useless or irrelevant. Furthermore, some analysts don't necessarily see the advantage and gain in time of modeling a separate dysfunctional architecture for safety analyses. However, as argued by Rauzy and Haskins in [START_REF] Rauzy | Foundations for model-based systems engineering and model-based safety assessment[END_REF], systems description and safety analysis models are different by nature and efforts to unify such models in one single super model remains unrealistic. Based on this argument, the right direction is to keep a consistent separate MBSA model but make its construction easier. Nevertheless, the separate model construction can be challenging for large systems (what details, what modeling strategy etc.).

To ease the construction of dysfunctional models to conduct safety assessment, efforts have been done mostly using generic libraries of system elements that exhibit some safety properties. An example is the Safety Architecture Pattern(SAP) approach proposed by Kheren in [START_REF] Kehren | Architecture Patterns for Safe Design[END_REF]. In this approach, a library of SAP (components that highlights useful system's attributes from a safety point of view) are developed and coded using the AltaRica language. The generic library is then reused to easily prototype safety-oriented systems architectures that can be reused to perform safety analysis using tools such as the OCAS Workshop [START_REF] Bieber | Safety Assessment with Altarica[END_REF]. Nevertheless, the proposed approaches are mostly aeronautic systems oriented. The developed libraries are often dedicated to avionics systems (such as pumps, electrical motors, valves, or control units). While these libraries can be used for modeling physical systems at system level in automobile, they are less suitable for modeling dedicated safety architectures of the embedded software. Moreover, although there are ongoing works that aim to apply the MBSA approach to automobile at system level, less focus is being put on embedded software. However, in the automotive context, ISO 26262 recommends conducting safety analyses not only at system level, but also at software level [2]. In this context, if the model-based approach is to be used for safety analysis at software architecture level as the commended by ISO26262, safety analyses can be made easier if patterns or libraries of safety related components (such as safety mechanisms used in software) can be developed drawing from the same principles as those described in the case of systems SAP-oriented approach.

III. Methodological proposal

The goal of the methodological proposal is to construct a fault library of reusable software safety mechanisms that are commonly found in safety related software components, drawing form the SAP library-oriented approach described in the state of the art and given the need for improving safety analysis practices at software level in the automotive context. Our choice has been to use the dedicated model approach coupled with dedicated languages such as AltaRica as described in [START_REF] Sirgabsou | Investigating the use of a modelbased approach to assess automotive embedded software safety[END_REF]. However, as stated earlier, building a dysfunctional model can be challenging especially for complex systems. To address this concern, our solution has been to focus on selecting and including in the dysfunctional model only components that are safety-related as described in our previous work [START_REF] Sirgabsou | L'ingénierie dirigée par les modèles pour assurer la sécurité des logiciels embarqués en automobile[END_REF]. Even then, the failure propagation logic of these components must be manually written by the modeler (which can be time consuming for large systems).To address this concern, our first hypothesis is that making the MBSA model construction easier can both benefit its adoption by companies and improve the quality of safety analysis. A second hypothesis is that limiting the MBSA model to safety related components is sufficient to carry out meaningful safety analysis and can improve both efficiency and the quality of safety analysis. Therefore, the position of this paper is to make less painful for safety analysist the construction of MBSA models by proposing a set of predefined reusable libraries of software fault models to ease the dysfunctional model building process and improve the quality and consistency of the analyses. Its scope is limited to the context of automotive safety analysis at software level consistently with ISO 26262. Our methodology proceeds in 3 basic steps. The first one is to identify the safety mechanisms and their related software failure modes. The second step is to write the failure propagation logic through the safety mechanisms based on their functions and the identified failure modes; and to store the components in a library. The last step consists in reusing the elements stored in the library to build a dysfunctional model and conduct safety analyses.

In the first step we proceed by identifying the software failure modes and associated safety mechanisms found in automotive software architectures. In automobile, these failure modes and associated safety mechanisms are well defined by IO 26262 [13, Annex E, Annex D]. They are further detailed by two annexes in AUTOSAR [START_REF]Overview of Functional Safety Measures in AUTOSAR[END_REF] [START_REF]Safety Use Case Example[END_REF]. These failure modes are clustered into 4 categories. They include "data integrity, initialization & configuration data" "data exchange", "timing & control flow" and "data processing". The "data Integrity" category summarizes all failure modes related to corrupted memory or initialization and configuration data related to the corruption of software data at one memory address such as the corruption of memory content, memory partitioning fault or memory access fault. The "Data Exchange" category covers failures related to data transmission between sender and receiver such as between different ECUs (Electronic Control Units) or software components. The "Timing and Control Flow" category covers failure modes related to the timing of execution and scheduling. To model these fault categories, we start by describing the software components' behavior through a generic abstract template using states and transitions. An example of such a template consisting of a software component with generic failure modes is provided in in Figure 1. It shows 4 states (Inactive, Nominal, Erroneous and Failed), representing the execution status of the software component linked by several possible transitions (Minor fault, Major fault, Recovery etc.). The inactive state is the idle or initial state preceding the initialization where the software component is not solicited or does not provide its function. From this state, the 3 outgoing transitions labeled "Activation", "Major fault" and "Minor fault" will lead the component to the "Nominal", "Failed", or "Erroneous" states respectively. In the nominal state, the software component executes and delivers its intended function. From this state, the component can revert to the Inactive state as indicated by the transition label "Deactivation" or move to the "Failed" or "Erroneous" states if a major or minor fault occurs as indicated by the transitions. In the erroneous state, the software component provides erroneous results while, in the "Failed" state, it fails to provide its intended function. From the generic pattern, more specific patterns related to specific fault categories as described by ISO 26262 can be easily derived.

As an illustrative example, let us consider a piece of software that reads some critical data from memory, performs some critical calculations, and stores the result back to memory. Based on its function, the failure modes of this software component could include memory access related faults (such as read and write errors) as well control flow and timing related failure modes (such as execution failure, or untimely execution). Depending on the exact safety requirement, safety mechanisms to prevent the failure of this software component's function could include a protection against unwanted writing and a watchdog timer. Based on this information, the elements that will be necessary to model the safety related behavior of this software component, are the failure mode related to the memory access and the two safety mechanisms that will be modeled though states and transitions within this software component.

Figure 1. Failure modes of a generic software component

The second step focuses on writing failure propagation logic of the identified safety mechanisms knowing their function and considering the associated failure modes identified in the first step. To do so, we first need to study the safety mechanism and identify their basic behavior, what their inputs are, and the result they produce in normal or faulty execution. Once done, we can write the failure propagation logic of the identified safety mechanisms. To this end we use Failure Truth Tables (FTT) that we introduced in our early work [START_REF] Sirgabsou | Investigating the use of a modelbased approach to assess automotive embedded software safety[END_REF]. FTTs are dysfunctional failure propagation tables consisting of discrete input and output variables whose possible values are defined depending of the failure behavior of the components function. FTTs can be used to capture the dysfunctional logic of a function based on its inputs. Depending on the expression of the safety requirements that the identified mechanisms are to fulfill, logical operators such as "or, and" and program control structures such as "if-then-else" can be useful to express the function or the combination of several mechanisms within a software component. Therefore, we also need to write and add to the library, the failure propagation logic through theses operators and control structures. This will result in a library of safety mechanism, operators and control structures that will be used in the last step to construct the dysfunctional model of a system.

The last step consists in using the elements stored in the library to construct the dysfunctional model. This step requires having a tool that offers library support. Using the elements stored in the library in an appropriate MBSA tool, one can easily model the dysfunctional architecture of a systems through drag and drop. This is possible since the safety mechanisms self-contain their propagation logics as well as the associated fault modes. Therefore, there is no need to write the failure propagation logic code in this step.

IV. Implementation in SimfiaNeo and case study

The objective of the case study was to first develop a library of safety mechanisms and used it to model the dysfunctional architectures. To this end, we are using SimfiaNeo [START_REF] Machin | SimfiaNeo, Complex Systems, yet Simple Safety[END_REF] , an MBSA modeling tool based on the AltaRica language and developed by APSYS-Airbus. It offers a graphical modeling interface based on Eclipse and implements the dataflow version of the AltaRica language. SimfiaNeo allows to graphically build the AltaRica model of a system and to directly perform various safety analyses based on minimal cuts or in the form fault trees and FMEAs directly from the AltaRica model. Furthermore, the tool in its latest version offers library features that make it suitable for our proposal. It allows the modeling, storing and instantiation of custom components in a library.

In the context of Advanced Driver Assistance Systems (ADAS) and Autonomous Driving (AD), we aim to apply the proposed methodology on a practical case study, the longitudinal control software component. The longitudinal control is a function of the ADAS technology. It is a software component whose purpose is to ensure speed and braking control in autonomous driving mode. It is built around the ACC (Adaptive Cruise Control), a speed and distance control system that calculates how fast the vehicle can travel while remaining in a safe situation with respect to certain predefined events (turns, traffic jams, stop signs, etc.). To ensure its safe activation, the longitudinal control relies on an internal monitor (a subcomponent called supervisor that manages its states) and a failsafe controller (a subcomponent that places the longitudinal control in some predefined safe states when some faults occur). Safety requirements associated with the longitudinal control are documented in a safety concept through safety goals that are declined to Unwanted Software Events (UWSE) at the software level. In this case study, we focus on one single UWSE related to the occurrence of an unintentional acceleration above the permissible acceleration limit (defined by ISO 22179) during travel (v ≠ 0 km/h) entailing longitudinal control. Other constraints also exist (e.g., the user must be able to deactivate the ACC at all times).

Pattern prototyping and dysfunctional model construction with SimfiaNeo

Using our described approach, we modeled the basic structures of software components using the previously described fault categories and safety mechanisms. First, we declared the necessary AltaRica domains in SimfiaNeo. To this end, we declare 5 domains with different literals encompassing data integrity, Data exchange, Safety mechanism, generic state, and generic data. Using these 5 domains we modeled bricks of components representing the elements of a software safety architecture including elements such as generic software components, generic safety mechanisms, Boolean operators, and program control structures.

Description the fault patterns

In automobile, software fault modes and associated safety mechanisms are clustered into 4 categories. Including include "data integrity, initialization & configuration data" "data exchange", "timing & control flow" and "data processing". In this subsection we present a fault pattern related to the data exchange as an illustrative example and a generic behavioral pattern for safety mechanisms.

As described by ISO 26262 in the software scope, the "Data Exchange" category covers failure modes related data transmission between sender and receiver such as between different ECUs or software components. The pattern presented on Figure 2 shows how the execution of a receiving software component subjected to this category can be affected. Like the pattern shown in Figure 1, it has 4 states: "Init," "Nominal", "Erroneous", and "Failed". From the initial state (blue state on Figure 2), the function will either move to the "Nominal" or "Erroneous" states as indicated by the outgoing transitions depending on a successful or unsuccessful data transmission initialization. In the nominal sate, the software component executes normally and fulfils its function. If a data transmission initialization fault has led the function to the "Erroneous" state, an execution of a safety mechanism can bring the function to nominal if successful or to the "Failed" state (red state on Figure 2) if unsuccessful. The function can also move to from the nominal state to the "Erroneous" state with the occurrence of inconsistencies of the transmitted data (such as corruption, incorrect data value, out of range data values or incorrect sequence of data). In the "Failed" state, the software component fails to provide the expected function due to missing or loss of transmitted data or due to the safety mechanism failure to recover from the "Erroneous" state. As it can be seen in Figure 2, this pattern is built on the generic pattern presented in Figure 1. However, it differs by the specificity of its failure modes expressed in the transitions that are specific to "Data Exchange" fault category. Based on this pattern, another related pattern was derived to cover the specificity of other data exchange related faults such as delayed data transmission that will cause the function's execution to be delayed. It such case, the "Erroneous" state was further split into several states depending on the specificity of the software component.

Figure 2. Data exchange fault pattern

The next pattern aims to capture the behavior of a generic software safety mechanism. The 4 states (Nominal Inactive, Nominal Active, Misleading and Failed) represent the execution state of the safety mechanism. In the "Nominal inactive" state, the safety mechanism is in its nominal execution state with no fault detected. When a fault is detected, it moves to the "Nominal active" state. In this state, the safety mechanism is successful in reacting and correcting the effect of the fault. From this states, an erroneous or failed reaction will lead the safety mechanism to "Misleading" or "Failed" states respectively.

Figure 3. Generic safety mechanism fault pattern

The behavior of the safety mechanism is completed with writing failure propagation logic knowing their function and considering the associated fault modes identified earlier. To do so, we first need to study the safety mechanism's function and identify their basic behavior, what their inputs are, and the result they produce in normal or faulty execution. Once done, we can write the failure propagation logic of the identified safety mechanisms. Depending on the expression of the safety requirements that the identified mechanisms are to fulfill, logical operators such as "or, and" and program control structures such as "if-then-else" can be useful to express the function or the combination of several mechanisms within a software component. While these operators are already part of the AltaRica semantics and can be expressed in assertions, modeling them in virtual bricks allows to graphically associate the components without rewriting the propagation logic. Therefore, we also need to write and add to the library, the failure propagation logic through theses operators and control structures. This will result in a library of safety mechanism, operators and control structures that will be use in the last step to construct the dysfunctional model of a system.

Modeling

Having opted for an approach based on a dedicated model, we must first identify the information necessary for its construction, starting with the software architectural design documents and the Technical Safety Concept (TSC) resulting from the system level safety assessments. The TSC is an aggregation of safety requirement specifications (often in textual and tabular format) from the system, as well as their allocations to hardware and software components and associated information (text, diagrams or sketches), which justify that safety measures and mechanisms are in place. Based on the TSC and the definition of the items, we can identify the safety-relevant components and interfaces to model, as well as the requirements and safety mechanisms to evaluate in the context of the dysfunctional architecture. In this way, we can represent in the dysfunctional architecture only those components that impact the safety goals, which are high-level safety requirements resulting from the preliminary risk analysis at the vehicle level (see ISO 26262-1 3.139). This will also avoid overloading the dysfunctional model with elements unnecessary for safety.

We used the readily available model bricks to model the patterns and their states in the tool, assigning the previously created domains to them. An overview of the modeled system along with the fault patterns is presented in Figure 4. Depending on the function of the pattern, different domains were used. Through the creation of AltaRica events in SimfiaNeo, we then modeled the transition firing conditions using the events identified in the state machines previously presented in the methodological proposal section. An AltaRica transition is characterized by a guard (condition to fulfill before triggering the effect), an effect (action resulting from the state change), and potentially a distribution (exponential, Dirac etc.) associated to the event. For each event, the SimfiaNeo tool allows to specify a probability that will be used during calculations. However, given in the context of software faults, such values are irrelevant and a probability of 1 was used instead. For each pattern, we wrote fault propagation logic linking the corresponding inputs (if they exist), internal states and outputs. We described the states of the inputs and outputs using the AltaRica domain that we named "Generic Data" and that incorporates four states: Nominal, Lost, Delayed, Erroneous. Using these states, we were then able to model the dysfunctional information flows between components trough AltaRica expressions. Using the elements stored in the library in an appropriate MBSA tool, one can easily model the dysfunction architecture of a systems through drag and drop. This is possible since the safety mechanisms self-contain their propagation logics as well as the associated fault modes. Therefore, there is no need to write the failure propagation logic code in this step. We constructed the model presented on Figure 4 using the elements stored in the library as shown by the library icon in some of the components (e.g., component identified as ETH at the bottom left of the models). As an examples the components labeled CAN and ETH in Figure 4 were both modeled thorough an instantiation of the data exchange fault pattern described earlier as well as some internal components in the "vehicle-status-input" component that receive these data. Similarly, the "memory" component shown on Figure 4 as well as some subcomponents in the "longitudinal controller" component that read data from the memory were modeled though the instantiation of the data integrity fault pattern. Through these reusable libraries, we are able to model the dysfunctional behavior and failure propagation without having to manually rewrite their AltaRica code.

Safety Analyses

Using the SimfiaNeo Tool, we performed various safety analyses including step by step simulations, minimal cuts, FTA and FMEA based on the dysfunctional model constructed using the previously developed fault patterns. Having completed the modeling of the longitudinal control software component and the components that interact with it, the objective was to perform safety analyses from the dysfunctional model. For this purpose, we set up AltaRica observers on the outputs we were interested in (as shown in Figure 4). An AltaRica observer is an indicator that can be associated with a failure condition or feared event that we wish to capture. For example, let us consider, the UWSE that we have chosen for our case study related to an "unintended Acceleration > ISO 22179 acceleration limit while travelling (v ≠ 0 km/h) requested by the longitudinal control feature". In our model, we identified that the acceleration target and request in the speed controller subcomponent (Speed-Ctrl) are limited to 0.2G until vehicle speed vehicle is above 10 km/h. We also identified that the final value of the acceleration target is transmitted to the engine through the engine management command 'PWTWheelTorqueCmd' (Powertrain Wheel Torque Command). Thus, any erroneous value of 'PWTWheelTorqueCmd' can result in the violation of the safety goal and the occurrence of the UWSE. Therefore, the observers predicate to capture the occurrence of this UWSE can simply be 'PWTWheelTorqueCmd =Erroneous'. Having added the expression of the observer to the model, the objective was to verify, by means of the simulation, FMEAs, minimal cuts and fault trees, whether we could Domains Fault patterns Observer determine the events or failures that could lead to the transmission of this erroneous command that can result in the violation of the chosen UWSE.

Simulations

Having run a series of simulations, we observed the propagation of failures to observers through the visualization of components in different colors (red: in the presence of a failure; orange: in error state; green: in nominal state), as shown in Figure 5. However, hierarchical components (consisting of blocks containing several subcomponents) appeared without color during the simulation. This step-although not overly formal-allows the model to be verified as it is being built. The analyst can then use it to quickly evaluate the dysfunctional architecture by visualizing how all the components and observers react to the presence of one or more failures at specific locations. The simulation can be used to confirm and demonstrate (for communication purposes) the feared scenarios identified with the classic methods (FMEA and fault trees) that we will discuss in the following subsections.

Failure Mode and Effects Analysis

We used SimfiaNeo to generate the FMEA tables from the dysfunctional model we had built. The FMEAs list all the events resulting in the violation of a safety goal (or of a created observer), doing so for each component of the model. Figure 6 shows an excerpt from an FMEA, containing a certain number of elements typically found in these tables. The first column (Event) lists the events causing the violation. In the following columns, we can find the Local Effect (effect of the event on the output of the initial component), the Intermediate Effect (effect of the event on all intermediate components between the initial component and the final observer), and the Final Effect (effect on the output of the model; in here, the effect on the observers described previously). For instance, in relation to our chosen UWSE, the excerpt from Figure 6 shows how faults in the vehicle status input component related to vehicle speed can affect other components and the final observer. SimfiaNeo can export this document as an Excel spreadsheet, allowing for better data processing and sharing. This is an important asset of the tool, considering that MBSA tools are not necessarily used by many but every engineer manipulates Excel files. Note, however, that Figure 6 represents only a very small excerpt from the initial FMEA, which has more than 18,000 lines. Therefore, if the goal is to obtain usable details or to make a synthesis, this representation is not ideal.

To analyses the usefulness of this FMEA, a comparison with a manually performed FMEA would have been interesting. However, in the context of current practice in our case, there are no software FMEAs performed using Figure 5. Simulation the traditional approach. In contrast to fault trees that focus on one feared event, FMEAs are systematic and constitute a great way of showing that all failure modes have been accounted for within the system. This remains a difficult task for the safety analyst especially in the software context where failure modes can be plethoric. Despite the absence of a comparison with a manually performed FMEA, we argue that our approach allows performing this type of analysis that is otherwise difficult to perform manually.

Minimal cut set and fault tree

The generation of the minimal cuts is achieved through the configuration of the cut set and sequence calculation engine for a given observer. Thus, for an observer we can choose the maximum order of the cuts, the filter type (minimal cut or minimal sequence) and the generation type (combination, permutation or stochastic) which will be used during the cut set or sequence calculations. The maximum order corresponds to the maximum number of primary events in a sequence. To choose the maximum order, we experimented with values ranging from 2 to 5. We observed that SimfiaNeo load on the processor and memory consumption remained relatively constant (respectively close to 30% and 700 Megabyte) regardless of the maximum order value, while the computation time increased exponentially from under 2 minutes for order 2 to 7 hours for order 5. Meanwhile the maximum order in the resulting minimal cut set remained equal to 3 even if we consider sequences of size 4 or 5. We chose order 3 for our case study-the higher the order, the longer the generation of the cut will take. An order of 3 is therefore a good tradeoff between computation time and accuracy. The choice of the filter type is also important; we have chosen the "minimal cuts" option since it makes fault tree generation possible. Lastly, the choice of the generation type specifies the combinatorial or stochastic sequence (based on random simulations) used during the generation of the cut.

Figure 7. Minimal cut

As an example, we considered the chosen UWSE linked to the transmission of an erroneous torque command to the engine. For this UWSE, we generated a minimal cut by choosing "order 3" as value of the maximum order, the "minimal cut" filter and "permutation" as the generation type. The generated minimal cut is shown in Figure 7. It shows combinations (of order 1, 2) of basic events that could cause the specified UWSE, as well as their associated probabilities (added by default). We can see that the cut highlights the events and the hierarchical components, enabling traceability of the components at high level (as shown in Figure 7). For dysfunctional models where several subcomponents have identical nomenclature, this traceability allows to clearly identify the origin of each event.

During the execution of these calculations, however, several compilation errors occurred, some of them due to the presence of loops in the failure propagation chain. This is a known issue related to the dataflow version of the AltaRica language. To solve this problem, we modified the assertions of the failure propagation involved in these loops. In the case of a redundant evaluation of a variable (where one of the assertion is part of the loop), removing the redundant evaluations of the involved variable allowed to break the loop. For this purpose, we considered a loop and identified the self-dependent variables in the chain of assertions constituting this loop. One possibility was to remove this variable from one of the assertions if it was already considered in another assertion. If this was not possible without modifying the validity of the assertion, the second possibility was to remove the dependency link and successively assign to the state variable all possible values and perform the calculations with each scenario. In the latter case, it was necessary to manually change the value of the variable to include the scenarios which were excluded by assigning it a fixed value. In both cases, the dysfunctional logic of the assertion remains valid.

For a defined feared event, SimfiaNeo allows the generation of FTAs from the equivalent minimal cut. Through its tree structure and logical combinations, FTAs illustrate how basic events (located at the bottom of the tree) can lead to the feared event (at the top of the tree). In other words, FTAs highlight the causal chain between the basic events at the component level (at the bottom of the tree) and the high-level feared event (at the top of the tree) through a tree structure represented in graphic form. In SimfiaNeo, it was possible to generate an equivalent reduced fault trees from minimal cuts for a defined feared event. Nevertheless, we did not identify any added value through this generation as the minimum cuts in our opinion present the same information in a more concise and readable format. Furthermore, generating FTAs from minimal cuts can be considered counter intuitive as in practice safety engineers use the reverse process (they use FTAs to compute minimal cuts).

V. Discussion

The results obtained from the application to the case study shows that using fault patterns and the adequate tool, the dysfunctional model construction is made easy and safety analyses can benefit from this alternative new analysis method. The specificity of the case study demonstrates that the use of software-oriented fault patterns can benefit the application of MBSA in the automotive software safety context. Furthermore, tools such as SimfiaNeo relieve the safety expert of manual calculations while allowing him to concentrate on modeling. The benefits are reflected in terms of reusability of the models. Once the fault patterns are built, they can be reused to build the dysfunctional model and make it possible to conduct analyses with different parameters for many UWSEs based on the same model. In the traditional approach, the safety analyst spends much of their effort interpreting various design documents to manually construct classical model's safety models such as FTAs or FMEAs. The analyst would need to manually construct as many FTAs as there are feared events. If the design evolves, they will need to individually update all the fault trees and the FMEAs. Through our proposal, the dysfunctional model is easy to construct. If the system design evolves, the dysfunctional model can be updated, and updated safety analyses can be automatically derived. In addition, representing the behavior of the system without ambiguity is possible through the formal semantics of AltaRica, turning it into a possible candidate in a certification context. All these elements make this safety analysis method an interesting alternative that has the potential not only improve current practices but to contribute to the adoption of the model-based approach. Nevertheless, we noted some limitations on the method. One of the difficulties brought about by a dedicated dysfunctional model is maintaining its consistency with the design model when the latter evolves; this problem was already true when working on analyses based on fault trees and is not addressed by our proposal. Implementing additional measures is therefore necessary to guarantee consistency. Another limitation, related to the AltaRica dataflow, is the inability to natively manage loops; our solution was to modify some assertions in order to remove these loops. Finally, as the case study involves a system of reduced complexity, the scaling up of our methodological proposal is yet to be evaluated and we will need to perform a comparison with safety analyses results obtained through the traditional manual approach to fully evaluate the proposed approach.

VI. Conclusion

This paper made a methodological proposal based on fault patterns that can be used to build a dysfunctional model, and from which it is possible to derive classic safety models. Using the SimfiaNeo tool and the AltaRica language, the methodology was applied on a case study, building a dysfunctional model of a software from which we were able to generate FMEAs and minimal cuts. These results are encouraging and demonstrate the usefulness of patterns to facilitate MBSA model construction. More generally, they show that it is possible to apply an MBSA approach to evaluate software safety, especially in automotive applications. They also highlight the benefits of generating safety analyses from a dysfunctional model (time saving and reusability). Building on these results, the study must now continue to evaluate the complexity of the systems for which the methodology and the tooling can be reasonably applied. Since the proposal of this paper is based on a dedicated dysfunctional model, it will also be essential to supplement the method with a mechanism that ensures consistency between the design models and the safety models.

Figure 4 .

 4 Figure 4. Pattern prototyping and dysfunctional model construction with SimfiaNeo

Figure 6

 6 Figure 6. FMEA


A Bressy
email: adele.bressy@enpc.fr

P Martinache

C Paijens

D Tedoldi

B Bonnaud

R Moilleron

Isothiazolinones: from their domestic uses to their release into the environment

Keywords: Biocide, combined sewer overflow, indoor dust, mass loads, river, wastewater

Isothiazolinones are ubiquitous in our daily life; they are used as preservatives in cosmetics, household products and in building materials. These biocidal substances may be emitted inside house, inducing a potential risk for health, and released into the environment through wastewater or runoff, thus inducing a potential impact on aquatic ecosystems. The overall objective of this communication was to explore the isothiazolinone contamination of the continuum House-Sewer-River and to assess the induced risks for health and the aquatic life. The experimental results showed the ubiquity of isothiazolinones in indoor dust, urban water discharges and in the Seine River. Annual mass loads discharged into the Seine River were estimated using Monte Carlo method at both the conurbation and annual scales, they showed that wastewater discharges constitute the major entry pathway.

CONTEXT AND OBJECTIVES

Despite intense use of isothiazolinones in household products [START_REF] De | Development and validation of a UPLC-MS/MS method for the quantification of isothiazolinones in the composition and emissions from consumer products[END_REF] and in construction materials [START_REF] Bollmann | Dynamics of biocide emissions from buildings in a suburban stormwater catchment -Concentrations, mass loads and emission processes[END_REF] and their potential harmful effects for humans and for aquatic receiving environments [START_REF] Kresmann | Ecotoxicological potential of the biocides terbutryn, octhilinone and methylisothiazolinone: Underestimated risk from biocidal pathways?[END_REF][START_REF] Kim | A Commonly Used Biocide 2-N-octyl-4-isothiazolin-3-oneInduces Blood-Brain Barrier Dysfunction via Cellular Thiol Modification and Mitochondrial Damage[END_REF], few studies have documented their presence along the continuum extending from the indoor environment to the river (Paijens et al. 2020a). Indeed, isothiazolinones are poorly integrated into monitoring programs. Most studies addressed small-scale emissions in runoff at building scale. This work focused on five isothiazolinones and aimed at (i) collecting data on their occurrence in indoor environment to assess Human exposure; (ii) quantifying their release into the receiving environment through wastewater and combined sewer overflows (CSOs); and (iii) assessing the contamination of surface waters and the associated risk.

METHODOLOGY

Sampling campaigns

This study focused on the continuum House-Sewer-River on the Paris conurbation example. Samples were collected at different levels of this continuum: (i) close to the use of isothiazolinones in the house, i.e., in indoor dust and greywater; (ii) at the level of urban discharges (wastewater treatment plant -WWTP, and CSOs); and (iii) in the aquatic receiving environment.

Sampling inside house: indoor dust and greywater

Indoor dust and greywater were sampled in order to estimate the indoor contamination by biocides linked to their uses in building material and daily products (households, cosmetics etc.). For that purpose, indoor dust was collected with a vacuum cleaner in 12 dwellings from the Paris conurbation. Greywater (cooking water, shower, floor washing, sink) was sampled 6 times in a building of collective catering.

Sampling inside the sewer system

Six sampling campaigns were conducted during dry weather in the influent and effluent of the Seine centre WWTP (900,000 inhabitant equivalent, physico-chemical lamellar unit and biofiltration). Eight sampling campaigns were performed during wet weather in the biggest Parisian CSO located in Clichy. For each campaign, composite samples (24 h for WWTP or discharge duration for CSO) were collected by automatic samplers equipped with Teflon® pipes and maintained at 4°C before analysis.

Sampling in the Seine River

The sampling campaigns targeted (i) the Marne and the Seine Rivers before their confluence upstream the Paris conurbation, and the Seine River downstream (10 sampling campaigns covering all hydrological characteristics of the year). At each sampling location, a boat was used to collect 5-L water samples in the middle of the river at a depth of approx. 1 m, thus avoiding the biases related to the uneven flow conditions and incomplete mixing close to the banks. These samples were composited from 5 sub samples kept in cleaned glass bottles and stored at 4°C.

Biocides analysis

Five isothiazolinones were selected on the basis of their intensive use in urban areas, known exposure of aquatic populations, ecotoxicity and technical feasibility (Paijens et al. 2020b): methylisothiazolinone (MIT), chloromethylisothiazolinone (CMIT), benzisothiazolinone (BIT), octylisothiazolinone (OIT), dichloro-octylisothiazolinone (DCOIT).

Dust samples were sieved at 100µm. Water samples were filtrated through glass fiber filters (0,7 µm, GF/F Whatman). Particles from water and from dust (~ 100 mg) were solvent extracted by microwaves. Dissolved fractions were solid phase extracted. All the extracted fractions were analyzed by ultra-performance liquid chromatography coupled with tandem mass spectrometry (Paijens et al. 2020b).

Data treatment

The method used for the assessment of biocide loads at the annual scale relies on Monte Carlo simulations and aimed at propagating the inherent variability and uncertainties of the measured variables [START_REF] Paijens | Urban pathways of biocides towards surface waters during dry and wet weathers: Assessment at the Paris conurbation scale[END_REF]. For WWTP and CSO results, theoretical distributions were first fitted to the experimental data, to characterize the temporal variability of concentrations. These distributions were then used to randomly reconstruct a 1-year series of daily/event-mean biocide concentrations, to which a Gaussian error term was added. These concentrations, multiplied by the daily flow rates measured during the study period, provided a point estimate of the annual load of each biocide. This procedure was repeated 100,000 times in order to achieve a statistical distribution of the biocide loads, thus enabling the associated uncertainty to be characterized.

RESULTS

Isothiazolinone emissions at home

Table 1 summarizes the frequencies of quantification and contents of isothiazolinones in indoor dust. MIT, BIT and OIT were quantified in all the 12 samples, CMIT in 9 samples and DCOIT in 2 samples out 12. These frequencies of quantification are consistent with two levels of concentrations. MIT, BIT and OIT were measured with respective median contents of: 490 ng/g, 850 ng/g and 260 ng/g. CMIT and DCOIT were measured with respective median contents of: 49ng/g and 34 ng/g. Greywater samples are currently being analysed.

Occurrence of isothiazolinones in the sewer system MIT, BIT and DCOIT were systematically quantified in WWTP influent and effluent samples as OIT in influent (Table 1), while CMIT was never quantified and OIT two times over six in effluent. Three levels of concentrations could be observed: (i) concentrations of MIT were generally higher than 100 ng/L, (ii) concentrations of BIT ranged between 20 and 60 ng/L, and (iii) concentrations of OIT and DCOIT were lower than 5 ng/L. In CSOs, isothiazolinones were quantified in almost all samples between <0.5 ng/L and 300 ng/L, except DCOIT exhibiting concentrations < 10 ng/L.

In order to discriminate the origins of biocides in CSO discharges (wastewater vs. stormwater), we looked for a trend between concentrations in CSO and proportions of stormwater (estimated from electrical conductivity). The concentrations of MIT and BIT decreased with increasing proportion of stormwater in the CSO, suggesting a dilution of wastewater by stormwater. These compounds are indeed widely used as preservatives in household products. No trends were observed for CMIT, OIT, DCOIT. These compounds thus appear to be contributed in equivalent proportions by wastewater and stormwater, consistent with their use in domestic products and as preservatives in building materials.

Occurrence of isothiazolinones in the Seine River

For the Seine River, the frequencies of quantification were found to be lower. MIT and BIT downstream were always detected, while BIT upstream, OIT and DCOIT were quantified between 5/10 and 7/10, and less than 3/10 for CMIT. Median total concentrations were generally below 10 ng/L, except for MIT (13-21 ng/L). Frequencies of quantification and median concentrations are higher downstream the Paris conurbation than upstream. The risk quotient (RQ ratio between the river concentration and the predicted no effect concentration (PNEC)) was calculated for each substance and each campaign. RQ was high (>1) for 70% of the samples for MIT (PNEC=50 ng/L). CMIT, OIT, DCOIT represented a low to medium risk (0.01 < RQ < 1). Annual biocide loads in WWTP effluents were found to range in the following order: DCOIT (1-10 kg/yr) < BIT (10-100 kg/yr) < MIT (> 100 kg/yr). The order appeared to be approximately the same for CSOs, but the mass loads differed by one or two order(s) of magnitude. Therefore, at the annual and city scales, WWTPs seem to constitute the major pathway of isothiazolinones into receiving waters because of greater discharged volumes. However, contrary to WWTP effluents, which are continuously discharged into the river over the year, CSOs are punctual releases and are generally discharged at higher flow rates than WWTP effluents. They can hence have a higher local and oneoff impact on aquatic organisms.

Table 1 .

 1 Quantification frequencies (F) and median concentrations (min/max) of isothiazolinones in the continuum House-Sewer-River

		Indoor dust		WWTP influent	WWTP effluent		CSO	Seine River upstream	Seine River downstream
			Median		Median		Median		Median		Median		Median
		F	conc. min/max	F	conc. min/max	F	conc. min/max	F	conc. min/max	F	conc. min/max	F	conc. min/max
			ng/g		ng/L		ng/L		ng/L		ng/L		ng/L
	MIT	12/12	490 59/1,360	6/6	620 (350 / 860)	6/6	150 39/350	8/8	70 9.8/290	9/10	13 <1/150	9/10	21 <1/170
	BIT	12/12	890 130/3,500	6/6	320 (210 / 660)	6/6	24 20/55	8/8	32 4.6/64	5/10	3.0 <1.2/6.8	9/10	7.0 <1.4/11
	CMIT	9/12	49 12/150	0/6	<13	0/6	<3.7	4/8	14 <9.5/160	3/10	<4.7	2/10	<14
	OIT	12/12	260 47/1,300	6/6	5.9 (1.1 / 9.0)	2/6	<1.8 <0.7/2.0	8/8	27 18/46	5/10	0.9 <0.3/1.2	7/10	1.1 <0.3/1.7
	DCOIT	2/12	34 31/36	5/6	4.5 (<0.6 / 10)	6/6	1.5 0.5/4.6	6/8	4.6 <0.5/6.3	7/10	1.1 <0.3/2.9	7/10	1.6 <0.4/43

AKNOWLEDGMENT

This work was supported by the Central Laboratory of the Police Prefecture (LCPP), the Biocid@Home research program (grant ANR-20-CE34-0001-01), the OPUR research program funded by the Seine-Normandy Water Agency, the City of Paris, the greater Paris sanitation authority (SIAAP), the Hauts-de-Seine Departmental Council, the Seine Saint-Denis Departmental Council and the Val-de-Marne Departmental Council. The authors want to thank the teams of the greater Paris sanitation authority (SIAAP) for their assistance for sampling.


Elsa Suberbielle

Daniel Gonzales-Dunia

Frédéric Pont
email: frederic.pont@toulouse.inserm.fr

Daniel Gonzalez-Dunia

pressure-resistant electrode

Keywords: HPLC, high performance liquid chromatography, MS, mass spectrometry, HPCF, high performance chromatofocusing, RPHPLC, reverse phase high performance liquid chromatography, NPS, non-porous silica, ACN, acetonitrile, TFA, tri-fluoro-acetic acid

High reproducibility of two-dimensional liquid chromatography using pH-driven fractionation with a

similarities, further stressing the need of careful pH adjustment and precise electrode calibration. Combined with good laboratory practice, 2D-LC using the PF2D system is a reproducible fractionation platform well suited for differential proteomics.

4

Introduction

One of the main challenges in proteomics centers on the ability to generate a reproducible fractionation of the protein samples. Indeed, the discovery of novel biomarkers, together with studies of disease pathogenesis, such as for cancer, rely on differential proteomics i.e., on the accurate comparison between control and pathological situations [START_REF] Monteoliva | Differential proteomics: an overview of gel and non-gel based approaches[END_REF]. Although the resolution and speed of mass spectrometers have been significantly improved these recent years, protein fractionation still represents a limiting step in a proteomic study. Indeed, it has been shown that the improvement of the fractionation step significantly increases the number of proteins that can be subsequently analyzed by a mass spectrometer, thereby resulting in a dramatic increase of the dynamic range of an analysis [START_REF] Sheng | Multidimensional liquid chromatography separation of intact proteins by chromatographic focusing and reversed phase of the human serum proteome: optimization and protein database[END_REF]. To date, there is no single fractionation strategy that has demonstrated the capacity to cover the whole proteome, but several complementary approaches are now available. Protein fractionation has been traditionally performed using 2-DE. However, this franctionation method exhibits several limitations. It often restricts the analysis to the most abundant proteins, has demonstrated a relative lack of reproducibility and has trouble resolving proteins of extreme hydrophobicity, mass or isoelectric point (pI) [START_REF] Issaq | The role of electrophoresis in disease biomarker discovery[END_REF][START_REF] Skalnikova | Protein fingerprints of anti-cancer effects of cyclin-dependent kinase inhibition: identification of candidate biomarkers using 2-d liquid phase separation coupled to mass spectrometry[END_REF]. To circumvent these limitations and to provide an alternative to 2-DE, liquid-based 2D liquid chromatographic (2D-LC) systems have been developed [START_REF] Lubman | Two-dimensional liquid separations-mass mapping of proteins from human cancer cell lysates[END_REF]. The recent availability of an automated system based on 2D-LC, the Proteome-Lab TM PF2D from Beckman Coulter, facilitates the effective implementation of a 2D-LC experiment. The PF2D system separates proteins in the first dimension according to their pI using chromatofocusing, followed by a fractionation according to hydrophobicity, using reversed phase chromatography in the second dimension. The pI-based fractionation, as performed during the first dimension of the PF2D (and in 2-DE), has been shown to be especially useful to study proteins post-translational modifications, such as phosphorylation [START_REF] Zhu | Protein pI shifts due to posttranslational modifications in the separation and characterization of proteins[END_REF]. Several recent proteomic studies comparing the different available fractionation strategies have shown that the different technologies are highly complementary, with a remarkably low number of common proteins identified by all technologies [START_REF] Mcdonald | Expanding the subproteome of the inner mitochondria using protein separation technologies: one-and two-dimensional liquid chromatography and two-dimensional gel electrophoresis[END_REF][START_REF] Li | Comparison of alternative analytical techniques for the characterisation of the human serum proteome in hupo plasma proteome project[END_REF]. However, it seems clear that 2D-LC allows the identification of more proteins than 2D-gels and may be more suited for novel protein discovery [START_REF] Mcdonald | Expanding the subproteome of the inner mitochondria using protein separation technologies: one-and two-dimensional liquid chromatography and two-dimensional gel electrophoresis[END_REF][START_REF] Li | Comparison of alternative analytical techniques for the characterisation of the human serum proteome in hupo plasma proteome project[END_REF]. Automated 2D-LC fractionation, such as that provided by the PF2D system, is relatively new and the reliability and reproducibility of this fractionation procedure remains to be formally assessed and quantified accurately. It was shown that PF2D could generate relatively reproducible UV profiles, based on the visual examination of the chromatograms [START_REF] Park | Comparative proteomic analysis of b. cenocepacia using two-dimensional liquid separations coupled with mass spectrometry[END_REF][START_REF] Soldi | Proteome profile of human urine with two-dimensional liquid phase fractionation[END_REF][START_REF] Yan | A comparison of drug-treated and untreated hct-116 human colon adenocarcinoma cells using a 2-d liquid separation mapping method based upon chromatofocusing pI fractionation[END_REF]. However, a certain level of variability was due to the dynamic nature of the pI fractionation in the first dimension fractionation. As a result, the pH limits of the fractionation gradient could not be fixed, leading sometimes to small shifts in fractionation. In this work, we present a systematic assessment of the reproducibility of PF2D fractionation, using an improved methodology for the first dimension fractionation. We used a software patch, ensuring that all fractionation experiments would start at the exact same pH value. In addition, we used a prototype pH electrode allowing the accurate pH measurement under pressure. We tested the reliability of this improved system by performing several rounds of PF2D fractionation using aliquots of the same protein extract. We used a semi-automated software tool for peak-to-peak comparison between 2D-LC fractionation experiments, facilitating the systematic analysis of more than 3000 peaks. We provide herein the first formal demonstration of the high level of reproducibility of PF2D fractionation, which underscores the interest of this fractionation method for studies of differential proteomics.

Experimental

Cell cultures -Vero E6 cells (ATCC CRL-1586) were grown to confluence in DMEM (Gibco-Invitrogen), supplemented with 10% heat-inactivated fetal calf serum (PAA) and 1X Bufferall (Sigma-Aldrich). Cells (140 x 10 Chromatofocusing was performed on an HPCF 1-D column (250 x 2.1 mm, Beckman Coulter). This first dimension HPLC module was equipped with a 5 ml sample loop.

The signal was recorded at 280 nm. The pH gradient was generated using Start Buffer (pH 8.5) and Elution Buffer (pH 4), both included in the ProteomeLab TM PF2D kit. The chromatofocusing column was first equilibrated for 130 minutes with Start Buffer at pH 8.5 at a flow rate of 0.2 ml/min, before being loaded with 1.5 mg of the desalted protein extract. The flow-through was collected and after a stable baseline was established (35 minutes), a linear pH gradient was initiated by infusing the elution buffer for 95 minutes with a constant flow rate of 0.2 ml/min. The proteins with a pI < 4 were finally eluted by washing the column with 1 M NaCl. A software patch (SP1 Beckman Coulter) was used to ensure that each fraction collection experiment would start at a pH value of 8.3. Fraction collection from the first dimension was controlled with an in-line pH meter equipped with a prototypic pressure-resistant pH electrode (replacement kit p/n A48657). Fractions were collected every 10 min, except during the pH gradient portion of the run, from pH 8.5 to 4.0, when fractions were collected at 0.3 pH unit-intervals.

The second dimension analysis used a nonporous RPHPLC using a C18 column (4.6

x 33 mm, Beckman Coulter) packed with 1.5 µm non-porous silica and kept at 50 °C in a heated column jacket. Eighteen fractions from the first dimension were injected in the NPS-C18 column and eluted using a water/acetonitrile gradient at 0.75 ml/min.

The injection volume was between 50 and 500 µl depending on the protein concentration in the first dimension fraction. Solvent A was 0.1% TFA in water and solvent B was 0.08% TFA in ACN. The gradient consisted in 100% solvent A for 2 minutes, 0-100% solvent B for 30 min and 100% solvent B for 4 minutes. The UV signal was recorded a 214 nm.

Chromatogram analysis -Chromatograms were integrated using the 32Karat software (Beckman Coulter, Fullerton, CA, USA). Data containing the surfaces and retention times of the peaks were exported as a text file and a peak-to-peak analysis was performed using the 4.020 version of the GC-LC-concordance software (Spectrochrom, Bouc Bel Air, France. http://www.spectrochrom.com/). Briefly, this software converts the chromatograms into histograms and performs pair-wise comparisons. The software automatically finds the best quadratic equation to model 1 0 the peak positions using the retention time tolerance and the equation order specified by the operator. A surface ratio between concordant peaks can also be calculated and used as a constraint to exclude peaks. The retention time tolerance was set to 5% and the surface of the peaks analyzed was set between 0.1% and 100% of the largest peak area. The percentage of concordance between the chromatogram A and B was calculated as follows:

% concordance = (number of peaks in chromatogram A concordant with chromatogram B x 100) / (total number of peaks in chromatogram B)

In addition to the peak-to-peak analysis using the GC-LC-concordance software, each concordant peak was visually examined and validated. When required, some peaks were manually integrated, in particular when a shoulder was not correctly detected by the automated analysis.

Results and discussion

Sample preparation: the procedure used for sample preparation is an important parameter that can drastically affect reproducibility and is particularly important during the design of a differential proteomic study. For our repetitive 2D-LC maps, we prepared aliquots of the same cellular extract that were stored in the lysis buffer at -80 °C, to limit the risks of protein degradation. The simultaneous preparation of the cellular extracts for all repetitive PF2D rounds allowed us to ensure that any residual proteolysis occurring during the lysis procedure would be similar for all experiments, thereby limiting the addition of another variable in the analysis. However, the last steps of desalting/gel filtration of the sample were performed prior to each fractionation experiment, to match the conditions of a differential proteomic study, in which these steps are mandatory and could be responsible for a decrease in reproducibility.

Chromatofocusing fractionation: One of the limitations of the chromatofocusing step in the PF2D system is the difficulty to precisely control the pH gradient from one run to another. Although the fraction collector is triggered by an online pH meter, there is still a significant level of variability, due to the difficulty to measure accurately a pH value under pressure. Indeed, the pH electrode has been designed to detect accurately pH changes, but not necessarily to record precise pH values. As a result, the repeatability of this fractionation step is completely dependent on the repeatability of the pH gradients. In addition, the pH limits of the fractions are not exactly similar between maps, which further complicates the analysis. In this work, we have implemented two significant improvements to the PF2D chromatofusing step. First, we used a novel prototypic electrode (p/n A48657) that has been designed not only to accurately measure changes in pH, but also to give accurate pH measurements under pressure. With this new electrode, the pH value under pressure was 8.51 ± 0.01 at the beginning of the gradient for all three experiments (Figure 1). This pH value measured at a calibrated bench pH meter was 8.55, indicating that the accuracy of the measurements under pressure is in the range of 0.04 pH units.

Second, we used a software patch, which enables the fractions to start precisely at pH 8.30 for each experiment. Consequently, the pH limits for each fraction are exactly the same between experiments (Figure 1), provided that the collection volume is not limiting. As a matter of fact, we increased the collecting time to 10 minutes, because we observed that the volume could be limiting for some fractions of the gradient with the usually recommended collecting time of 8.5 min. Indeed, reproducibility was improved when comparing a collecting time of 10 minutes and 8.5 minutes (data not shown).

Chromatogram concordance analysis: A total of 54 chromatograms were generated, representing more than 3000 peaks. In order to assess accurately the reproducibility of the chromatograms with such a large amount of data, we performed a peak-topeak analysis using a modified version of the GC-LC-concordance software. This software enables the comparison of two histograms obtained after peak integration of the chromatograms (Figure 2). The software uses a mathematical model and the retention time error to determine the theoretical position of a peak (see methods).

Second dimension repeatability: the repeatability of the second dimension is an important parameter to test because it gives the theoretical limit of the repeatability of the 2D-LC experiment. To test for the second dimension repeatability, the same series of first dimension fractions were used several times for a second dimension fractionation. The use of a C18 NPS column in association to a UV detector at 214 nm enables the detection of proteins in the nanogram up to the microgram range [START_REF] Dauly | Protein separation and characterization by np-rp-hplc followed by intact maldi-tof mass spectrometry and peptide mass mapping analyses[END_REF]. For this reason the chromatograms were analyzed with a dynamic range of 10 3 .

When testing the reproducibility of a simple chromatogram containing about 30 peaks, the peak concordance between chromatograms was 100% (Figure 3A).

The most complex chromatogram of a 2D-LC experiment is generally obtained with the fraction of basic proteins (pI > 8.5) not retained by the column. We injected this fraction in triplicate onto the NPS column and we obtained chromatograms of about 80 peaks (Figure 3B). With these complex chromatograms, the peak concordance between chromatograms was 96.76 ± 0.98 % (n = 3) . The lower concordance rate for complex chromatograms may be due to different factors. First, it is important to note that most of the peaks are not well resolved in these complex chromatograms. Indeed, MS analysis of 2D-LC peaks has shown that most of the peaks are in fact composed of a mixture of proteins. A large peak can contain 10, 20 or even more proteins. Therefore, the peak integration of such chromatograms is a critical step, since a small change in a shoulder slope can be responsible for a different integration of the same peak between two chromatograms.

As a result, the same peak will be split into two peaks in one case and not in the other, creating a mismatch immediately detected during the software processing of the chromatograms. Second, this decreased concordance rate may be the result of local deformations of the chromatograms, due to slight variations in the ACN gradient between experiments. An example of these local deformations is shown on Fig 3B, where we observed that the first peaks of the top chromatogram had a slightly higher retention time than the 2 other chromatograms. While not globally affecting the rest of the chromatogram, this nevertheless impacted the resolution of peaks located at around 12 min (Fig 3B).

In any event, the concordance rates remained remarkably high, and our results indicate that the second dimension is repeatable enough to detect small changes in peak intensities.

2D proteins maps: Three 2D maps were obtained after subjecting aliquots of the same protein extract to the whole PF2D fractionation procedure, including sample preparation. The maps were obtained with a starting protein amount of 1.5 mg. In some fractions, the protein amount was barely sufficient to obtain a good UV signal, stressing the necessity of using at least 1.5 mg of protein to obtain a satisfactory map when working with a whole proteome. A visual inspection of the maps revealed good similarities and the differential maps displayed peaks of small intensities, indicating a good reproducibility of the experiments (data not shown).

A peak-to-peak analysis of all chromatograms was performed, leading to an accurate measurement of the reproducibility. Results of the map-to-map reproducibility are displayed in Table I.

When performing the overall analysis of the reproducibility between all three experiments, we concluded that the mean peak concordance was 90.19 ± 4.26 % (n= 54). It should be stressed that these values reflect faithfully the reproducibility of the whole fractionation procedure, including not only the two HPLC steps of 2D-LC but also the last steps of sample preparation, in particular the gel filtration step, and the potential impact of sample storage.

In order to explain the slight variations in reproducibility observed between the three performed experiments, we calculated the residual sum of squares (RSS) for each couple of curves of the pH gradients. The RSS is an evaluation of the distance between two pH curves and we found that the RSS(#1,#2) = 6.7, the RSS(#2,#3) = 8.19 and the RSS(#1,#3) = 13.62. Therefore, it appears that the most reproducible maps have also the most reproducible pH curves, further stressing the importance of this initial fractionation step and the interest of the improvements tested herein. In any event, it remains clear that the reproducibility of a fractionation experiment using PF2D depends on several, non exclusive parameters: the reproducibility of the sample preparation and the effect of sample storage, the reproducibility of the pH curves during the fractionation and the reproducibility of the second dimension gradient. Our results indicate that the reliability of fractionation using PF2D is well suited to perform an accurate differential proteomic study. It should nevertheless be kept in mind that this two-dimensional fractionation method will not permit the complete fractionation of all the proteins from a global cell extract and that the UV intensity of a peak reflects the signal of a mixture of proteins. However, the high level of reproducibility of the 2D maps, as demonstrated in this study, enables a quick identification of the fractions in which an abundant protein is responsible for a UV change. In addition, when a subproteome is analyzed by 2D-LC, the resulting peaks are sharper and less coelution is observed [START_REF] Cussac | Proteomic analysis of anaplastic lymphoma cell lines: identification of potential tumour markers[END_REF]. In most cases, a subsequent quantification strategy by MS will nevertheless be required to evaluate the amount of each protein individually, either at the level of the intact protein [START_REF] Zhao | Comparative proteomic analysis of barrett's metaplasia and esophageal adenocarcinomas using 2-d liquid mass mapping[END_REF] or at the peptide level after tryptic digestion [START_REF] Old | Comparison of label-free methods for quantifying human proteins by shotgun proteomics[END_REF].

Conclusions

Here, we provide the first comprehensive assessment of protein fractionation using PF2D, by performing a systematic peak-to-peak measurement of 2D-LC reproducibility. We took advantage of recent hardware and software improvements in the chromatofocusing step of 2D-LC, in particular the use of a new pH electrode giving accurate measurements under pressure combined with a software modification enabling a fractionation based on fixed pH values. Our results demonstrate the high level of reproducibility of the PF2D system, indicating its suitability for differential proteomic studies. Top: protein UV (214nm) profiles obtained after two injections of the same first dimension fraction (pH 6.8-7.1) on a C18 NPS column. Inset: detail of the peak integration performed by the 32karat software for the analysis. The inset table shows the results of peak concordance.

Bottom: protein UV (214nm) profiles obtained after three injections of the same first dimension fraction (pH>8.5) on a C18 NPS column. The inset table shows the results of peak concordance. Results are expressed as mean ± SD (n = 18).

Tables

6 )

 6 were harvested by scraping, washed once in Phosphate Buffered Saline (Invitrogen) and recovered by centrifugation for 5 min. at 430 x g at room temperature. Cell pellets were then resuspended in lysis buffer (ProteomeLab TM PF2D, Beckman) containing a cocktail of protease inhibitors (Sigma-Aldrich) and lysis was performed according to the ProteomeLab TM PF2D human cell lysis protocol. Protein concentration was determined using the 2D quant protein assay (Amersham, Biosciences, Piscataway, NJ, USA). Aliquots of cell extracts were stored at -80°C until ready for use. Liquid chromatography -Prior to chromatofocusing, an aliquot of the cell extract was thawed, desalted on a PD-10 Sephadex TM G-25 gel filtration column with a 5 kDa cutoff (Amersham Biosciences, Piscataway, NJ, USA) and eluted using the chromatofocusing Start Buffer. 1.5 mg of protein extract was used in each 2D-LC experiment. 2D-LC was performed using the ProteomeLab TM PF2D Protein Fractionation System (Beckman Coulter, Fullerton, CA, USA), which consists of two HPLCs, two UV detectors, an auto sampler and a fraction collector. The first dimension fractionation of PF2D consists in chromatofocusing, based on charge. After collection of the fractions from the first dimension in the collector module, each of them is automatically introduced into the second dimension reversed phase chromatography column, which separates proteins based on their hydrophobicity. The fractions can be finally collected into 96-deep-well plates.

Figure 1 .

 1 Figure 1. Reproducibility of the pH gradients using the optimized PF2D system.

Figure 2 .

 2 Figure 2. Comparative chromatogram analysis using the GC-LC-concordance

Figure 3 .

 3 Figure 3. Representative examples of the repeatability of 2D-LC fractionation.

  Figures

Figure 3 .

 3 Figure 2.

TABLE 1 :

 1 Mean peak concordance between 2D maps (%)

		Peak concordance	Map #2	Map #3
	Map #1	w/o	91.59 ± 3.99	89.24 ± 3.16
		# of concording peaks	1062	1067
	Map #3	w/o	89.76 ± 5.24	NA 1
		# of concording peaks	1074	NA 1

Acknowledgements

This work was supported by grants from INSERM, from Région Midi-Pyrénées, from the Genopôle and Canceropôle Grand-Sud-Ouest and from an ANR-06-MIME-005-01 grant to DGD. ES is a recipient of a "poste d'accueil" from INSERM.

DGD is supported by the CNRS.

The authors would like to acknowledge Beckman Coulter and especially Jerald S.

Feitelson, Michael Simonian and Jeff Chapman for their support and for providing the prototype electrode and software patches. Raymond Loyer, from Spectrochrom, is acknowledged for providing the GC-LC-concordance software and for his useful collaboration in the tuning of the mathematical parameters. 1 9 300