Farid Osman Salem

Farid Naït-Abdesselam

Ahmed Mehaoua

Medina Hadjem

Naït Abdesselam

Early Detection of Myocardial Infarction Using WBAN

Keywords: Wireless Body Area Networks, Mycardial Infarction, ECG, ST elevation, Anomaly detection, CUSUM

published or not. The documents may come

I. INTRODUCTION

Cardiovascular diseases are the leading causes of death worldwide, and Myocardial Infarction (MI), commonly known as heart attack, is the most serious one. Only a quarter of patients are able to recognize the symptoms of MI immediately. According to specialists, emergency aid is most effective if it is given over the next 4 hours after MI. In fact, it is during this period that we can prevent serious complications and damage of heart by early detection of MI in order to alert doctors in real time, especially for patients with high risk.

Recent technological advances in wireless networks, sensors integration and miniaturization allow fundamentally modernizing and changing the way healthcare services are deployed and delivered. Specific Wireless Sensor Networks (WSNs) for medical applications, known as Wireless Body Area Networks (WBAN), or Personal Area Networks (PAN) have been developed and deployed.

These networks consist of several sensors operating on or inside the human body, and transmit collected data to a gateway device that acts as a base station (e.g. smart phone, tablet, PDA, etc.) for real time processing, in order to send medical information and alarms to healthcare providers if an anomaly is detected. An example of a WBAN architecture is shown in figure 1 Using WBAN in order to detect MI is an important issue, especially with expected benefits in term of life saving and cost reduction. According to the third global definition of Myocardial Infarction [START_REF] Thygesen | Task Force for the Universal Definition of Myocardial Infarction[END_REF], ECG is an integral part of the diagnostic process of MI and should be acquired and interpreted correctly within 10 minutes after first symptoms.

The Electrocardiogram (ECG) is a waveform that represents the propagation of electric potentials through the heart muscle with respect to time. The propagation of these potentials results in the quasi-periodic contraction of the heart muscle. Each part of the cardiogram refers to a depolarization or a re-polarization of some region in the heart. The cardiogram consists of five major waves, also known as deflections in the cardiology literature, the P , Q, R, S, and T waves. Figure 2(a) illustrates the basic features and intervals of a one-cycle ECG waveform.

The ECG provides a non-invasive method for investigating heart function. Standard ECG measurements utilize 12 leads or views of the electrical activity of the heart (as shown in figure 3). However, ECG measurements using wireless sensors are generally for ambulatory applications and will typically utilize a subset of these leads. Several ECG sensors were proposed and designed by researchers like in [START_REF] Cao | Design and evaluation of a novel wireless three-pad ecg system for generating conventional 12lead signals[END_REF]- [START_REF] Zhang | A batteryless 19 muw mics/ism-band energy harvesting body sensor node soc for exg applications[END_REF]. There are also a plenty of devices in the market that propose portable ECG monitoring, called Holter ECG which can capture ECG during a long period (24h or 48h) like CardioDay proposed by getemed [START_REF]Getemed medizin -cardioday[END_REF]. Myocardial Infarction is an acute ischemic heart disease characterized by a necrosis (death) of a portion of the heart muscle because of deprivation from oxygen. MI causes a serious disturbance of the cardiovascular system that leads to a direct threat for life.

MI is typically characterized by an elevation in the ST segment of ECG (as shown in figure 2(b)) which is normally iso-electric for healthy subjects (described in figure 2(a)). ST segment elevation is generally one of the first symptom of MI and is usually accompanied by chest pain. But in order to be more specific to MI (or suspicious of MI), the ST elevation must be significant in amplitude (up to 0.2 mV) and prolonged in time (several minutes) as indicated in [START_REF] Thygesen | Task Force for the Universal Definition of Myocardial Infarction[END_REF].

The objective of this paper is to propose an efficient method to detect the deviation in "ST segment" time series based on prolonged ECG recording or real time capture of ECG using WBAN technology.

The rest of this paper is organized as follows. Section II surveys related work. Section III briefly reviews related techniques and presents our approach for early detection of MI. Section IV presents our experimental results. Finally, Section V concludes the paper.

II. RELATED WORK

Various ECG monitoring systems were proposed in the literature. MEDiSN [START_REF] Andjong Hyun | Medisn: Medical emergency detection in sensor networks[END_REF] is a WSN for monitoring patients physiological data including ECG, in hospitals and during disaster events. MEDiSN comprises Physiological Monitors (PMs) and Relay Points (RPs) that self-organize into a multi-hop wireless backbone for carrying physiological data. RECAD [START_REF] Zhou | A real-time continuous cardiac arrhythmias detection system: Recad[END_REF] is a realtime continuous arrhythmias detection system based on WSNs technology. It uses Ambulatory Wireless ECG Sensor (AWES) to provide all time cardiac monitoring services.

Many ECG monitoring devices are available on the market, like CardioNet [START_REF]Cardionet[END_REF], a Mobile Cardiac Outpatient Teleme-try™ (MCOT™) which offers a real time monitoring of ECG and wireless transmission of collected data to the monitoring center. Another example, Medtonic Reveal TX [START_REF]Medtonic[END_REF] which is a subcutaneous arrhythmia detection device.

In order to detect ECG anomalies, classification algorithms are generally proposed and used. The purpose of these algorithms is to classify measurements into 2 classes : normal or abnormal. In particular for MI classification, several methods were explored. In [START_REF] Chang | Myocardial infarction classification with multi-lead ecg using hidden markov models and gaussian mixture models[END_REF], Pei-Chann Chang et al. proposed a classification method of MI based on Hidden Markov Model (HMM) and Gaussian Mixture Models (GMM). Their accuracy was 82.5% with specificity of 79.82% and sensitivity 85.71%. Muhammad Arif et al. [START_REF] Arif | Detection and Localization of Myocardial Infarction using K-nearest Neighbor Classifier[END_REF] used K-Nearest Neighbor (KNN) classifier and got an accuracy of 98.3%, sensitivity of 97% and specificity of 99.6%. In [START_REF] Dhawan | Detection of Acute Myocardial Infarction From Serial ECG Using Multilayer Support Vector Machine[END_REF], Akshay Dhawan et al. used Multilayer Support Vector Machine (SVM) and Genetic Algorithm (GA) to detect MI. They obtained a sensitivity of 86.82% and a specificity of 91.05%. E. S. Jayachandran et al. [START_REF] Jayachandran | Analysis of Myocardial Infarction Using Discrete Wavelet Transform[END_REF] obtained an accuracy of 95% using Discrete Wavelet Transform (DWT).

Yang et al. in [START_REF] Yang | Recurrence quantification analysis and principal components in the detection of myocardial infarction from vectorcardiogram signals[END_REF] transformed linearly a 12-lead ECG signals into 3-lead vector-cardiogram using Dower transformation, and then radial basis neural networks are used for classification. The detection accuracy of MI cases was 97% and accuracy of normal cases was 75%. In [START_REF] Haraldsson | Detecting acute myocardial infarction in the 12-lead ecg using hermite expansions and neural networks[END_REF], Henrik Haraldssona et al. decomposed 12-lead ECG using hermite Basis functions and the resulting coefficients were used as inputs to a Bayesian ANN Classifier that were trained to detect MI. The accuracy obtained for MI cases was 94% and 93.3% for patients without MI. Zheng et al. [START_REF] Zheng | Supervised Classification Models to Detect the Presence of Old Myocardial Infarction in Body Surface Potential Maps[END_REF] used SVM, Naïve Bayes (NB) and Random Forecast (RF) methods. Accuracies obtained for these classifiers were 81.9% for NB, 82.8% for SVM and 84.5% for RF.

These algorithms are generally quite complex and are not adequate for an implementation on a WBAN whose main restrictions are limitation in terms of reduced processing capacity, limited storage and power. The main contribution of this paper is the proposition of:

• A Real-time method for early detection of MI • Autonomy of patients and remote capture of ECG using WBAN • A low power consumption algorithm adapted to WBAN using CUSUM Method III. PROPOSED APPROACH Our proposed system is composed of 3 modules :

• ECG sensors that can capture ECG and transmit it wirelessly like those proposed in [START_REF] Cao | Design and evaluation of a novel wireless three-pad ecg system for generating conventional 12lead signals[END_REF]- [START_REF] Zhang | A batteryless 19 muw mics/ism-band energy harvesting body sensor node soc for exg applications[END_REF] • A Cell phone playing the role of a gateway that will receive the digitized ECG wirelessly, pre-processes it, extracts ECG features (P , Q, R, S, T ) and analyzes the extracted features to detect reliable ST elevation in order to send an alert to healthcare professionals along with the ECG captured and the position of the patient We consider a real deployment scenario where sensors send ECG measurements to a smart phone in real time. The main steps performed by our system consist of : pre-processing the received signal, extraction ECG features and detection of significant ST elevation that indicates a potential MI. ECG Signal received from sensors is digitized but some steps are necessary in order to extract ECG parameters. The diagram in figure 4 details these steps. First, Wavelet decomposition is used in order to clean and denoise the original signal. The process of wavelet decomposition consists of taking the samples at a much lower frequency than the original signal, in order to reduce noise and preserve the QRS complex of the ECG, this technique is detailed in [START_REF] Akansu | Emerging applications of wavelets: A review[END_REF]. Then we carry out a baseline wandering removal, which is generally caused by patient movements, dirty or loose sensors and a variety of other reasons. The baseline wandering consists of changes in iso-electric line position and causes an artefactual data when measuring ECG. The removal is achieved by subtracting the first sample, which is generally the reference point, from the rest of signal. Finally, we extract the main ECG features by first detecting the R peaks and their positions in the signal, then we can extract the other peaks (P , Q, S, T ) respectively and their onsets and offsets based on the R peaks. ST segment amplitude can be calculated from S-offset and T -Onset.

The result of the 2 steps described above is a time series of ST segment amplitudes per cycle of ECG. We consider only positive amplitudes of ST in our system, because ST elevation is more specific to MI than a ST depression [START_REF] Thygesen | Task Force for the Universal Definition of Myocardial Infarction[END_REF]. As described in section I, the diagnostic of MI is complex and we can't consider any ST elevation as a sign of MI, this is why we use a change detection algorithm in order to detect new and significant ST elevation changes in ECG. Moreover, we define a window size and a minimum number of deviations that need to be detected within this window before raising an alarm, the aim is to ensure that the detected deviations are sufficiently prolonged to avoid false alarms.

Our proposed algorithm is based on Cumulative Sum (CUSUM) algorithm, which is a sequential analysis statistical tool that is particularly suited for the identification of deviations with the lowest delay.

To detect change point in the sequence of observations ST = (ST 1 , ST 2 , . . . , ST n ), the CUSUM algorithm uses the log-likelihood ratio of the observation ST i :

S n = n i=1 s i
where

s i = log p θ1 (ST i ) p θ0 (ST i ) (1)
Where p θ (ST i ) is the probability density function (PDF) of ST with parameter θ. CUSUM uses hypothesis testing and assumes that data follows a PDF with parameter θ 0 before the change occurrence (hypothesis H 0 ), and with parameter θ 1 after the change (hypothesis H 1 ).

H 0 : θ = θ 0 and H 1 : θ = θ 1 (2)
In fact, the log likelihood ratio s i is more likely to be negative in hypothesis H 0 (normal patient state) and positive after the change H 1 (sick or abnormal). To detect change, the relevant information lies in the difference between the value of the loglikelihood ratio and its minimum value:

g n = S n -m n and m n = min 1≤j≤n (S j ) (3)
A threshold h is used to reject hypothesis H 0 (if g n ≥ h) and to raise an alarm. The value of the threshold must be chosen subject to low false alarm probability. If ST is normally distributed with variance σ 2 in both hypothesis, and with mean µ 0 in H 0 , and µ 1 in H 1 , g n becomes:

g n = g n-1 + µ 1 -µ 0 σ 2 ST n - µ 1 + µ 0 2 + ( 4
)
µ 1 cannot be known before the change, and we approximate its value by µ 1 = α × µ 0 . Equation 4 becomes:

g n = g n-1 + µ 0 (α -1) σ 2 ST n - µ 0 (α + 1) 2 + (5)
Where {y} + = max (0, y). Our CUSUM detection method is described in algorithm 1, where we consider the sequence ST = (ST 1 , ST 2 , . . . , ST n ) as the extracted ST elevation values, and we use a threshold h ≥ 2. In order to reduce false alarms, we introduce two new parameters to CUSUM algorithm:

• A sliding windows w (as shown in figure 5) that contains the last values of ST elevation to calculate the number of deviations

• The parameter k that indicates the minimum number of deviations (A i = 1 in figure 5) that must be detected in the sliding window w before raising an alarm

IV. EXPERIMENTAL RESULTS

In order to evaluate our proposed approach for the early detection of MI, we use the EDB medical database from the Physionet [START_REF] Goldberger | PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals[END_REF]. This database consists of 90 annotated ambulatory ECG recordings from 79 subjects. These subjects have various heart anomalies (vessel disease, hypertension, coronary artery disease, ventricular dyskinesia, and myocardial infarction). Each data trace is two hours in duration and contains two signals , each sampled at 250 samples per second with 12-bit resolution over a nominal 20 millivolt input range. The sample values were rescaled after digitization with reference to calibration signals in the original analog recordings, in order to obtain uniform scale of 200 ADC units per millivolt for all signals.

Figures 6(a) shows the original ECG signal of a patient with MI for a period of 4 seconds before filtering & denoising to illustrate ST elevation, Although the signal is not denoised yet we can visually see that the ST segment is elevated relatively to the PR segment which indicates the iso-electric level.

Figure 6(b) shows the result of wavelet decomposition and baseline wonder removing. The frequency bands of the original signal were separated in four levels. The second level decomposition of the signal was considered as the ideal ECG signal for the features extraction because it is the least noisy. Figure 6(c) is the final result of filtering and base line correction.

Figure 7 shows the variation of ST elevation amplitudes extracted for a complete ECG record (2 hours) of a patient with MI. Figure 8 shows the raised alarms by the modified version of CUSUM algorithm, when the threshold h ≥ 2, the windows size w = 100 samples and the minimum deviations k = 3. In this case, 41 deviations were detected but only 5 alarms are raised for 3310 ST elevation extracted from the MI ECG. If we change the threshold h = 1, the number of deviations increases to 96 and the number of alarms to 9. Also, if we decrease the minimum changes k to 1 for a threshold h = 2, the number of alarms doubles to 10. Furthermore, the size of the windows w is also an important parameter, if the size of w is too large, we may miss important new ST elevation episodes, and if it is too small, we increase the sensitivity of the algorithm to new changes and then increase false alarms. The choice of the right values of parameters h, k and w is then a tradeoff between false alarms and miss detection.

To evaluate the performance of our proposed ST detection algorithm, we apply our algorithm on a subset of 50 records of the EDB Physionet database which contains MI ECG and other cardiac abnormalities ECG. Each record is 2 hours corresponding approximately to 7200 heart beats per record. We use the Receiver Operating Characteristic (ROC) curve in order to represent the fraction of True Positive Rate (TPR) vs. the False Positive Rate (FPR) with various values of the threshold h. The ROC curve is presented in figure 9.

Detection Rate = T P T P + F N

False alarm rate = F P F P + T N

Where T P , F P , T N , F N are respectively, the number of true positives, false positives, true negatives and false negatives. The results obtained show that the optimal rates of DR and F AR are DR=73% and F AR=5% respectively.

V. CONCLUSION

In this paper, we proposed a new approach for early detection of Myocardial Infarction (MI) in real time manner. The proposed approach is based on the detection of deviations in ST segment elevation in the ECG. We combined Wavelet decomposition and Cumulative Sum Method (CUSUM) for a low power detection system adapted to WBAN. The Wavelet decomposition technique is first used to denoise and filter the original ECG signal and extract the values of ST segments, then the CUSUM algorithm is applied on these extracted values to detect significant deviations of ST elevation in order to raise alarms. To reduce the probability of false alarms, we adapted the CUSUM algorithm by introducing a sliding window with fixed size, to hold the number of deviations that must be detected before raising an alarm. Finally, we applied our proposed approach on a real medical database (the Physionet EDB database) using MI and other cardiac problems ECG.

Our experimental results show that our proposed approach can achieve a detection rate of 73% with a false alarm rate of 5%.

Fig. 1 .

 1 Fig. 1. An example of Wireless Body Area Network

Fig. 2 .

 2 Fig. 2. Normal ECG and ST elevation ECG

Fig. 4 .

 4 Fig. 4. ECG signal processing and features extraction

Fig. 5 .

 5 Fig. 5. Sliding window containing the number of deviations

  (a) MI ECG signal from the edb database (b) MI ECG signal after Wavelet Decomposition (c) ECG after filtering and base line correction

Fig. 6 .

 6 Fig. 6. ECG signal pre-processing

  • A remote healthcare professionals (hospital or medical center) that will receive alerts and patients data in real time and take appropriate decisions about the patient condition

		ECG Wavelet Decomposition
		Baseline wander Removal
			ECG features
		R peaks detection
	P peaks	Q peaks	S peaks	T peaks
	P onsets	Q onsets	S onsets	T onsets
	P offsets	Q offsets	S offsets	T offsets
				ST segment amplitudes

  Calculate the initial value of variance σ 2 ; Calculate the initial value of mean µ; for i = 1 →n do Save the value of actual G in tmp; Calculate the new value of G; if (G ≥ h) then The detected deviation is stored (A i = 1); Calculate the sum of A i in the sliding window; if (sum(A i ≥ k) then Calculate µ using ST i values with A i = 0 ; Remove the first element in the window; Shift whole values by one case to the left; Insert the new ST i value in the last case; end for

	Algorithm 1 ST elevation detection algorithm
	Set the windows size w;
	G=0;
	Raise an alarm ;
	G=0 ;
	else
	G=tmp;
	end if
	else if (G <= 0) then
	G=0;
	end if


Benoit Perthame
email: benoit.perthame@ljll.math.upmc.fr

Cristobal Quiñinao
email: cristobal.quininao@college-de-france.fr

Jonathan Touboul
email: jonathan.touboul@college-de-france.fr

COMPETITION AND BOUNDARY FORMATION IN HETEROGENEOUS MEDIA: APPLICATION TO NEURONAL DIFFERENTIATION

Keywords: Morphogen gradients, Reaction-diffusion systems, Traveling waves, Asymptotic analysis, Boundary formation Mathematics Subject Classification: 35B25, 35B36, 35K57, 82C32, 92C15

We analyze an inhomogeneous system of coupled reaction-diffusion equations representing the dynamics of gene expression during differentiation of nerve cells. The outcome of this developmental phase is the formation of distinct functional areas separated by sharp and smooth boundaries. It proceeds through the competition between the expression of two genes whose expression is driven by monotonic gradients of chemicals, and the products of gene expression undergo local diffusion and drive gene expression in neighboring cells. The problem therefore falls in a more general setting of species in competition within a non-homogeneous medium. We show that in the limit of arbitrarily small diffusion, there exists a unique monotonic stationary solution, which splits the neural tissue into two winner-take-all parts at a precise boundary point: on both sides of the boundary, different neuronal types are present. In order to further characterize the location of this boundary, we use a blow-up of the system and define a traveling wave problem parametrized by the position within the monotonic gradient: the precise boundary location is given by the unique point in space at which the speed of the wave vanishes.

Introduction

In this paper we undertake a rigorous mathematical analysis of the boundary formation in a model of developing tissue. Our motivation can be traced back to the work of Alan Turing in the middle of last century, that lead to his celebrated theory of instabilities [START_REF] Turing | The chemical basis of morphogenesis[END_REF]. In his paper, Turing proposed, before substantial knowledge about the development and maturation of living systems was acquired, that the determination of territories was the result of the competition between different chemical substances, he called morphogens, that were reacting together and diffusing, in the presence of a third specie which acts as a catalyst on the expression of both species. In a certain regime of diffusion, these equations yield what we now call Turing patterns, that define a partition of the tissue into differentiated areas (expressing one or the other chemical specie), whose exact shape and location are unpredictable and depend on the initial condition.

In contrast to this indeterminacy of the boundary location in Turing's model, morphogenesis in living systems is an extremely reliable process. Actually, precision of the boundary location is crucial from an evolutionary perspective, in that it ensures proper transmission of essential hereditary patterns. Notwithstanding this qualitative distinction, several years after introduction of Turing's model, biological experiments validated Turing's intuition: transcription factors (called homeoproteins) expressed in cells during development have been shown to have self-activating and reciprocal inhibitor properties as in Turing's theory, but moreover, where shown to have the property to exit the cellular nucleus and membrane and enter the neighboring cells nucleus where it exerts its transcriptional properties [START_REF] Prochiantz | Can transcription factors function as cell-cell signalling molecules?[END_REF][START_REF] Layalle | Engrailed homeoprotein acts as a signaling molecule in the developing fly[END_REF]. However, in contrast to the initial Turing model, the catalyst chemical specie show a specific spatial organization: it forms one-dimensional monotonic gradients of concentration [START_REF] Brunet | The topological role of homeoproteins in the developing central nervous system[END_REF]. This arrangement of catalysts along gradients lead to the development of the french flag model (FFM) [START_REF] Wolpert | Positional information and the spatial pattern of cellular differentiation[END_REF]. This model assumes no diffusion of genetic material, but only all-or-none response to specific thresholds of the catalyst gradient, therefore yielding boundary at a specific location in space. However, this model is relatively sensitive to noise and necessitates to introduce finer mechanisms in order to ensure robustness and accuracy of the boundary location [START_REF] Meinhardt | Space-dependent cell determination under the control of a morphogen gradient[END_REF][START_REF] Gierer | A theory of biological pattern formation[END_REF].

Combining both phenomena of non cell-autonomous activity (small diffusion of transcription factors, acting as Turing morphogens) and graded expression of a catalyst (FFM-like model) lead to a recently developed minimalistic model of boundary formation [START_REF] Quininao | ] local homeoproteins diffusion: combining turing morphogen and french flag mechanisms ascertains boundaries during development[END_REF] reproducing in a parsimonious way both reliability and accuracy of boundary location. This model is given by nonlinear parabolic equations with spatially-dependent coefficients. Simulations indeed showed that in the absence of diffusion, there is no clear separation in two regions, but even very small diffusions disambiguate the differentiation process and lead to a clear definition of the boundary. The object of the paper is to rigorously understand this stabilization in the regime of small diffusions. The mathematical problem we shall be analyzing is actually much more general than the problem of neurodevelopment that motivates the study. Indeed, systems characterized by the competition of two species that are self-activating and reciprocal inhibitor are ubiquitous in life science, and extend to spatially extended population models, large-scale systems of bacterias and social interactions. The particularity of the model we shall analyze, and which may occur in different situations in the cited domains, is precisely the presence of the nonspatially homogeneous catalyst, producing predictable and reproducible patterns.

Due to the ubiquity of such competing systems in life science, we shall propose here a general model supporting reliable pattern formation, and relevant to biology. To this purpose, we complete this introduction by briefly exposing details on neuronal differentiation, before introducing the model we shall be investigating and summarizing our main mathematical results.

1.1. Biological motivation. Let us make more precise the model we have in mind in our developments. The central question we shall address the emergence of reliable boundaries in the developing nervous system. The neural tube indeed develops into a complex functional and anatomical architecture endowed with complex connectivity patterns [START_REF] Raff | The shape of life: genes, development, and the evolution of animal form[END_REF]. The size and shape of functional areas in the cortex is of primary importance: it conditions acquisition of functions, and disruptions are associated to severe conditions, including neuropsychiatric and cognitive disorders [START_REF] Uhlhaas | The development of neural synchrony and large-scale cortical networks during adolescence: relevance for the pathophysiology of schizophrenia and neurodevelopmental hypothesis[END_REF][START_REF] Garey | When cortical development goes wrong: schizophrenia as a neurodevelopmental disease of microcircuits[END_REF]. In the beginning of this century, biologists analyzed developmental genes transcription factors, and showed that these are endowed with non cell-autonomous activity (they belong the homeoprotein family), thanks to two short peptidic sequences present in their DNA-binding domain [START_REF] Joliot | Transduction peptides: from technology to physiology[END_REF]. These transcription factors have the capability to exit the nucleus of the cells, leave the intracellular medium and penetrate the nucleus of neighboring cells where they exert they transcriptional activity. This direct signaling was experimentally demonstrated in vivo during development in the zebrafish [START_REF] Lesaffre | Direct non-cell autonomous Pax6 activity regulates eye development in the zebrafish[END_REF][START_REF] Wizenmann | Extracellular engrailed participates in the topographic guidance of retinal axons in vivo[END_REF], or involved in plasticity of adult networks [START_REF] Beurdeley | Otx2 binding to perineuronal nets persistently regulates plasticity in the mature visual cortex[END_REF][START_REF] Sugiyama | Experience-dependent transfer of Otx2 homeoprotein into the visual cortex activates postnatal plasticity[END_REF][START_REF] Spatazza | Homeoprotein signaling in development, health, and disease: a shaking of dogmas offers challenges and promises from bench to bed[END_REF][START_REF] Brunet | The topological role of homeoproteins in the developing central nervous system[END_REF]. The spatial extension and rate of this process are very low: transcription factors can diffuse and reach at most three cell ranks [START_REF] Layalle | Engrailed homeoprotein acts as a signaling molecule in the developing fly[END_REF], and since the diffusion is passive, important loss reduce the effective number of transcription factors involved. Notwithstanding, it was shown recently [START_REF] Quininao | ] local homeoproteins diffusion: combining turing morphogen and french flag mechanisms ascertains boundaries during development[END_REF] in an elementary model of neurodevelopment that even very low diffusion had major effects on the outcome of the differentiation process. Indeed, in the absence of diffusion, there is an ambiguity in the differentiation in a specific region of the neural tissue, which yield imprecise boundaries that are not reproducible, and sensitively depend on initial condition and possible heterogeneity or noise, but in the presence of small diffusion, the location of the boundary is highly reliable, and the differentiation yields a smooth boundary.

Understanding this dramatic regularization is precisely the object of the present paper. This problem falls in the frame of the competition of two diffusing species A and B that are reciprocal inhibitor and self-activating, with saturation and spatially heterogeneous production rates H A (x, A, B) and H B (x, A, B) (depending on the cell location x). In the neurodevelopment problem, transcription factors expressed by two genes G A and G B constitute our two competing species, and the space heterogeneity corresponds to the graded concentration of morphogens. For simplicity, we shall restrict here our analysis to a one-dimensional case 4 in which the differentiating tissue is along the interval [0, 1]. A schematic version of the model is plotted is Figure 1. (1.1)

   -εd A ∆A = A H A (x, A, B), 0 < x < 1, -εd B ∆B = B H B (x, A, B),
with Robin type boundary conditions stated below.

Here H A and H B are maps from [0, 1] × R + × R + on R, assumed to be of class C 2 . Based on our description of the phenomena, we assume that, for 0

< x < 1, A > 0, B > 0, (1.2)
         H A (x, 0, 0) > 0, H B (x, 0, 0) > 0, ∂ x H A (x, A, B) < 0, ∂ x H B (x, A, B) > 0, ∂ B H A (x, A, B) < 0, ∂ A H B (x, A, B) < 0,
which can be interpreted as follows: on the one hand, the morphogen gradients do not vanish and vary monotically, on the other hand the system expresses competition between species A and B.

Because, we are interested in the limit ε = 0, the solutions of (1.1) in the absence of diffusion are useful. We assume that there exists two solutions (F A (x) > 0, 0) and (0,

F B (x) > 0) (1.3) H A x, F A (x), 0 = 0, H B x, 0, F B (x) = 0,
and that they are respectively stable for x ∈ (0, x a ) and for x ∈ (x b , 1), with x a > x b , i.e., there exists a bistable zone. It means that the linearized matrix at (F A (x), 0) have negative eigenvalues for x ∈ (0, x a ). The same holds at (0, F B (x)) for x ∈ (x b , 1). Moreover, we assume

(1.4) H B (x, F A (x), 0) > 0 for x > x a > x b , H A (x, 0, F B (x)) > 0 for x < x b < x a .
The first inequality, for instance, can be interpreted as follows: for x > x a , A loses stability because resource concentration of B overcomes inhibition from A. 4 Generalization to higher dimensions in situations where geometry of the space and the spatial variations along gradients are sufficiently simple can be handled in the same manner. In [START_REF] Quininao | ] local homeoproteins diffusion: combining turing morphogen and french flag mechanisms ascertains boundaries during development[END_REF], we propose a two-dimensional extension of this property.

Finally, we assume that there exists a unique additional solution (A * (x) > 0, B * (x) > 0) in the interval (x b , x a ) which is a saddle, i.e.

(1.5)

         H A (x, A * (x), B * (x)) = 0, for x b < x < x a , ∂ A H A (x, A * , B * )∂ B H B (x, A * , B * ) -∂ B H A (x, A * , B * )∂ A H B (x, A * , B * ) < 0, ∂ A H A (x, A * , B * ) < 0, ∂ B H B (x, A * , B * ) < 0
that simply express the negativity of the determinant of the Jacobian matrix at this point:

A * ∂ A H A (x, A * , B * ) A * ∂ B H A (x, A * , B * ) B * ∂ A H B (x, A * , B * ) B * ∂ B H B (x, A * , B * ) < 0.
In order to complete the definition of our system (1.1), we need to specify the boundary conditions considered. We are interested in solutions in which the system decomposes the domain into two separate areas in which A or B dominate. In the limit where ε going to zero, it is therefore natural to consider that the system is subject to Dirichlet boundary conditions, but as the diffusion coefficient increases, loss of transcription factor through the boundary becomes increasingly prominent. These mechanisms correspond to Robin (also called third type) boundary conditions:

(1.6)    A(0) - √ ε ∂ ∂x A(0) = F A (0), A(1) + √ ε ∂ ∂x A(1) = 0, B(0) - √ ε ∂ ∂x B(0) = 0, B(1) + √ ε ∂ ∂x B(1) = F B (1).
At this level of generality, assumptions (1.2)-(1.5) may appear formal. These are actually very natural, and we refer to Section 5.1 for a basic example where they are satisfied. They formulate in a general fashion the elements of our problem: the first assumption expresses the existence of two stable differentiated states at both ends of the differentiating tissue in the absence of diffusion, whose domain of stability may overlap. In other words, in the absence of diffusion, levels of concentration of morphogen are sufficient to support differentiated states at the boundaries of the interval, and there exists generically an overlap between these two regions. Within this overlap (in the bistable regime), a saddle fixed point naturally emerges between the two solutions due to the properties of planar vector fields, and in our system, at this fixed point, concentrations of A and B perfectly balance the concentrations of morphogen.

The main result that we will be demonstrating in the present manuscript is the fact that in the presence of small diffusion, a clear boundary between two differentiated domains exists and is unique, and may be characterized univocally. In detail, we shall demonstrate the following:

Theorem 1.1. Under assumptions 1.2-1.5, there exists a classical stationary solution (A ε , B ε ) of (1.1) which satisfies

(1.7) d dx A ε (x) < 0, d dx B ε (x) > 0,
and is obtained as t → ∞ in the corresponding parabolic equation. Moreover

(i) As ε → 0, (A ε , B ε ) converges a.
e. towards a pair (A 0 , B 0 ). These maps are discontinuous at some point x * ∈ [x b , x a ] and have disjoint supports This qualitative result falls in the class of free boundary problems, a well developed asymptotic theory in the frame of homogeneous elliptic or parabolic semilinear equations and systems [START_REF] Carr | Invariant manifolds for metastable patterns in ut = ǫ 2 uxxf (u)[END_REF][START_REF] Fleming | PDE-viscosity solution approach to some problems of large deviations[END_REF][START_REF] Evans | A PDE approach to certain large deviation problems for systems of parabolic equations[END_REF][START_REF] Barles | Wavefront propagation for reaction-diffusion systems of PDE[END_REF][START_REF] Barles | Front propagation for reaction-diffusion equations arising in combustion theory[END_REF][START_REF] Souganidis | Front propagation: theory and applications[END_REF]. As reviewed in [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF], these results generally rely on the definition and analysis of viscosity solutions of the resulting Hamilton-Jacobi equation. The second point of the theorem involves a traveling wave with inhomogeneous speed. A vast literature, in particular in the domain of mathematical ecology of competing populations, have been interested in related questions (see e.g. [START_REF] Aizik | Traveling wave solutions of parabolic systems[END_REF][START_REF] Xin | Front propagation in heterogeneous media[END_REF]). Most applications of this theory are related to front propagations and rules to compute their speeds and invasion properties in homogeneous or heterogeneous environments [START_REF] Berestycki | Front propagation in periodic excitable media[END_REF]. Here, we ask a distinct question concerned with the determination of the precise point where a transition between two stable states occurs.

supp(A 0 ) = [0, x * ] and supp(B 0 ) = [x * , 1]. (ii)
Theorem 1.1 will therefore show the existence of monotonic solutions. The monotonicity is a consequence of analogous properties of the equilibria in the absence of diffusion, which can be readily proved under the current assumptions. Similarly, the monotonicity of equilibria A * (x) and B * (x) can be characterized. This is the object of the following:

Lemma 1.2.
Under assumption (1.2), the functions defined in (1.3) and (1.5) satisfy

(1.8) d dx F A (x) < 0 for x ∈ [0, x a ), d dx F B (x) > 0 for ∈ (x b , 1],
(1.9)

d dx A * (x) > 0 and d dx B * (x) < 0, x b < x < x a .
Proof. Since F A (x) is a fixed point of the system in the absence of diffusion, we have:

d dx H A (x, F A (x), 0) = ∂ x H A (x, F A (x), 0) + ∂ A H A (x, F A (x), 0) d dx F A (x) = 0,
and therefore

d dx F A (x) = - ∂ x H A (x, F A (x), 0) ∂ A H A (x, F A (x), 0) . Assumption (1.4) ensures that ∂ A H A (x, F A (x), 0) < 0 readily implies that d dx F A (x) < 0 for x in [0, x a ). By a similar argument, d dx F B (x) > 0 for any x in (x b , 1]. Hypotheses (1.
2) and (1.5) also constrain the monotonicity of A * and B * . Indeed, since the vector function (H A , H B ) is constant along the curve (x, A * (x), B * (x)), we have

d dx A * (x) = ∂ B H A ∂ x H B -∂ x H A ∂ B H B ∂ A H A ∂ B H B -∂ B H A ∂ A H B , d dx B * (x) = ∂ x H A ∂ A H B -∂ A H A ∂ x H B ∂ A H A ∂ B H B -∂ B H A ∂ A H B .
Using the assumptions (1.2) and (1.5), we conclude the inequalities (1.9).

The manuscript is devoted to the demonstration of Theorem 1.1, and to the development of an application to a specific model of neuronal differentiation. We shall start by proving the existence of a monotonic solution of the elliptic system (1.1), (1.6) by analyzing the long-time properties of the associated parabolic system. The proof of the existence of monotonic solutions and the characterization of the boundary combines stability and monotonicity arguments, WKB asymptotics and a suitable dilation of the spatial variable. The proof proceeds as follows: the limit where ε → 0 is investigated in section 3 and we will show existence and uniqueness of the boundary point x * for small diffusions, and in section 4, we characterize the boundary point x * as the value when a certain traveling wave problem has zero speed, completing the proof of Theorem 1.1. Section 5 puts in good use this theory on a simple model of neuronal differentiation.

Analysis of the parabolic problem

We start with the parabolic problem associated with (1.1)

(2.1)    ∂ t A -εd A ∆A = A H A (x, A, B), 0 < x < 1, t ≥ 0, ∂ t B -εd B ∆B = B H B (x, A, B),
completed again with the Robin boundary conditions (1.6).

We show that for a well chosen pair of initial conditions, solutions to the parabolic problem (2.1)-(1.6) are monotonic in time. Since all coefficients are regular, solutions are classical and therefore bounded.

From here, the existence of steady states is granted.

Because F B is an increasing function in (x b , 1] we can expect that any non-negative solution for the second equation of (2.1) is upper bounded by F B (1). Under the change of variables B = F B (1) -B, system (2.1) becomes

(2.2)    ∂ t A -εd A ∆A = AH A x, A, F B (1) -B , 0 < x < 1, t ≥ 0, ∂ t B -εd B ∆B = -(F B (1) -B)H B x, A, F B (1) -B ,
with the respective boundary conditions

(2.3)    A(0) - √ ε ∂ ∂x A(0) = F A (0), A(1) + √ ε ∂ ∂x A(1) = 0, B(0) - √ ε ∂ ∂x B(0) = F B (1), B(1) + √ ε ∂ ∂x B(1) = 0. hypothesis 1.2, 1.
3 and 1.4 imply that the pair (0, 0) (respectively (F A (0), F B (1))) is a sub-solution (resp. super-solution) of the steady state problem related to (2.2)-(2.3). Therefore, taking (0, F B (1)) as initial condition in (2.1) we have the existence of a regular solution (A ε (t, x), B ε (t, x)) such that:

(2.4) 0 ≤ A ε (t, x) ≤ F A (0) and 0 ≤ B ε (t, x) ≤ F B (1), t ≥ 0, 1 ≤ x ≤ 1.
Lemma 2.1. Then for all t ≥ 0 and x ∈ [0, 1], we have

∂ t A ε (t, x) ≥ 0 and ∂ t B ε (t, x) ≤ 0. Proof. Defining u := ∂ t A ε and v := ∂ t B ε , we have ∂ t u -d A ε∆u = uH A + A ∂ A H A u + A ∂ B H A v, ∂ t v -d B ε∆v = vH B + B ∂ A H B u + B ∂ B H B v,
multiplying the first equation by u -:= min{0, u}, the second one by v + := max{0, v} and integrating over [0, 1] we get 1 2

d dt u 2 -+ d A ε |∂ x u -| 2 -d A εu -∂ x u - 1 0 = u 2 -H A + A ∂ A H A + A ∂ B H A u -v, 1 2 d dt v 2 + + d B ε |∂ x v + | 2 -d B εv + ∂ x v + 1 0 = v 2 + H B + B ∂ B H B + B ∂ A H B v + u.
Time continuity of (A ε (t, x), B ε (t, x)) together with initial conditions imply that for any x:

u -(0, x) = 0 and v + (0, x) = 0.

Thus, there exists C > 0 such that

d dt (u 2 -+ v 2 + ) ≤ C (u 2 -+ v 2 + ),
with zero initial condition. We conclude using Grönwall's lemma.

2.1. Monotonicity in space. We have shown that the monotonicity property of the maps H A and H B in space implies monotonicity of F A (x) and F B (x), solutions of the zero diffusion problem at location x. This is also true of the maps (A ε , B ε ) solutions of the parabolic equation (2.1). In detail, we show that monotonic initial conditions ensure monotonic solutions (A ε , B ε ) in space for all times. This property has two remarkable implications: time dependent solutions belong to the bounded variation class and also their respective steady states.

Lemma 2.2. For any ε > 0 fixed, let us consider any solution (A ε , B ε ) of (2.1)-(1.6) with initial conditions A(0, x) decreasing and B(0, x) increasing. Under assumption (2.4), we have for all t ≥ 0

∂ ∂x A ε (t, x) ≤ 0 and ∂ ∂x B ε (t, x) ≥ 0, 0 ≤ x ≤ 1.
Proof. The proof proceeds as that of Lemma 2.1: we define u := ∂ x A ε and v := ∂ x B ε , we have

∂ t u -d A ε∆u = uH A + A ∂ x H A + A ∂ A H A u + A ∂ B H A v, ∂ t v -d B ε∆v = vH B + B ∂ x H B + B ∂ A H B u + B ∂ B H B v,
multiplying the first equation by u + and the second one by v -, integrating over [0, 1] and using that boundary conditions (1.6) and (2.4), we get 1 2

d dt u 2 + + d A ε |∂ x u + | 2 -d A εu + ∂ x u + 1 0 = u 2 + H A + A ∂ A H A + A ∂ B H A u + v + A ∂ x H A u + ≤ C u 2 + -C u + v -, 1 2 d dt v 2 -+ d B ε |∂ x v -| 2 -d B εv -∂ x v - 1 0 = v 2 -H B + B ∂ B H B + B ∂ A H B v -u + B ∂ x H B v -≤ C v 2 --C u + v -,
where we have also used that ∂ x H A ≤ 0 and ∂ x H B ≥ 0. It is then easy to see that there exists C > 0 such that:

d dt (u 2 + + v 2 -) ≤ C (u 2 + + v 2 -),
with zero initial condition, and to conclude the proof using Grönwall's lemma.

We have therefore constructed a pair (A ε (t, x), B ε (t, x)) such that (2.4) is satisfied. We can apply Lemma 2.2 and find that for any time A ε (t, x) is decreasing and B ε (t, x) increasing. Moreover, Lemma 2.1 together with (2.4) imply that pointwise in space A ε (t, x) (resp. B ε (t, x)) converges to A ε (x) (reps. B ε (x)) solution to (1.1), together with the boundary condition (1.6), in the weak sense.

Bootstrap method allows us to conclude that

A ε (x), B ε (x) ∈ C 2 (0, 1) ∩ C 0 ([0, 1]
) which proves the first part of Theorem 1.1 and (1.7).

2.2.

Positivity of the solutions. We now consider the pair (A ε , B ε ) solution of the stationary problem (2.1)-(1.6) for ε > 0 fixed. We now provide finer estimates of sub-solutions in order to control A ε (0) and B ε (1) away from zero.

Proposition 2.3. There exists ε 0 > 0 such that for any ε < ε 0 , A ε is strictly positive and A ε (0) is, uniformly in ε, larger than some δ A > 0. The same holds for B ε and B ε (1).

Proof. The proof consists in finding a strictly positive sub-solution for

(2.5) -εd A d 2 dx 2 φ A = H A x, φ A (x), B ε (x) φ A , i.
e., the equation for A ε when B ε is fixed. To this purpose, we analyze a completely solvable linear problem related to (2.5), whose solution constitutes a sub-solution of (2.5) and is defined and strictly positive up to the boundary. This solution can thus be used to find a lower bound for A ε (0).

Consider the following linear equation

(2.6) -εd A d 2 dx 2 φ A = min 0≤s≤FA(0) H A 1, s, F B (1) φ A
with boundary conditions inherited from (1.6):

(2.7)

φ A (0) - √ ε d dx φ A (0) = F A (0), φ A (1) + √ ε d dx φ A (1) = 0.
Clearly, the solution takes the form

φ A (x) = α ε e x √ µ/ε + β ε e -x √ µ/ε , with d A µ = - min 0≤s≤FA(0) H A 1, s, F B (1) > 0.
Using (2.7), one can find the exact values of α ε and β ε as a function of the system parameters

α ε β ε = ( √ µ -1) ( √ µ + 1) e -2 √ µ/ε and β ε = F A (0) √ µ + 1 1 - ( √ µ -1) 2 ( √ µ + 1) 2 e -2 √ µ/ε -1
.

Taking ε → 0, we immediately compute

α ε e √ µ/ε → 0 and β ε → β := F A (0) √ µ + 1 ,
thus, for any ε > 0 small enough, φ A becomes positive and 0 < min

0≤x≤1 φ A (x) ≤ max 0≤x≤1 φ A (x) ≤ |α ε |e √ µ/ε + β ε ≤ F A (0).
Then, using that H A is decreasing in both x and B, we obtain

-εd A d 2 dx 2 φ A = min 0≤s≤FA(0) H A 1, s, F B (1) φ A ≤ H A x, φ A (x), B ε (x) φ A .
Therefore, φ A is a sub-solution to (2.5) comprised between 0 and F A (0). Since A ε is a solution to the same problem with the same bounds and φ A (0) is converging to β, the existence of δ A > 0 follows.

Asymptotic analysis as ε vanishes and front position

We now consider the monotonic stationary solutions (A ε , B ε ) for ε > 0 defined in Theorem 1.1. Thanks to Proposition 2.3, we know that for any x ∈ [0, 1]

(3.1) 0 < A ε (x) ≤ F A (0) and 0 < B ε (x) ≤ F B (1).
We are now in a position to demonstrate the convergence of the pair (A ε , B ε ) as ε → 0 towards a pair (A 0 , B 0 ) that are discontinuous at the same point x * and are characterized by point (i) of Theorem 1.1. The proof proceeds as follows. First, using the monotonicity of (A ε , B ε ) we find the existence of A 0 and B 0 , and we characterize those limits as a family of critical points of (1.1) indexed by x. That characterization gives us three possibilities for the support of A 0 . Using a WKB change of variables and the monotonicity properties of critical points (characterized by lemma 1.2), we discard two of them. This allows to conclude on the existence of a unique x * ∈ [x b , x a ] with the properties stated in Theorem 1.1.

3.1.

The limit as ε vanishes. We recall that by monotonicity and L ∞ bounds, the total variations of A ε and B ε are uniformly bounded in ε. Classical theory of Bounded Variation functions (see for instance [13, Theorem 4, p.176]) ensures that there exists a subsequence ε k and BV-functions A 0 , B 0 such that, almost everywhere and in all L p (0, 1), 1 ≤ p < ∞,

(3.2)    A ε k -→ A 0 , 0 ≤ A 0 (x) ≤ F A (0), d dx A 0 ≤ 0, B ε k -→ B 0 , 0 ≤ B 0 (x) ≤ F B (1), d dx B 0 ≥ 0, .
Those limits satisfy, almost everywhere,

(3.3)    A 0 H A x, A 0 (x), B 0 (x) = 0, B 0 H B x, A 0 (x), B 0 (x) = 0.
This means that at each point x, (A 0 , B 0 ) is one of the four nonnegative equilibrium points; (0, 0) and those three given by hypothesese (1.3), (1.5). Because A 0 is decreasing, three possible scenarios arise:

(a) There exits x * such that (A 0 (x), B 0 (x)) = (F A (x), 0), for x < x * and (A 0 (x), B 0 (x)) = (0, F B (x)), for x > x * . (b) There exists x * -< x * + such that (A 0 (x), B 0 (x)) = (F A (x), 0), for x < x * -, (A 0 (x), B 0 (x)) = (0, 0), for x * -< x < x * + .

(c) There exists x * ≥ x b such that (A 0 (x), B 0 (x)) = (F A (x), 0), for x < x * and (A 0 (x), B 0 (x)) = (A * , B * ) for x > x * close enough to x * .

Notice that neither (a) nor (b) exclude the possibility that A 0 is identically zero. Indeed, at this stage, x * (or x * -) could be the origin. Our aim now is to show that only the first scenario is possible for some x * ∈ [x b , x a ] proving part (i) of Theorem 1.1. Scenario (c) can be readily discarded. Indeed, if (c) holds, then the relationship (1.9) would be in contradiction with the monotonicity of A 0 (x) in a neighborhood of x * .

WKB change of unknown. We define ϕ ε

A := -√ ε log(A ε ), which is well defined thanks to Proposition 2.3. Furthermore,

d dx ϕ ε A = - √ ε d dx A ε A ε and d 2 dx 2 ϕ ε A = - √ ε d 2 dx 2 A ε A ε - | d dx A ε | 2 A 2 ε ,
and we find that ϕ ε A is solution of the eikonal equation

d dx ϕ ε A 2 - √ ε d 2 dx 2 ϕ ε A = -H A (x, A ε , B ε ), with d dx ϕ ε A (0) = F A (0) A ε (0) -1, d dx ϕ ε A (1) = 1.
The same constructions can be made for ϕ ε B . If we prove that the family (ϕ ε A ) has some regularity, then we can take let ε go to 0 in ϕ ε A and ϕ ε B . That is the object of the following: Lemma 3.1. The sequence (ϕ ε A ) is uniformly Lipschitz-continuous with respect to ε. Therefore, after extracting a subsequence,

ϕ ε k A -→ ε k →0
ϕ 0 A , a Lipschitz continuous, non-decreasing viscosity solution of

(3.4) d dx ϕ 0 A 2 = -H A (x, A 0 , B 0 ).
The same construction for B ε provides us with a function ϕ B , Lipschitz continuous, non-increasing viscosity solution of

(3.5) d dx ϕ 0 B 2 = -H B (x, A 0 , B 0 ).
Proof. Since A ε ≥ 0 and d dx A ε ≤ 0 we get directly that d dx ϕ ε A ≥ 0. We are going to prove that there exists C ε0 , independent of ε, such that

0 ≤ sup x∈[0,1] d dx ϕ ε A (x) ≤ C ε0 . Consider y, one argmax of d dx ϕ ε A (y). If 0 < y < 1, then d 2 dx 2 ϕ ε A (x) = 0 and d dx ϕ ε A 2 = -H A (x, A ε , B ε ),
which is uniformly upper-bounded because H A is continuous and evaluated on (0, 1) × (0, F A (0)) × (0, F B (1)). The upper bound follows.

If y = 0, Proposition 2.3 tells us that A ε (0) is bounded from below by some positive constant δ A independent from ε. Then, we may conclude again because

d dx ϕ ε A (0) = F A (0) A ε (0) -1 ≤ F A (0) δ A -1 < ∞.
If y = 1, we immediately conclude thanks to the boundary condition and thus, we have proved the uniform Lipschitz estimate.

The Ascoli-Arzela theorem allows us to take a subsequence of ϕ ε A which converges uniformly and we conclude thanks to the usual theory of viscosity solutions [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF][START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF]. Note that the viscosity procedure only allows to control the limsup or liminf of the right hand sides of (3.4), (3.5), and this information sufficient for the conclusion we want to draw.

A direct consequence of Lemma 3.1 is that scenario (b) cannot hold. Indeed, in that case,

d dx ϕ 0 A 2 = -H A (x, 0, 0) < 0, ∀x ∈ (x * -, x * + ) which is contradictory.
The only possible scenario is therefore (a). In order to conclude the proof, we are left showing that x * ∈ [x b , x a ]. It suffices to show that A 0 (x) becomes positive when x → 0 and the same with B 0 (x) when x → 1. Proof. Let us assume that there exists y ∈ (x a , 1) such that A 0 (y) > 0. We have shown that we are necessarily in scenario (a), which implies that B 0 (y) = 0 and by monotonicity

B 0 (x) = 0, A 0 (x) = F A (x) for 0 ≤ x ≤ y.
Using the fact that (F A (x), 0) is linearly unstable for x ∈ (x a , 1] and that ϕ 0 B is a viscosity solution of (3.5), we have

d dx ϕ 0 B 2 = -H B (y, F A (y), 0) < 0,
which is impossible, hence A 0 ≡ 0 on (x a , 1). The same argument ensures us that B 0 ≡ 0 in (0, x b ). One can therefore define the intervals I a and I b by maximality as the supports of A 0 and B 0 .

Characterization of the Front

Now that we have proved the existence of a boundary x * , we can turn to the characterization of this point. To this purpose, we start defining the point,

x * ε such that A ε (x * ε ) = B ε (x * ε )
, which, by monotonicity, is unique. We also know, by compactness and unique limit, that x * ε → x * when ε → 0.

We perform the change of variables y = (x -x * ε )/ √ ε, and define a ε (y) = A ε (x * ε + √ εy) and b ε in the same way. System (5.1) becomes

         -d A d 2 dy 2 a ε (y) = a ε (y)H A x * ε + √ εy, a ε (y), b ε (y) , -d B d 2 dy 2 b ε (y) = b ε (y)H B x * ε + √ εy, a ε (y), b ε (y) , a ε (0) = b ε (0).
Because a ε and b ε are uniformly bounded, by elliptic regularity they are uniformly bounded in C 2 and, after extraction of a subsequence (by uniqueness, as we will show, in fact the full sequence converges), we may take the limit as ε → 0 (which we know is well defined, bounded, Lipschitz-continuous). We find that this limit, denoted (a 0 , b 0 ), is solution of (4.1)

         -d A d 2 dy 2 a 0 (y) = a 0 (y)H A x * , a 0 (y), b 0 (y) , ∂ y a 0 (y) ≤ 0, -d B d 2 dy 2 b 0 (y) = b 0 (y)H B x * , a 0 (y), b 0 (y) , ∂ y b 0 (y) ≥ 0, a 0 (0) = b 0 (0).
This solution is characterized as follows:

Theorem 4.1. The limits satisfy a 0 = 0, b 0 = 0 and there exists a unique value x * such that the system (4.1) has a non-trivial solution. This solution is the unique traveling wave defined as

(4.2)                -c(x) ∂ ∂y a(y; x) -d A ∂ 2 ∂y 2 a(y; x) = a(y; x)H A x, a(y; x), b(y; x) , y ∈ R, -c(x) ∂ ∂y b(y; x) -d B ∂ 2 ∂y 2 b(y; x) = b(y; x)H B x, a(y; x), b(y; x) , lim y→-∞ a(y; x) = F A (x), lim y→+∞ a(y; x) = 0, lim y→+∞ b(y; x) = F B (x), lim y→-∞ b(y; x) = 0,
with speed zero, that is c(x * ) = 0, and connecting

(F A (x * ), 0) to (0, F B (x * )).
Proof. The proof is split into three steps. First we show that functions a ε and b ε cannot converge both at the same time to the zero function. Then, using that a 0 and b 0 converge at -∞ to solutions of (3.3), we show that limit conditions of (4.2) are satisfied. Finally, thanks to a monotonicity argument on the speed c(x), we show that (a 0 , b 0 ) are in fact the unique traveling wave solutions of (4.2) such that c(•) = 0.

1st step. The pair (a ε , b ε ) does not converge to the zero function.

Indeed, for any interval (y -, 0) with y -< 0, integrating by parts the equation on a ε after dividing it by a ε , we compute

1 d A 0 y - H A x * ε + √ εy, a ε (y), b ε (y) dy = - 0 y - | d dy a ε | 2 a 2 ε dy - d dy a ε a ε 0 y - ≤ - d dy a ε (0) a ε (0) .
Moreover, Lemma 3.1 tells us that, for ε < ε 0 ,

d dx ϕ ε A = - √ ε 1 A ε d dx A ε < C ε0 .
This implies directly that for any y ∈ R

- d dy a ε (y) a ε (y) = - 1 A ε (x * ε + √ εy) d dy A ε (x * ε + √ εy) = - √ ε dAε dx (x * ε + √ εy) A ε (x * ε + √ εy) ≤ C ε0 .
Taking the limit ε → 0, using the continuity of H A and that (a ε , b ε ) → (a 0 , b 0 ) uniformly, we find

(4.3) 1 d A 0 y - H A x * , a 0 (y), b 0 (y) dy ≤ C ε0 .
If (a 0 , b 0 ) ≡ (0, 0), then the left hand side becomes |y -|H A (x * , 0, 0)/d A which goes to ∞ when y -→ -∞. Therefore, one of them, say a 0 is positive in some interval and by the strong maximum principle, a 0 (y) > 0 for any y ∈ R. By the condition a 0 (0) = b 0 (0), then b 0 is also positive.

2nd step. The pair (a 0 , b 0 ) satisfies the conditions at infinity in (4.2).

We treat for instance the limit at -∞. Again by elliptic regularity and thanks to (4.3), d 2 dy 2 a 0 (y) and

d 2
dy 2 b 0 (y) vanish at -∞. Therefore the limits of a 0 and b 0 are steady state solutions with a 0 (-∞) > b 0 (-∞).

The case when this steady state is (A * (x * ), B * (x * )) is discarded by stability hypothesis (1.5) and saturation hypothesis (1.2). Indeed, we can rewrite the system defining ǫ A = A(x * ) -a 0 and ǫ B = b 0 -B(x * ). These functions are always positive and have non negative derivatives. Moreover, both they and their first derivatives, go to zero when y → -∞. We can write

d 2 dy 2 ǫ A (y) ǫ B (y) ≈ -∂ A H A /d A ∂ B H A /d A ∂ A H B /d B -∂ B H B /d B ǫ A (y) ǫ B (y) ,
where the matrix is evaluated at (x * , A * (x * ), B * (x * )) and we have neglected the terms of the type ǫ 2 A , ǫ 2 B and ǫ A ǫ B (which do not play a role in the analysis of the signs when y → -∞). Integrating between -∞ and any value y ≪ -1 we get

0 ≤ d dy ǫ A (y) ǫ B (y) ≈ -∂ A H A /d A ∂ B H A /d A ∂ A H B /d B -∂ B H B /d B y -∞ ǫ A y -∞ ǫ B , which is only possible when ∂ A H A ∂ B H B -∂ B H A ∂ A H B ≥
  c(x) ∂ ∂y Φ a (y; x) -d A ∂ 2 ∂y 2 Φ a (y; x) = M 11 Φ a + M 21 Φ b , Φ a > 0 c(x) ∂ ∂y Φ b (y; x) -d B ∂ 2 ∂y 2 Φ b (y; x) = M 12 Φ a + M 22 Φ b , Φ b < 0.
We now consider the x-derivative: z a (y) = ∂ ∂x w a (y; x) and z b (y

) = ∂ ∂x w b (y; x) satisfying    -c(x) ∂ ∂y z a (y; x) -d A ∂ 2 ∂y 2 z a (y; x) = M 11 z a + M 12 z b + c ′ (x)w a + a ∂ x H A , -c(x) ∂ ∂y z b (y; x) -d B ∂ 2 ∂y 2 z b (y; x) = M 21 z a + M 22 z b + c ′ (x)w b + b ∂ x H B .
Integrate in y against the test function Φ and add the two lines, it remains

0 = c ′ (x) [w a Φ a + w b Φ b ] <0 dy + <0 [aΦ a ∂ x H A + bΦ b ∂ x H B ] dy, 0 < x < 1,
thus c ′ < 0. The uniqueness of x * follows directly.

This result concludes the proof of theorem 1.1. We now use this result on a simple model of differentiating neuronal tissue.

Application

5.1. Model. As discussed in [START_REF] Quininao | ] local homeoproteins diffusion: combining turing morphogen and french flag mechanisms ascertains boundaries during development[END_REF], a classical illustration of neurodevelopment is provided by the compartmentalization of the neural tube in response to the diffusion of the ventral and dorsal morphogens Sonic Hedgehog (SHH) and Bone Morphogenetic Protein (BMP), respectively [START_REF] Ribes | Distinct sonic hedgehog signaling dynamics specify floor plate and ventral neuronal progenitors in the vertebrate neural tube[END_REF][START_REF] Ulloa | Morphogens and the control of cell proliferation and patterning in the spinal cord[END_REF]. In this system, a continuous gradient activates ventral and dorsal genes, transcription factors are reciprocal inhibitor and self-activitor and diffuse through boundaries. This is well-known to result in the clear definition of territories that express distinct transcription factors subsets [START_REF] Ashe | The interpretation of morphogen gradients[END_REF][START_REF] Dessaud | Dynamic assignment and maintenance of positional identity in the ventral neural tube by the morphogen sonic hedgehog[END_REF][START_REF] Dessaud | Interpretation of the sonic hedgehog morphogen gradient by a temporal adaptation mechanism[END_REF][START_REF] Kiecker | Compartments and their boundaries in vertebrate brain development[END_REF].

We analyze a simplified version of the model proposed in [START_REF] Quininao | ] local homeoproteins diffusion: combining turing morphogen and french flag mechanisms ascertains boundaries during development[END_REF], which includes:

• Epigenetic phenomena: the more a specie has been expressed, the more it is likely to be expressed. This phenomenon scales the production rate with a coefficient α i (A, B). • The presence of morphogens with a graded concentration along the neural tissue, F i (x), i ∈ {A, B}, • The self-activation of transcription factors • and the saturation effects, limiting the production rate of each species proportionally to the total concentration within a cell.

• Eventually, diffusion of homeoproteins will be considered, through a small diffusion parameter ε ≪ 1.

We will show that these four mechanisms regulating the gene expression (response to gradients, selfactivation, reciprocal inhibition and saturation) precisely correspond to our theoretical assumptions.

Assuming that the number of cells is large, we consider a space-continuous description of the system, and we denote by A(x) and B(x) the concentrations of transcription factors at location x on the neural tissue. The system described above readily translates into the system of parabolic equations:

∂ t A -εd A ∆A = α A A F A (x) + A) -β A A(A + B), 0 < x < 1,
and a similar equation for B. In this equation, we considered epigenetic phenomena to have linear effects: α A (A, B) = α A A. Therefore, the term α A A is the transcriptional intensity, β A is the saturation parameter, and we assume 0 < α A < β A because saturation will overcome necessarily the self-activation. The parameter d A incorporates the relative level of diffusion of the parameter A compared to that of B (at least one of these constants can be incorporated in the ε). We shall assume that the system is subject to Robin type boundary conditions (1.6).

It is not hard to rescale the system so as to write the stationary solutions in the form:

(5.1)

-εd A ∆A = A F A (x) -A -s A B , 0 < x < 1, -εd B ∆B = B F B (x) -B -s B A ,
where, for simplicity of notation, we use the same terms F A (x) and F B (x) to represent the rescaled action of external morphogen gradients. We introduce the parameters s i as positive constants taking into account the relation between α i and β i :

s A = β A β A -α A > 1 and s B = β B β B -α B > 1.
In the limit ε goes to 0, we look for a decreasing solution A connecting the value F A (0) with 0. The morphogen gradients are monotonic and smooth, assumed to be twice differentiable, defined on the closure of the domain, strictly positive and monotonic. Summarizing, there exists δ > 0 such that for any 0 ≤ x ≤ 1

(5.2)

F A (x) > δ, d dx F A (x) < 0, F B (x) > δ, d dx F B (x) > 0.
We have mentioned that diffusion is extremely small. Non-trivial differentiation at these levels of diffusion would require that steady states for ε = 0 are non-trivial as well. This is why we shall assume that:

(

5.3) ∃ (x a , x b ) ∈ I, x b < x a such that F A (x b ) = s A F B (x b ), F B (x a ) = s B F A (x a ).
A first remark is that combining assumptions (5.2) and ( 5.3) we get that (5.4)

F B (x) < F B (x a ) = s B F A (x a ) < s B F A (x), for x ∈ [0, x a ), F A (x) < F A (x b ) = s A F B (x b ) < s A F B (x), for x ∈ (x b , 1].
We have already noticed that both saturation coefficients s A and s B are greater than 1. For the sake of generality, we make the weaker assumption:

(5.5)

s A s B > 1.
Of course, in these notations, the parabolic system reads:

(5.6)

   ∂ t A -εd A ∆A = A F A -A -s A B , 0 < x < 1, t ≥ 0, ∂ t B -εd B ∆B = B F B -B -s B A ,
with the Robin boundary conditions (1.6). If (1.2)-(1.5) are met for these H A and H B , then Theorem 1.1 allow us to say that starting with monotonic initial conditions, then solution to (5.6) defines a unique point x * as a boundary between the two functional areas considered, disambiguating the boundary location.

To start with, note that assumption (1.2) is valid thanks to (5.2) and that they fit the interpretation for neurodevelopment. They are trivially checked in our case since the maps H A (x, A, B) and H B (x, A, B) are linear. We are therefore left to characterizing the equilibria of the system and their stability.

Lemma 5.1. The properties (1.3)-(1.5) are valid for our model (see Fig. 2). In details, under assumptions (5.2), (5.3) and (5.5) and in the absence of diffusion, we have

(i) (F A (x), 0) is a stable fixed point for x ∈ [0, x a ), (ii) (0, F B (x)
) is a stable fixed point for x ∈ (x b , 1], (iii) and there exists an additional solution, which is saddle, in (x a , x b ).

Proof. First two fixed points are trivial solutions, and their stability is obtained by investigating the eigenvalues of the Jacobian matrix at these points (i) At F A (x), 0 , the Jacobian matrix reads

-F A (x) -s A F A (x) 0 F B (x) -s B F A (x) ,
and (5.4) ensures us that this point is stable only on the region [0, x a ). (ii) The pair 0, F B (x) which is analogous to the previous point and stable on (x b , 1]. (iii) Because of hypothesis (5.5), there is an extra fixed point (A * , B * ) given by

A * = s A F B -F A s A s B -1 , B * = s B F A -F B s A s B -1 .
From (5.4) and (5. Remark 5.2. Let us eventually notice the following fact explaining the topology of the phase plane for x ∈ (x b , x a ). The space R + × R + is partitioned into the attraction basin of (F A (x), 0) and that of (0, F B (x)), in addition to lower-dimensional invariant manifolds. The attraction basins of the fixed point are separated by the one-dimensional stable manifold of the saddle fixed point (A * (x), B * (x)), which is an invariance manifold serving as a separatrix between those trajectories converging to F A , 0 and 0, F B .

By a direct application of Theorem 1.1, the system has a unique differentiated solution in the limit of small diffusion. But when considering only cell-autonomous mechanisms, the bistable region x ∈ (x b , x a ) induces an indeterminacy in the differentiation between two domains: cells may choose independently to differentiate into type A or type B, yielding irregular and non-reproducible boundaries depending on the initial condition. This phenomenon is illustrated in Figure 2, right panel: in the absence of diffusion, the region within the interval (x b , x a ) has an unpredictable behavior that depends on space, while in the presence of even a very small diffusion, ambiguity disappears and a unique steady state emerges (see Figure 2). In that sense, a small diffusion suffices to stabilize the transition. From an evolutionary viewpoint, endowing developmental transcription factors with non diffusion properties is a simple mechanism ensuring dramatic stabilization and robustness of the differentiation process. These numerical simulations further open some new perspectives. Indeed, we observe that the convergence towards the monotonic differentiated solutions seem to occur even when we relax the initial condition monotonicity hypothesis of Theorem 1.1. Moreover, with random initial conditions, we numerically observe that for small times, A converges rapidly to F A in [0, x b ) and B to F B in (x a , 1], before the appearance of two abutting traveling fronts that develop toward the center of the coexistence zone, whose speed decreases as the solution converge. Proving that the theorem persists for general initial conditions remains an open problem.

Figure 1 . 1 . 2 .

 112 Figure 1. Our model describes the dynamics of two species in competition (A, pink and B, violet) responding to monotonic resource distributions (bottom line), with reciprocally inhibitory activity and subject to local diffusion.

  The point x * is characterized by the relation c(x * ) = 0 where c(•) represents the speed of propagation of a traveling wave problem parametrized by x (see equation (4.2)).

Lemma 3 . 2 .

 32 There exists two non empty intervals, namely I b and I a = [0, 1] \ I b , such that B 0 ≡ 0 in I b and A 0 ≡ 0 in I a . Moreover, [0, x b ) ⊂ I b and (x a , 1] ⊂ I a .

2 ∂y 2 2 ∂y 2

 2222 0 contradicting the saddle characterization of (A * (x), B * (x)). 3rd step. Finally because the system is competitive, the positive solutions are unique and, in the case at hand, traveling waves with speed 0. We recall why the speed c(•) is monotonic. Considering the derivatives w a (y) = ∂ ∂y a(y) < 0, w b (y) = ∂ ∂y b(y) > 0 they satisfy    -c(x) ∂ ∂y w a (y; x) -d A ∂ w a (y; x) = M 11 w a + M 12 w b , -c(x) ∂ ∂y w b (y; x) -d B ∂ w b (y; x) = M 21 w a + M 22 w b . The signs M 12 := ∂ B H A < 0 and M 21 := ∂ A H B < 0 are compatible with the Krein-Rutman theory, and by consequence the dual problem has a signed solution 

  [START_REF] Berestycki | Front propagation in periodic excitable media[END_REF], we get that (A * , B * ) is admissible (i.e. both coordinates are non-negative) only in the region [x b , x a ]. Monotonicity properties are trivial from the explicit expression, and the stability is governed by the eigenvalues of the Jacobian matrixJac(A * , B * ) = -A * s A A * s B B * B * ,which has negative determinant (as a consequence of assumption (5.5)). Therefore, its eigenvalues are real with opposite signs, i.e. the point (A * , B * ) is a saddle fixed point, completing the proof.

Figure 2 .

 2 Figure 2. Morphogenesis model with exponential morphogen gradients F A (x) and F B (x) (not shown) and s A = s B = 2. (Left) Equilibria of the system in the absence of diffusion together with their stability (thick solid line: stable, thin solid line: repulsive, dashed: saddle). (Right) Numerical simulations of stationary states of the system (5.6)-(1.6) shows (top) the ambiguity of boundary location for ε = 0 and (bottom) the disambiguation for small diffusion ε = 10 -6 .


Julien Schorsch

Vincent Laurain

Marion Gilson

Hugues Garnier

Refined IV-based method for LPV partial differential equation model identification

This paper presents a direct identification method for linear parameter varying models described by partial differential equations in an input-output setting. The continuous space-time model is firstly rewritten as a multiple-input singleoutput model. The continuous filtering operations are reformulated as a discrete convolution product and a refined instrumental variable technique is developed to efficiently estimate the model parameters. The performance of the proposed method is then illustrated via a representative simulation example.

I. INTRODUCTION

The linear parameter varying (LPV) partial differential equation (PDE) identification case is considered in this paper where the parameters are assumed to be function of a measurable deterministic time-space varying signal. This modeling concept allows for a wide representation of physical processes.

Several identification methods for PDE models are available in the literature (see e.g. [START_REF] Kubrusly | Distributed parameter system identification A Survey[END_REF] for survey) and two main approaches can be distinguished. The indirect approach including two steps involves the transformation of the continuous time-space (CTS) model into a discrete timespace (DTS) model and the identification of the DTS model parameters [START_REF] Ali | Identification of Box-Jenkins models for parameter-varying spatially interconnected systems[END_REF]. The direct approach aims to estimate the coefficients of the CTS model directly from sampled signals ( [START_REF] Chochol | An original differentiation tool for identification on continuous structures[END_REF] and [START_REF] Schorsch | Instrumental variable methods for identifying partial differential equation models[END_REF]). However, only a few papers deal with the PDE identification problem into the LPV framework. The use of LPV models has been applied to the identification of linear time invariant (LTI) PDE system by approximating the PDE at each working point with a DT-LPV model [START_REF] Belforte | LPV approximation of distributed parameter systems in environmental modelling[END_REF].

A method has been already proposed to identify a continuous space-time LPV partial differential equation model [START_REF] Schorsch | Identification of LPV partial differential equation models[END_REF]. The approach is based on an iterative procedure develop to estimate efficiently the model parameters into the errorprediction minimization framework. However, as it will be shown in the numerical example developed in this paper, in the situation where the output is corrupted by an additive colored noise, this method delivers biased estimates.

The proposed direct approach uses here an instrumental variable (IV) based method which provides consistent parameter estimates whatever the kind of noise [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF] or [START_REF] Young | Recursive Estimation and Time-series Analysis: An Introduction for the Student and Practitioner[END_REF]. The Simplified Refined Instrumental Variable for Continuoustime systems (SRIVC) method was recently extended to continuous-time LPV ordinary differential equation (ODE) models [START_REF] Laurain | Direct identification of continuous-time linear parameter-varying input/output models[END_REF] and to linear PDE models [START_REF] Schorsch | Instrumental variable methods for identifying partial differential equation models[END_REF]. It is here adapted for the LPV-PDE models.

An advantage of the SRIVC algorithm is to make use of a filter which reduces the noise influence on the estimated parameters. In the ODE model case, this filter follows the same dynamic as the system. However, in the LPV-PDE case, the filter should be dependent on an external variable ρ, if it is chosen to follow the same dynamics as the system. The influence of this dependency is twofold. First, this dependency is injected on the filtered noise which decreases the quality of the estimated parameters. Secondly, the classical algebraic properties (commutativity, associativity,...) are not guaranteed. The proposed solution is based on the reformulation of the single-input single-output (SISO) LPV-PDE model to a multiple-input single-output (MISO) PDE model. From this reformulation a new filter is emphasized. It is not anymore dependent on ρ, keeps algebraic properties and follows a dynamic part of the considered system.

However, the resulting filter is written in terms of partial derivative operators and therefore it cannot be directly defined as a rational transfer function. Usually, the transfer functions associated with PDE models are irrational functions which involve an infinite number of poles and zeros [START_REF] Curtain | Transfer function of distributed parameter systems : A tutorial[END_REF]. Nevertheless, the solution developed here, is to discretize the filter by using classical numerical schemes and to write the filter operations as a convolution product.

The paper is organized as follows: in Section 2, the general class of CTS-LPV-PDE models is introduced. The reformulation of the system is provided Section 3 including a description of the proposed filter. Section 4 deals with the writing of this filter by using the discretization of the differential operator. A description of the proposed estimator is given in Section 5. Finally, an numerical example is used in Section 6 to illustrate the performance of the proposed approach.

II. PROBLEM FORMULATION

A. Deterministic model

Let Ω a close set and consider the continuous model described by the following PDE:

A(∂, ρ)χ(t, x) = B(∂, ρ)u(t, x) on Ω×]0, +∞[ χ(0, x) = φ(x) on Ω (1)
x ∈ Ω is the space variable, note that for clarity of the presentation, x is supposed to be one-dimensional. t ∈ [0, T ] is the time variable. χ(x, t) is the model output and u(x, t) is the input. φ(x) is a continuous space function which represents the model behaviour for t = 0. For (i t , i x ) ∈ N 2 , ∂ it,ix denotes the partial differential operator:

∂ it,ix χ(t, x) = ∂ it t ∂ ix x χ(t, x) = ∂ it+ix χ(t, x) ∂t it ∂x ix , (2)
where ∂ it t = ∂ i t ∂t i t is the i t -th time partial derivative and ∂ ix x = ∂ ix ∂x ix is the i x -th space partial derivative. To complete Equation (1) boundary conditions needs to be specify. A large number of boundary conditions can be chosen, depending on the problem formulation, the number of variables involved and the nature of the equation.

A(∂, ρ) and B(∂, ρ) are bivariate polynomials in ∂ of degrees n t , n x and m t , m x respectively:

             A(∂, ρ) = p nt,nx + n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 a(ρ) it,ix ∂ it,ix (3a) B(∂, ρ) = mt,mx jt,jx=0 b(ρ) jt,jx ∂ jt,jx (3b)
In order to represent a larger behavior range, the model is assumed to be parameter varying. a it,ix (ρ) and b jt,jx (ρ) are analytic functions with dependence on an external deterministic scheduling variable ρ(t, x) = ρ, which is called static dependence. (n t , n x ) are the orders of the system (assumed to be known). In this work, it is assumed that n x ≥ m x , m t ; n t ≥ m x , m t ; and A(∂, ρ) is a monic polynomial in ∂ t variable.

B. Process model

Assuming ρ a priori known, the A(∂, ρ) and B(∂, ρ) polynomials can be parametrized as:

                     A(∂, ρ, θ) = ∂ nt,nx + n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 a it,ix (ρ)∂ it,ix a it,ix (ρ) = a 0 it,ix + nα l=1 a l it,ix f l (ρ) i t = 0, • • • , n t i x = 0, • • • , n x (4)                    B(∂, ρ, θ) = mt,mx jt,jx=0 b jt,jx (ρ)∂ jt,jx b jt,jx (ρ) = b 0 jt,jx + m β l=1 b l jt,jx g l (ρ) j t = 0, • • • , m t j x = 0, • • • , m x . (5)
In this parametrization, {f l } nα l=1 and {g l } m β l=1 are meromorphic functions of ρ, with static dependence. They can be chosen for example as linearly independent functions. The associated model parameters θ are stacked columnwise:

θ = [a 0,0 • • • a nt,nx-1 b 0,0 • • • b mt,mx ] T ∈ R n θ (6)
where

a it,ix = a 0 it,ix a 1 it,ix • • • a nα it,ix ∈ R nα+1 (7) b jt,jx = b 0 jt,jx b 1 jt,jx • • • b m β jt,jx ∈ R m β +1 (8)
and n θ = ((nt + 1)(nx + 1) -1)(nα + 1) + (mt + 1)(mx + 1)(m β + 1).

C. Noisy model

In this paper, stochastic PDE class is considered. These systems are a generalization of PDE taking into account a random noise term. Thus due to measurement inaccuracies, the measure χ(t, x) is not available but a disturbed version of it is obtained.

Moreover, in practical situations, such as those encountered in environmental science, the measured data are not continuous in time and in space. Given the sampled nature of the data, an obvious assumption is that the sampled output is contaminated by a discrete time-space noise.

The data-generating system takes the following Output-Error structure:

           A(∂, ρ, θ)χ(t, x) = B(∂, ρ, θ)u(t, x) on Ω×]0, T ] χ(0, x) = φ(x) on Ω y(t k , x ℓ ) = χ(t k , x ℓ ) + e(t k , x ℓ ) (9)
where the observed noisy output y(t k , x ℓ ) and χ(t k , x ℓ ) are the sampled signals at the time t k and the observation point

x ℓ for k ∈ {1, • • • , N } and ℓ ∈ {1, • • • , L}. e(t k , x ℓ
) is a white noise, i.e. a zero mean, random spatial array, with no temporal or spatial correlation.

D. Objective

Finally, the objective is to estimate the parameter vector (6) of the model ( 9) based on the measured sampled input, output data and the scheduling variables

Z = {u(t k , x ℓ ), y(t k , x ℓ ), ρ(t k , x ℓ )} L,N ℓ=1,k=1 .

III. THE LPV-PDE REFORMULATION

A. Reformulation

The LPV-PDE model identification problem cannot be formulated directly as a prediction error minimisation because of the non-commutativity of the LPV polynomials. This problem is well-known when dealing with LPV models. In [START_REF] Laurain | Direct identification of continuous-time linear parameter-varying input/output models[END_REF] (for ODE models) and [START_REF] Schorsch | Identification of LPV partial differential equation models[END_REF] (for PDE models), the proposed solution is to rewrite the LPV model as a MISO LTI model for which the PEM can be formulated. The idea here is to follow the same concept and to reformulate the LPV model in order to express it as a linear regression form:

y (nt,nx) (t k , x ℓ ) = ϕ(t k , x ℓ )θ + F (∂, θ)e(t k , x ℓ ) (10)
Such a reformulation can be expressed as:

χ(nt,nx) (t k , x ℓ ) + n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 a 0 it,ix χ(it,ix) (t k , x ℓ ) F (∂,θ)χ(t k ,x ℓ ) + n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 nα l=1 a l it,ix f l (ρ)χ (it,ix) (t k , x ℓ ) χ(i t ,ix ) l (t k ,x ℓ ) = mt,mx jt,jx=0 m β l=0 b l jt,jx g l (ρ)u (jt,jx) (t k , x ℓ ) u (j t ,jx ) l (t k ,x ℓ ) (11)
where g 0 (ρ) = 1 and the superscript (j t , j x ) for a signal, like u (jt,jx) , denotes the j th t partial derivative in time and the j th

x partial derivative in space and

F (∂, θ) = ∂ nt,nx + n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 a 0 it,ix ∂ it,ix (12)
Note that in this way, the main advantage is to transpose the time-space variation of the coefficients onto the signals now expressed as

           χ(it,ix) l (t k , x ℓ ) = f l (ρ)χ (it,ix) (t k , x ℓ ) (13a) {i t , i x , l} ∈ {0 • • • n t , 0 • • • n x , 1 • • • n α } u (jt,jx) l (t k , x ℓ ) = g l (ρ)u (jt,jx) (t k , x ℓ ) (13b) {j t , j x , l} ∈ {0 • • • m t , 0 • • • m x , 0 • • • m β }
Therefore, the process part of the LPV model is rewritten as a MISO system with ((n t + 1)(n x + 1) -1)n α + (m t + 1)(m x + 1)(m β + 1) inputs {χ l it,ix } nt,nx,nα it=0,ix=0,l=1 and {u l jt,jx } mt,mx,m β jt=0,jx=0,l=0 . By using ( 11), ( 9) can be rewritten in terms of the output signal y(t k , x ℓ ) as

             F (∂, θ)y(t k , x ℓ ) = - n t ,nx (i t ,ix ) =(n t ,nx ) it,ix=0 nα l=1 a l it,ix χ(it,ix) l (t k , x ℓ ) + mt,mx jt,jx=0 m β l=0 b l jt,jx u (jt,jx) l (t k , x ℓ ) + F (∂, θ)e(t k , x ℓ ).
(14) which can be written under the regression form of [START_REF] Morton | Numerical Solution of Partial Differential Equations[END_REF] where

ϕ(t k , x ℓ ) = -y (0,0) (t k , x ℓ ) • • • -y (nt,nx-1) (t k , x ℓ ) -χ (0,0) 1 (t k , x ℓ ) • • • - χ(nt,nx-1) nα (t k , x ℓ ) u (0,0) 0 (t k , x ℓ ) • • • u (mt,mx) m β (t k , x ℓ ) ] T . ( 15
)
The regressor is made up with the input u(t k , x ℓ ), the noisy output y(t k , x ℓ ) and noise-free output χ(t k , x ℓ ).

B. Linear regression form

In the given context, contrary to the LPV formulation where the LPV dependency appears into the signal definitions, the filter F (∂, θ) commutes so that (10) can be rewritten as

y (nt,nx) f (t k , x ℓ ) = ϕ T f (t k , x ℓ )θ + e(t k , x ℓ ) (16)
where

y f (t k , x ℓ ), χ(it,ix) l,f (t k , x ℓ ) and u (jt,jx) l,f
(t k , x ℓ ) represent the outputs of a filtering operation.

y f (t k , x ℓ ) = 1 F (∂, θ) y(t k , x ℓ ) (17)
and

ϕ f (t k , x ℓ ) = -y (0,0) f (t k , x ℓ ) • • • -y (nt,nx-1) f (t k , x ℓ ) -χ (0,0) 1,f (t k , x ℓ ) • • • - χ(nt,nx-1) nα,f (t k , x ℓ ) u (0,0) 0,f (t k , x ℓ ) • • • u (mt,mx) m β ,f (t k , x ℓ ) ] T (18)
Moreover, the proposed filter does not depend on the deterministic scheduling variable ρ but only on the LTI part of the PDE. However, since the inverse of ∂ is not defined in the partial derivative context, the digital implementation of the filtering operation is discussed in the next Section.

IV. THE DIFFERENTIAL OPERATOR DISCRETIZATION

In contrary to ODE model which can be expressed by a transfer function model between the excitation signal and output signal, it is more complicated to write PDE under this form. However, for specific linear PDE systems, several approaches are possible. As an example, a method developed in [START_REF] Curtain | Transfer function of distributed parameter systems : A tutorial[END_REF] consists in using the Laplace transforms with respect to the time variable t and then to apply it on the system and the boundary conditions. PDE are then written as irrational transfer functions which can be expanded into infinite dimensional rational transfer functions. A different solution is to use fractional models to describe the PDE system (see e.g. [START_REF] Malti | Tutorial on system identification using fractional differentiation models[END_REF]). Another method consists in taking the Fourier transforms with respect to the space variable x. This method is used when the space domain is assumed to be infinite and permits to obtain an analytic solution of the PDE.

The method used here, to tackle the irrational function problem, is to consider a discretization along the dimension axis, leading to a so-called parameter distributed system. This approach consists in reducing the infinite dimensional system described by PDE to a finite dimensional model, using the input and output observations to construct an approximation of the solution. The continuous filter F (∂, θ) is then transformed into a discrete filter F (ð, θ), where ∂ is approximated by a discretized operator ð = (ð t , ð x ):

∂ t = ð t + O(T e ) ∂ x = ð x + O(h) (19)
where T s and h are the sampling periods in time and space and O(.) is the truncation error. This discrete filter can be expressed as a convolution product

[F (ð, θ)χ](t k , x ℓ ) = (f * χ)(t k , x ℓ ). (20)
Using the simplest partial derivatives approximations made by finite difference method in equation (20), the filter can be realized as a finite-dimensional rational transfer function.

Example:

Consider the heat equation

∂ 1,0 χ(t, x) -∂ 0,2 χ(t, x) = u(t, x) (21)
∂ 1,0 χ(t, x) can be computed by the following Euler forward difference approximation:

[ð 1,0 χ](t k , x ℓ ) = q +1 t -1 T s χ(t k , x ℓ ) (22)
where q t is the shift operator. In the same way, ∂ 0,2 χ(t, x) can be approximated by the central difference approximation:

[ð 0,2 χ](t k , x ℓ ) = q +1 x -2 + q -1 x h 2 χ(t k , x ℓ ). ( 23
)
where q x is the shift operator. Equation ( 21) can then be expressed into the following discrete form

(ð 1,0 -ð 0,2 )χ(t k , x ℓ ) = u(t k , x ℓ ), (24)
which can be written under the following transfer function:

       G(q t , q x ) = T s q +1 t q 0 x + Ts h 2 q 0 t q +1 x + (1 -2Ts h 2 )q 0 t q 0 x + Ts h 2 q 0 t q -1 x χ(t k , x ℓ ) = G(q t , q x )u(t k , x ℓ ) (25)
The stability of this transfer function can be proven by the classical methods used for studying the convergence of the numerical scheme [START_REF] Morton | Numerical Solution of Partial Differential Equations[END_REF]. Here the transfer function G(q t , q x ) is stable when

-1 + 2T s h 2 ≤ 0. (26)
V. ESTIMATION Under the assumption of a priori known F (ð, θ) and χ, the optimal estimator in the prediction error minimization sense is the least-squares (LS) estimator:

θLS =   N,L k,ℓ=1 ϕ f (t k , x ℓ )ϕ T f (t k , x ℓ )   -1 •   N,L k,ℓ=1 ϕ f (t k , x ℓ )y (nt,nx) f (t k , x ℓ )   (27)
In practice, the filters F (ð, θ) and the signals χ(it,ix) are unknown and therefore an iterative procedure is proposed in order to cope with this issue, similarly as it is achieved in the approaches developed for the ODE case in [START_REF] Lohnberg | Iterative Least-Squares Parameter Estimation for ARMA Pulse Response and Output Disturbance[END_REF] and [START_REF] Schorsch | Identification of LPV partial differential equation models[END_REF] for the PDE case. However, the LS method delivers biased estimates in the general practical situation where the noise structure is not known. It is then necessary to develop estimators which give unbiased estimates whatever the noise. In this paper, the proposed estimator is based on the IV method which is known to provide efficient estimates.

A. Basic IV estimator

The main idea behind the IV method is to modify the normal equations ( 27) so that they yield consistent estimates for arbitrary additive noise. This involves the introduction of an IV vector ζ f (x ℓ , t k ) whose elements are correlated with the noise-free regression vector φf (x ℓ , t k ), where

φf (t k , x ℓ ) = F (ð, θ)φ(t k , x ℓ ) (28) with φ(t k , x ℓ ) = -χ (0,0) (t k , x ℓ ) • • • -χ(nt,nx-1) (t k , x ℓ ) -χ (0,0) 1 (t k , x ℓ ) • • • - χ(nt,nx-1) nα (t k , x ℓ ) u (0,0) 0 (t k , x ℓ ) • • • u (mt,mx) m β (t k , x ℓ ) ] T , (29)
but uncorrelated with the additive noise. The basic IV parameter estimates are then given by the following IV normal equations:

θIV =   N,L k,ℓ=1 ζ f (t k , x ℓ )ϕ T f (t k , x ℓ )   -1 •   N,L k,ℓ=1 ζ f (t k , x ℓ )y (nt,nx) f (t k , x ℓ )   ( 30
)
provided that the matrix-inverse exists.

B. Consistency properties

By introducing ( 16) into (30) it comes:

θIV = θ +   N,L k,ℓ=1 ζ f (x ℓ , t k )ϕ T f (x ℓ , t k )   -1 •   N,L k,ℓ=1 ζ f (x ℓ , t k )e f (t k , x ℓ )   (31)
It can be deduced from (31) that θIV is a consistent estimate of θ if

E[ζ f (t k , x ℓ )ϕ T f (t k , x ℓ )] is not singular (32a) E[ζ f (t k , x ℓ )e f (t k , x ℓ )] = 0 (32b)
where E(.) denotes the mathematical expectation.

C. Refined Instrumental Variable algorithm

Theoretically, the IV method corrects the statistical inconsistency of the least squares solution but its practical implementation requires some method for generating an instrument that satisfies the conditions given by (32a) and (32b). Several different procedures for the ODE models have been suggested in the system and control literature. The most successful of these is the statistically optimal SRIVC algorithm which involves an iterative procedure where at each iteration, the auxiliary model used to generate the IV, as well as the associated filters are updated based on the parameter estimates at the previous iteration. The proposed LPV-SRIVC-PDE model identification algorithm can be summarized as shown below. The algorithm consists of three first steps for the initialization and the last five iterative steps are used to estimate the filter, the derivatives and finally the model parameters.

• Step 1 Initialization of θ 0 for l = 0 only:

θ0 = â0 0,0 • • • â0 nt,nx-1 0 • • • 0 b0 0,0 • • • b0 mt,mx 0 • • • 0 T
using for instance the LS-PDE method proposed in [START_REF] Schorsch | Instrumental variable methods for identifying partial differential equation models[END_REF]. In this first step the system is considered as a LTSI system.

• Step 2 Compute a first estimated filter on the basis of the estimates obtained in Step 1

Q(ð, θ0 ) = 1 F (ð, θ0 )
where F (ð, θ0 ) is given as [START_REF] Schorsch | Instrumental variable methods for identifying partial differential equation models[END_REF]. Use the filter in order to generate the estimated derivatives

{u (jt,jx) l,f } mt,mx,m β jt=0,jx=0,l=0 = Q c (ð, θ0 ){u (jt,jx) l } mt,mx,m β jt=0,jx=0,l=0 {y (it,ix) f } nt,nx it=0,ix=0 = Q c (ð, θ0 ){y (it,ix) } nt,nx it=0,ix=0
and the regressor ϕ f (t k , x ℓ ) built up from the filtered noisy output only and the input for this initialization step

ϕ f (t k , x ℓ ) = -y f (t k , x ℓ ) • • • -y (nt,nx -1) f (t k , x ℓ ) -y (0,0) 1,f (t k
, x ℓ ) • • • -y (nt,nx-1) nα,f (t k , x ℓ ) u (0,0) 0,f (t k , x ℓ ) • • • u (mt,mx) m β ,f (t k , x ℓ ) ] T • step 3 Compute an estimate θ1 : θ1 =   N,L k,ℓ=1 ϕ f (t k , x ℓ )ϕ T f (t k , x ℓ )   -1 •   N,L k,ℓ=1 ϕ f (t k , x ℓ )y (n t ,nx) f (t k , x ℓ )   Set τ = 1.
[End of initialization part.]

• Step 4 At any iteration τ +1, compute an estimate of the noise-free output χ(t k , x ℓ ) via numerical approximation of A(ð, ρ, θτ ) χ(t, x) = B(ð, ρ, θτ )u(t, x)

based on the estimates θτ at the previous iteration.

• Step 5 Compute the estimated filters

Q(ð, θτ ) = 1 F (ð, θτ )
and use the filters as well as χ(t k , x ℓ ) in order to generate the estimates of the derivatives:

{ χ(it,ix) l,f } = Q(ð, θτ ){ χ(it,ix) l } for i t = 0, • • • , n t ; i x = 0, • • • , n x and l = 1, • • • , n α {u (jt,jx) l,f } = Q(ð, θτ ){u (jt,jx) l } for j t = 0, • • • , m t ; j x = 0, • • • , m x and l = 1, • • • , m β {y (it,ix) f } = Q(ð, θτ ){y (it,ix) } for i t = 0, • • • , n t and i x = 0, • • • , n x .
• Step 6 Build the estimated filtered regressor φf (t k , x ℓ ) as:

φf (t k , x ℓ ) = -y (0,0) f (t k , x ℓ ) • • • -y (nt,nx-1) f (t k , x ℓ ) - χ(0,0) 1,f (t k , x ℓ ) • • • - χ(nt,nx-1) nα,f (t k , x ℓ ) u (0,0) 0,f (t k , x ℓ ) • • • u (mt,mx) m β ,f (t k , x ℓ ) ]
T φf (t k , x ℓ ) is composed of the filtered noisy output, the filtered noise-free estimated output and the input. And build the filtered estimated IV ζf (t k , x ℓ ) as:

ζf (t k , x ℓ ) = - χ(0,0) f (t k , x ℓ ) • • • - χ(nt,nx-1) f (t k , x ℓ ) - χ(0,0) 1,f (t k , x ℓ ) • • • - χ(nt,nx-1) nα,f (t k , x ℓ ) u (0,0) 0,f (t k , x ℓ ) • • • u (mt,mx) m β ,f (t k , x ℓ ) ]
T ζf (t k , x ℓ ) vector is composed of the filtered noisefree estimated output and the input and represents the instruments.

• Step 7 Compute the LPV-SRIVC-PDE estimate θτ+1 at the (τ + 1)th iteration:

θτ+1 =   N,L k,ℓ=1 ζf (t k , x ℓ ) φT f (t k , x ℓ )   -1 •   N,L k,ℓ=1
ζf (t k , x ℓ )y

(n t ,nx) f (t k , x ℓ )   •
Step 8 If θτ+1 has converged or the maximum number of iterations is reached, then stop, else increase τ by 1 and go to Step 4.

VI. NUMERICAL EXAMPLE

The performance of the described algorithm is evaluated in the case of the advection-diffusion equation (ADE). This equation is often used in the water resource quality analysis to describe the transport and dispersion of a solute into a river channel [START_REF] Tych | A matlab software framework for dynamic model emulation[END_REF], [START_REF] Bencala | Simulation of solute transport in a mountain pool-and-rifle stream: a transient storage model[END_REF]. The ADE takes the following form:

             ∂χ(t, x) ∂t = a 02 (ρ(t)) ∂ 2 χ(t, x) ∂x 2 -a 01 (ρ(t)) ∂χ(t, x) ∂x +b 00 (ρ(t))u(t, x) χ(t = 0, x) = 0 y(t k , x ℓ ) = χ(t k , x ℓ ) + e(t k , x ℓ )
(33) The noise-free signal χ(x ℓ , t k ) (Fig. 1) represents the solute propagation in the river and is obtained by discretization of the partial differential equation by a finite difference method (see e.g. [START_REF] Morton | Numerical Solution of Partial Differential Equations[END_REF]).

The temporal sampling is T s = 10s and N = 8640 temporal points are measured. The spatial sampling is h = 10m for L = 61 space points. These 61 points can represent 61 sensors distributed in the river. u(t, x) describes the source of the pollutant, and is represented by a PRBS signal at the space point x = 20m. In realistic conditions of simulation the a 02 (ρ(t)) diffusion and the a 01 (ρ(t)) advection coefficients are dependent of the river flow Q(t) (Fig. 1). The scheduling variable ρ is then taken as ρ(t) = Q(t) and the coefficients a 02 (ρ(t)), a 01 (ρ(t)) and b 00 (ρ(t)) are given by

     å02 (ρ(t)) = 532 + 7.3ρ(t) (34a) å01 (ρ(t)) = 52 + 0.72ρ(t) (34b) b00 (ρ(t)) = 1 (34c)
For model identification purposes, the output signal is supposed to be corrupted by a two dimensional, zero-mean and normally distributed, discrete-time noise signal, with a Signal-to-Noise Ratio (SNR) defined as

SN R = 10log 10 P χ P e (35)
where P e and P χ are the average power of the noise and deterministic output.

To provide representative results, Monte Carlo simulation analysis for 50 runs are used to illustrate the performance of the approach. Table I and II display the estimated parameters when the additive noise is white and colored respectively. In the colored case, the additive colored noise v(t k , x ℓ ) is done by:

v(t k , x ℓ ) = 1 1 -1.2q -1 t + 0.4q -2 t e(t k , x ℓ ). (36)
The parameters are estimated by the proposed approach denoted as SRIVC and the LPV-SILS-PDE developed in [START_REF] Schorsch | Identification of LPV partial differential equation models[END_REF] and denoted as SILS (for simplified iterative least squares).

From both Tables, it can be noted that the estimates of the proposed SRIVC method are unbiased with relatively small standard deviations. As expected, the estimates delivered by the SILS method are unbiased only when the output is corrupted by a white noise.

VII. CONCLUSION

A refined instrumental variable algorithm has been presented to estimate linear parameter-varying partial derivative equations. The implied continuous filtering operations are rewritten under a discrete transfer function form. Based on a numerical example, the proposed method delivered unbiased estimated parameters whatever the noise structure while the iterative least squares method gave biased estimates when the additive noise was colored. Future research investigates the case where few sensors are available for the reconstruction of the noise-free output.

Fig. 1 .

 1 Fig. 1. River flow ρ(t) and noise-free output χ(t, x)

TABLE I MC

 I SIMULATION RESULTS, WHITE NOISE, SN R = 25dB â0

			02	â1 02	â0 01	â1 01	b0 00
	method	True value	532	7.3	52	0.72	1
	SILS	mean	534.2	7.263	51.85 0.7263	1.001
	iter= 7	std	1.2	0.042	0.08	0.0035	0.001
	SRIVC	mean	533.4	7.321	51.98 0.7207	1.001
	iter= 8	std	1.2	0.0419	0.08	0.0035	0.001
			TABLE II
	MC SIMULATION RESULTS, COLORED NOISE, SN R = 25dB
			â0 02	â1 02	â0 01	â1 01	b0 00
	method	True value	532	7.3	52	0.72	1
	SILS	mean	392.6 9.7476	58.68 0.4313	0.937
	iter= 7	std	4.9	0.1729	0.34	0.0148	0.003
	SRIVC	mean	530.8 7.2991	52.02 0.7184	0.999
	iter= 8	std	5.3	0.1811	0.35	0.0154	0.004

The authors are with University of Lorraine, CRAN, UMR 7039,

rue Jean Lamour, F-54519 Vandoeuvre-les-Nancy, France; 2 CNRS, CRAN, UMR 7039, France Email adresses: julien.schorsch@univ-lorraine.fr; vincent.laurain@univlorraine.fr; marion.gilson@univ-lorraine.fr; hugues.garnier@univ-lorraine.fr