Brahim Belaoucha

Mouloud Kachouane

Théodore Papadopoulo

Multivariate Autoregressive Model Constrained by Anatomical Connectivity to Reconstruct Focal Sources

In this paper, we present a framework to reconstruct spatially localized sources from Magnetoencephalography (MEG)/Electroencephalography (EEG) using spatiotemporal constraint. The source dynamics are represented by a Multivariate Autoregressive (MAR) model whose matrix elements are constrained by the anatomical connectivity obtained from diffusion Magnetic Resonance Imaging (dMRI). The framework assumes that the whole brain dynamic follows a constant MAR model in a time window of interest. The source activations and the MAR model parameters are estimated iteratively. We could confirm the accuracy of the framework using simulation experiments in both high and low noise levels. The proposed framework outperforms the two-stage approach.

I. INTRODUCTION

EEG and MEG are two non-invasive modalities that allow us to measure brain activity with high temporal resolution. Obtaining the distributed sources activation from these measurements is an under-determined problem due to the small number of measurements with respect to the number of sources. To obtain a unique solution, different constraints can be added [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF]. These can be divided into three different categories: spatial, temporal, and spatiotemporal constraints.

There are several studies that show that the reconstruction improves when including spatiotemporal constraints [START_REF] Lamus | A spatiotemporal dynamic distributed solution to the MEG inverse problem[END_REF], [START_REF] Fukushima | MEG source reconstruction constrained by diffusion MRI based whole brain dynamical model[END_REF]. In this work, we assume that the source signals follow a multivariate autoregressive (MAR) model. Fukushima et al. [START_REF] Fukushima | MEG source reconstruction constrained by diffusion MRI based whole brain dynamical model[END_REF] proposed a MAR model in which matrix elements depend on the anatomical connectivity, but their method is time-consuming and hard to implement. Other studies, [START_REF] Giraldo | Estimation of dynamic neural activity using a kalman filter approach based on physiological models[END_REF], [START_REF] Long | Large scale kalman filtering solutions to the electrophysiological source localization problem-a meg case study[END_REF], assume that sources increase or decrease by the same factor or interact with only their direct neighbors. These assumptions significantly reduce the degrees of freedom of the dynamics. We provide more degrees of freedom by allowing source interactions constrained by anatomical connectivity obtained from dMRI information. This allows us to generalize our work to the whole brain. Focal reconstructions are important because they are consistent with the notion of functional specialization of small cortical regions. Focal Underdetermined System Solver (FOCUSS) [START_REF] Gorodnitsky | Neuromagnetic source imaging with FO-CUSS: a recursive weighted minimum norm algorithm[END_REF] is used to detect spatially localized sources. It is applied to each time sample and tends to set to zero sources with low activation at each time sample. This makes it unsuitable for MAR parameter estimation. Our framework estimates the source activation and the MAR parameters iteratively. This approach corrects the false activation and functional connections that 1

Université

Côte d'Azur, Inria, Nice, France belaoucha.brahim@etu.unice.fr brahim.belaoucha@inria.fr 2 USTHB, LRPE, ParIMed, Algiers, Algeria can be obtained by two-stage approaches when working with low signal-to-noise ratio (SNR).

We have developed a framework constrained by a temporal dynamic over the whole brain to reconstruct sources. The source dynamics represented by the MAR model is based on the underlying anatomical connections obtained from dMRI information. Through simulations at different noise levels, we compare the accuracy of the proposed framework in the source reconstruction and MAR model parameters estimation with the ones obtained from a non-dynamical sparse solution based on mixed-norm estimates (MxNE) [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. Our framework uses the MxNE estimates to initialize the source intensities then estimates the MAR entries iteratively.

II. DATA ACQUISITION AND PROCESSING

Structural and diffusion MRI data were taken from Wakeman et al. [START_REF] Wakeman | A multi-subject, multi-modal human neuroimaging dataset[END_REF]. Cambridge University Psychological Ethics Committee approved all experimental procedures involving human subjects [START_REF] Wakeman | A multi-subject, multi-modal human neuroimaging dataset[END_REF]. The T1 weighted images of size 256 × 256×192 were acquired by a Siemens 3T Trio with GRAPPA 3D MPRAGE sequence (TR = 2250 ms; TE = 2.99 ms; flip-angle = 9 • ; acceleration factor = 2) at 1 mm isotropic resolution. The diffusion weighted images of size 96×96×68 were collected by the same scanner at 2 mm isotropic resolution (64 gradient directions and b-value =1000 s/mm 2 ), with one b0 image. The cortical surface was extracted using Freesurfer [START_REF] Dale | Cortical surface-based analysis I: Segmentation and surface reconstruction[END_REF] from T1 and remeshed to 2003 vertices (sources). The transformation between the anatomical and diffusion space was obtained by registering the brain in the two spaces using FSL [START_REF] Smith | Advances in functional and structural MR image analysis and implementation as FSL[END_REF]. The MEG/EEG forward problem, based on Boundary Element Method (BEM), is obtained using OpenMEEG [START_REF] Gramfort | OpenMEEG: opensource software for quasistatic bioelectromagnetics[END_REF], [START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF].

III. STRUCTURAL CONNECTIVITY

The structural connectivity of the sources (nodes on the mesh) is obtained from the dMRI information. Each node in the cortical mesh is considered as a seed. For each seed, we compute the connectivity profile, i.e. its connectivity with every voxel in the diffusion image space. Connectivity profiles are obtained by running probabilistic tractography for all seeds. We assume a model of multiple fiber orientations for each voxel developed in FSL [START_REF] Belaoucha | Cortical surface parcellation via dMRI using mutual nearest neighbor condition[END_REF] with 5 × 10 3 samples. The connectivity between node i and j, c ij , is the number of samples that were drawn from i and arrived at j divided by the total number of samples. CM is the connectivity matrix with elements CM (i, j) = max(c ij , c ji ). The source i and j are considered to be structurally connected if CM (i, j) ≥ c th . We set c th to 0.1 to remove weak connections. We obtained 3935 anatomical connections.

IV. ALGORITHM

The EEG/MEG measurements (M ∈ R n×T ) are related to the brain sources (J ∈ R m×T ) by the following linear relationship;

M = GJ + ,
where n and m are the number of sensors and sources respectively, T is the length of the time window, is Gaussian observation noise and G ∈ R n×m is the lead field matrix. We assume that the sources (J t ) at time t are related to their past values by:

J t = AJ t-1 + ω t ,
where A is the MAR matrix and ω t is Gaussian generator noise. The elements of the MAR matrix are constrained by the structural connectivity obtained from dMRI information. The only non-zero elements of A represent anatomically connected source pairs including self-connections. The dynamics of the model depend on the eigenvalues of A. The sources can be obtained by minimizing the following functional;

E(J) = 1 2T M -GJ 2 2 + λU (J) , (1)
where U is some prior on sources and λ is a positive regularization parameter. In focal activation, most of the sources are inactive. This sparsity can be obtained by putting (J) = J 21 , where J 21 = i t J t (i) 2 . This norm, l 12 , induces sparsity over the time window i.e. the sites of active sources will be the same over the entire time window of interest. The solution of Eq. ( 1) is called the mixednorm estimate (MxNE). More details about the MxNE can be found in [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. To enhance the solution of the mixed-norm solution in low SNR values and to estimate the functional connectivity which is represented by the MAR model, we introduce an iterative framework. The MAR matrix, A, is estimated over a time window of interest by minimizing D:

D(A) = J 1, T -AJ 0, T -1 2 2 ,
where J 1, T and J 0, T -1 are the values of the sources in a time window of T samples. This can be rewritten as;

D(A) = J 1, T -J * A 2 2 , (2)
where A = vec(A) and J 1, T = vec(J 1, T ) and

J * =       J 0, T -1 0 T ×m • • • 0 T ×m 0 T ×m J 0, T -1 • • • 0 T ×m . . . . . . . . . . . . 0 T ×m 0 T ×m • • • J 0, T -1      
, where 0 T ×m is a T × m zero matrix and (• ) is the transpose operator. The solution of Eq. ( 2) is;

A = (J * J * ) -1 J * J 1,T , (3)
Because of the block of zeros in J * , it is possible to estimate separately each row A i of A:

A i = (J 0, T -1 J 0, T -1 ) -1 J 0, T -1 J 1, T (i) i ∈ {1, .
., m}, where J 1, T (i) is the i th row of J 1, T . To estimate the source intensities and their interaction iteratively, we replace J 1, T by AJ 0, T -1 in the data fit term of Eq. ( 1) and U (J) by U (J 0, T -1 ). Eq. ( 1) becomes:

E(J 0, T -1 ) = 1 2T M 1, T -GAJ 0, T -1 2 2 + λ J 0, T -1 21 (4)
The framework, iterative Source and Dynamics Reconstruction (iSDR), is explained in Algorithm 1. We start by initializing the MAR matrix A with the identity (I), J with the MxNE solution and the maximum number of iterations N . In step 4, we estimate the sources by optimizing Eq. ( 1). Then A is updated by Eq. ( 3) in steps 6-8. The algorithm stops iterating if there is a small change in the data fit error, RSS = M 1,T -GAJ 0, T -1 2 , or after N iterations. The difference between Eq. ( 4) and Eq. (1) (MxNE) is that the lead field matrix, G, is multiplied by the MAR matrix, A. The columns and rows of A corresponding to inactive sources are set to zero in step 7. This allows the framework to reduce the effect of inactive sources on the measurements, i.e. reduce the number of unknowns resulting in better estimation of the source activation. JT ← Ã JT -1

6:

for j=1:m do Update à 7:

à j ← (J 0,T -1 J 0,T -1 ) -1 J 0,T -1 J 1,T (j) 8:
end for 9:

RSS i ← M 1,T -G Ã J 2 10: if RSS i -RSS i-1 < ε then 11: break 12:
end if

13:

end for 14:

return J, Ã 15: end procedure V. RESULTS AND DISCUSSION We used a realistic three layers head model. 70 EEG electrodes were considered. To test the framework, we first obtain the structural connectivity from the dMRI as explained in section III. Then, we investigate the parameter estimation accuracy of the proposed framework by using simulated data generated from a simple MAR model. Four sources (S1, S2, S3, and S4) are located in the Occipital lobe, see Fig. 1 (a) and connected anatomically, see Fig. 1 (b). We activated S1 and S2 while S3, S4, and the other sources are inactive. We decided not to activate S3 and S4 to investigate the effect of the assumption that the sources functionally interact if they are anatomically connected. S1 and S3 are anatomically connected but functionally independent.

J t (S1) = 0.96J t-1 (S1) + 0.25J t-1 (S2) J t (S2) = -0.25J t-1 (S1) + 0.95J t-1 (S2)
with J 0 (S1) = -J 0 (S2) =5 nAm (nano Ampere), the simulated sources can be seen in Fig. 1 (c). The MAR matrix has two non-zero eigenvalues 0.955 ± 0.25j. The regularization parameter, λ, is fixed by a 3-fold cross validation [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].

We simulate sources in low and high SNR. The simulated measurements were obtained from the simulated sources (Fig. 1 (c)) multiplied by the lead field matrix, G, and we added a Gaussian observation noise [SNR: 40, ±5 dB] of zero mean. The source reconstruction of MxNE can be seen in Fig. 2 (a) to (c). Fig. 2 (d) shows the estimated eigenvalues using two-stage approach (MxNE then Eq. ( 3)) at SNR equal to 5 dB. Fig. 2 (e) to (g) represent the result of the proposed method at different noise levels.

Fig. 2 (h) to (j) show the true, in black, and the estimated, in red, eigenvalues by the proposed framework of the MAR matrix in polar coordinates at different noise levels and for 100 runs. The majority of the eigenvalues are estimated to be zero, this can be guessed clearly from the reconstructions of the iterative approach at the different noise levels, Fig. 2 (e) to (g), in which the framework detected the true inactive sources (the corresponding MAR coefficients equal to zero). Fig. 2 (k) and (l) represent the mean over 100 runs of the estimated functional interactions between the four sources at SNR equal to ±5 dB.

The reconstruction error between the estimated sources and the ground truth, J -J shown in Table 1. The mean error is computed over 100 runs at the different noise levels. The iSDR enhances significantly the mixed-norm estimates solution especially at low SNR settings, compare Fig. 2 (b) with (f) and (c) with (g) and see the last column of Table 1.

The iSDR provides accurate results because MxNE sets sources to zero in the whole time window of interest. This results in rows and columns of zeros in the MAR matrix, step 7 in Algorithm. In the next iteration, A, with the zero columns and rows, is multiplied by G. This is equivalent to removing some generators from the source space resulting in a better approximation. The number of iterations, N , is inversely proportional to noise level i.e. high N value for low SNR.

The two-stage approach found false-positive active sources in low SNR, Fig. 2 (b) and (c), which result in false functional connections. This is the limitation of two-stage approximation of the MAR model. The iterative framework could correct the false activation, see Fig. 2 (f) and (g). It could detect the correct connections between the sources for SNR ≥ 5 dB (Fig. 2 (k)) and it underestimates the MAR coefficients for 5 ≥ SNR ≥ -5 dB (Fig. 2 (l)) because of the noise but they are still close to the ground truth.

In this paper, we present an iterative framework for EEG and MEG source reconstruction constrained by a whole brain dynamics model based on dMRI information. In our simulation study, the proposed framework correctly reconstructed sources and estimated the MAR elements, while the non-dynamical two-stage approach based on MxNE solution detects false-positive sources resulting in false functional connections. The framework assumes that the source dynamics is constant over a time window of T samples. The optimal number of samples needed to estimate the MAR entries is proportional to the maximal number of anatomical connections per source. It is hard for the framework to detect sources that have small activation in time and duration. Like all initialization dependent approaches, the initialization has to be within some proximity to the true sources, in order to find the correct activation. This means that as long as the initial active set contains true active sources, we are able to recover correct activations and interactions.

Representing the source dynamics by linear autoregressive modeling would introduce errors to the estimated functional interactions for sources with non-linear relationships. In this work, we dealt with focal sources. Extended generators (constant activation per region) can be assumed on regions obtained from dMRI-based parcellation [START_REF] Fukushima | MEG source reconstruction constrained by diffusion MRI based whole brain dynamical model[END_REF], [START_REF] Belaoucha | Cortical surface parcellation via dMRI using mutual nearest neighbor condition[END_REF] in the proposed framework. Higher-order MAR model can be assumed.

Algorithm 1 iterative 3 :

 13 Source and Dynamics Reconstruction (iSDR) 1: procedure iSDR(G, M, N ) 2: Initialize: Ã ← I, (n, m) ← shape(G), ε for i=1:N do 4: J0,T -1 ← arg min J 0,T -1 (E) Estimate J 5:

Fig. 1 :

 1 Fig. 1: Location of the simulated sources on the inflated cortical surface (a) and their anatomical connections (b), only the connections between the four sources. The magnitude of the simulated sources (c), sampling frequency equal to 1 KHz.

Fig. 2 :

 2 Fig. 2: The source reconstruction using mixed-norm estimates at SNR equal to (a) 40 dB, (b) 5dB and (c) -5dB. The estimated eigenvalues of à using two-stage approach , i.e MxNE then Eq. (3) in (d). The source reconstruction obtained by the proposed method (iSDR) at SNR (e) 40 dB result after 4 iterations, (f) 5 dB result after 10 iterations and (g) -5 dB result after 14 iterations. The estimated eigenvalues (100 runs) by the proposed iterative method at different noise levels, (h) 40 dB, (i) 5 dB and (j) -5 dB. In (k) and (l), we show the mean over 100 runs of the estimated functional interactions between S1, S2, S3 and S4 for SNR equal to 5 and -5 dB respectively.

TABLE I :

 I The mean and standard deviation (std) of the reconstruction error over 100 runs at different noise levels (SNR ={40,5,-5} dB). The values are in nAm.

VI. CONCLUSION

We developed a source reconstruction framework constrained by a whole brain dynamical model based on dMRI information. The framework iteratively solves the source activation and MAR entries. In our study, the framework could correct the false activation in low SNR level generated by the two-stage approach and detect the correct functional interactions between sources. Application to real datasets and comparison of the solutions to the physiological literature is essential to validate this framework.


Said Kouachi
email: koashy@qu.edu.sa

A balance law for some strongly coupled reaction-di¤usion systems and an invariance Theorem

Keywords: AMS Classi…cation: 35K57, 35K45 Reaction di¤usion systems, Invariant regions, Global existence

The purpose of this paper is to prove an invariance theorem for m-components strongly coupled reaction-di¤usion systems with non-homogenous boundary conditions. A characteristics of these systems have been de…ned and a new balance law has been constructed. Global existence in time of solutions is deduced without conditions on the nonlinearities growth.

INTRODUCTION

We consider the reaction-di¤usion system @u i (t; x) @t m j=1 a ij u j (t; x) = f i (t; x; u); in R + ;

(1) with the boundary conditions @u i (t; x) @ = g i (x; u) on 1 and u i (t; x) = i (x) on 2 :

(2) and the initial data

u i (0; x) = u 0 i (x); in ; (3)
where is an open bounded domain of class C 1 in R n with bound-

ary @ = 1 [ 2 is a disjoint union with 1 = m i=1 1 i ; 2 = m i=1 2 i
and u = (u 1 ; :::; u m ) : The di¤usion matrix A = (a ij ) 1 i;j m with real entries, is supposed to be diagonalizable with positive eigenvalues: 0 < 1 < 2 < ::: < l ; each i is of multiplicity i 1 . The reaction term f = (f 1 ; :::; f m ) is a continuous function from R + R m into R m , the functions g = (g 1 ; :::; g m ) and = ( 1 ; :::; m ) on the boundary are continuous from @ R m into R m and @ @ denotes the outward normal derivative on @ . Systems of this form can be interpreted as models in several chemical, biological and dynamics of populations models where the quantity u i stands for the chemical or molecular concentration and population density of the i-th specie respectively(see [START_REF] Garcia-Ybarra | Cross transport e¤ects in premixed ‡ames[END_REF], [START_REF] De Groot | Non-Equilibrium Thermodynamics[END_REF], [START_REF] Jorne | The di¤usion Lotka-Volterra oscillating system[END_REF], [START_REF] Lee | On the general linear coupled system for di¤usion in media with two di¤usivities[END_REF] and [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF]). The technique based on the bounded invariant regions is among those used to solve the problem of global existence in time of strong solutions for systems such as these and their quasilinear generalizations and accordingly it has received much attention.

There exists a wide bibliography on invariant sets for nonlinear parabolic and elliptic systems. The reader is referred to the expository articles [START_REF] Amann | Invariant sets and existence theorems for semilinear parabolic and elliptic systems[END_REF], [START_REF] Chueh | Positively invariant regions for systems of nonlinear di¤usion equations[END_REF], [START_REF] Kuiper | Invariant sets for elliptic and parabolic systems[END_REF], [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF], [START_REF] Redlinger | Invariant sets for strongly coupled reaction di¤usion systems under general boundary conditions[END_REF], [START_REF] Smoller | Shock Waves and Reaction-Di¤usion Equations[END_REF], [25], [START_REF] Kouachi | Invariant regions and global existence of solutions for reaction di¤usion systems with a full matrix of di¤usion coe¢ cients and no homogeneous boundary conditions[END_REF], [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of Reaction-Di¤usion Equations with a balance law[END_REF], [START_REF] Kanel | Global solutions to a reaction di¤usion system[END_REF] and references there). Whereas all authors mentioned thus far have considered the balance law condition only for tow components systems, here we de…ne a characteristics and introduce a new balance law for the general system [START_REF] Amann | Invariant sets and existence theorems for semilinear parabolic and elliptic systems[END_REF]. Then we prove an invariance theorem which gives global existence of solutions when the constructed invariant sets are bounded. This work was motivated by the interesting results of [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] who studied the two component strongly coupled system on the form (1) satisfying the following balance law conditions 2 (x) + ( 1 a 22 ) jf 1 (t; x; u)j a 12 jf 2 (t; x; u)j

( 2 a 22 ) jf 1 (t; x; u)j + 1 (x) ; and (x) + ( 1 a 22 ) jg 1 (x; u)j a 12 jg 2 (x; u)j

( 2 a 22 ) jg 1 (x; u)j + (x) ;

(4) on R + R 2 and @ R 2 respectively, where (x) and (x) are measurable functions on and @ respectively. He proved global existence of solutions via bounded invariant regions technique under some conditions on the di¤usion matrix A (with positive distinct eigenvalues and a 12 > 0) and the following dissipative conditions u 1 f 1 (t; x; u) 0 and u 1 f 2 (t; x; u) 0; u 1 g 1 (x; u) 0 and u 1 g 2 (x; u) 0;

(5) on R + R 2 and @ R 2 respectively. As [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] when m = 2, by these results we generalize the works of [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of Reaction-Di¤usion Equations with a balance law[END_REF] and [START_REF] Kanel | Global solutions to a reaction di¤usion system[END_REF] to more general di¤usion matrices by eliminating the need for the di¢ cult-to-establish bound on the di¤usion terms: Also we eliminate for the reaction terms that need to satisfy a polynomial growth condition and a strict balance law. Following [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of Reaction-Di¤usion Equations with a balance law[END_REF], [START_REF] Kanel | Global solutions to a reaction di¤usion system[END_REF] and [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] we suppose that u p f k ; k = 1; :::; m don't change sign for some …xed p = 1; :::; m not on the hole space R n but for ju p j su¢ ciently large. Simultaneously permuting unknowns, reactions and rows-columns of the di¤usion matrix, if necessary, we may assume without lost of generality p = 1:

The paper is organized as follows. After clarifying the notations and stating the various de…nitions in the next section, we present the results and proofs in the third section where we begin by treating the homogenous case in the …rst subsection. In the second one we treat the non-homogenous case. The last section will be consecrated to some applications for the two and three components reaction di¤usion systems.

De…nitions and notations

In this section we present some necessary notations and de…nitions which are needed in our subsequent sections. For technical reasons due to the di¢ cult calculus of the explicit eigenvalues, namely when m 3, we suppose the di¤usion matrix A to be symmetric or at last Symmetrizable, then there exists a diagonal matrix S with positive entries (called Symmetrizer), such that the matrix S:A is symmetric and the system can be transformed by the change variable: v = S 1 2 u to an equivalent system with symmetric matrix S 1 2 AS 1 2 . Note that, from the Cauchy interlace Theorem for m m symmetrizable matrices, the (m 1) order principal minors have the same sign. As examples of such matrices we cite sign symmetric (i.e. a ij :a ji > 0 and a ij = 0 ) a ji = 0; 0 i; j m) Tridiagonal, Pentadiagonal, Heptadiagonal and particularly, Symmetric matrices. The case m = 2 does not require the symmetrizability of the di¤usion matrix and will be treated alone as a remark at the end of the next section. Since the di¤usion matrix is supposed to be diagonalizable, then it possess a full set of eigenvectors. If we denote by x i = x i1 ;

; x im t an eigenvector of the matrix A t (left eigenvector of A) associated with its eigenvalue i and X the m m matrix formed with the rows x i ; the dissipativity condition on the reactions (i.e. u 1 f k don't change sign, k = 1; :::; m) requires x i1 6 = 0; i = 1; :::; m ( i.e. the …rst column of the matrix X does not contain zeros). This condition is satis…ed by those eigenvectors associated to the simple eigenvalues i since the rank of the matrix A t i I which is also symmetrizable, is equal to m 1 and then, all its principal minors of order m 1 have the same sign and consequently non null, particularly those obtained by deleting the …rst row-column pair. When the eigenvalue i is of multiplicity i , the rank of the matrix A t i I m is equal to m i , then at last one of its principal minor of order m i is non null. Let us choose i integers 1 i 1 < ::: < i i m such that i 6 = 0; where i = det (a ir; is ) r;s6 =1;:::

; i i I m i ; (6)
and where I m i denotes the unit diagonal matrix of order m i : Let V i = x i 1 ; x i 2 ; :::; x i i be the eigenspace associated to the multiple eigenvalue i , then for each j = 1; :::; i ; i is a simple eigenvalue of the i -order submatrix A (i j ) obtained from A t by deleting the ( i 1) row-column pairs i 1 ; :::; i j 1 ; i j+1 ; :::; i i

A (i j ) = 0 B B B @ a l 1 l 1 a l 2 l 1 a l i l 1 a l 1 l 2 a l 2 l 2 a l i l 2 . . . . . . . . . a l 1 l i a l 2 l i a l i l i 1 C C C A ; (7)
where i = m i + 1 and l 1 < l 2 < ::: < l r 1 < l r =: i j < l r+1 < ::: < l (m i +1) are di¤erent from i 1 ; :::; i j 1 ; i j+1 ; :::; i i . Since the …rst component of the jth eigenvector x i j associated to i is supposed to be non null, then

i j ;1 = det (a lslr ) s; r6 =1 i I m i =: a l 2 l 2 i a l i l 2 . . . . . . a l 2 l i a l i l i i 6 = 0: (8)
where i j ;1 denotes the determinant of the (m i )-order principal submatrix obtained from A (i j ) i I i by deleting its …rst row-column pair is non null and we can take in (7) l 1 = 1:

In terms of matrices, the system (1) can be written

@u @t A u = f (t; x; u); in R + ;
with the boundary conditions analogous to (2) @u @ = g (x; u) ; on 1 and u(t; x) = (x) ; on 2 ;

and the initial data u(0; x) = u 0 (x); in ;

where g = (g 1 ; :::; g m ) t and = ( 1 ; :::; m ) t . If we denote by D = diag f 1 ; :::; m g and make the change of variables z = Xu, this system is equivalent to the following diagonal system

@z @t D z = Xf (t; x; X 1 z); in R + ; (9)
with the boundary conditions

@z @ = Xg (x; X 1 z) ; on R + 1 ;
and

z(t; x) = X (x) ; on R + 2 ;
and initial data z(0; x) = Xu 0 (x); in :

In order to see things more clearly, we shall write the diagonal system on its scalar form. But before this, let us collect here brie ‡y some notations used frequently in the following sections. We shall denote for all j = 1; :::; i the components of the jth eigenvector x i j associated to the eigenvalue i of multiplicity i by x i j ;l k ; k = 1; :::; i , the others are zeros. When the eigenvalue is simple, for simplicity in notation, we shall consider it as of multiplicity i = 1, then i = m and the corresponding components of its associate eigenvector are denoted in the sequel by the same way: x i j ;l k ; k = 1; :::; i .

To each eigenvector x i j ; j = 1; :::; i associated to the eigenvalue i of multiplicity i , corresponds an equation

@z i j @t i z i j = F i j (t; x; z) ; in R + ; (10)
with the boundary conditions

@z i j (t; x) @ = G i j (x; z) on R + 1 ;
and

z i j (t; x) = i j (x) on R + 2 ; (11)
and the initial data

z i j (0; x) = z i j (0; x) = z 0 i j (x); in ; (12)
where

z i j = i k=1 x i j ;l k u l k ; (13)
F i j (t; x; z) = i k=1 x i j ;l k f l k (t; x; X 1 z) ; (14)
G i j (x; z) = i k=1 x i j ;l k g l k x; X 1 z ; (15)
and

i j (x) = i k=1 x i j ;l k l k (x) :
Let us de…ne, for each simple eigenvalue, the linear form i (u) to be the determinant of the matrix obtained from (A t i I m ) by replacing the …rst row with the corresponding components of u. i.e.

i (u) = u 1 u 2 u m a 12 a 22 i a m2 . . . . . . . . . a 1m a 2m a mm i : (16)
When the eigenvalue i is of multiplicity i , we shall prove the existence of i linearly independent forms i j (u) ; j = 1; :::; i analogous to i (u), each one depends on i components of u. Indeed, as i is a simple eigenvalue of the i -order submatrix A (i j ) given by [START_REF] Jorne | The di¤usion Lotka-Volterra oscillating system[END_REF], then analogously to the case when the eigenvalue is simple, for each j = 1; :::; i ; we de…ne the linear form i j (u) to be the determinant of the i -order matrix obtained from A (i j ) i I i by replacing its …rst row with the corresponding components of u: It can be written as follows

i j (u) = u l 1 u l 2 u l i a l 1 l 2 a l 2 l 2 i a l i l 2 . . . . . . . . . a l 1 l i a l 2 l i a l i l i i ; j = 1; :::; i ; (17)
where the coe¢ cient i j ;1 of u 1 =: u l 1 is given by [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of Reaction-Di¤usion Equations with a balance law[END_REF]. We note that in the case of simple eigenvalues the de…nitions ( 16) and ( 17) coincide.

In other words, for each j = 1; :::; i ; the linear form i j (u) analogous to ( 16) is the determinant of the i -order submatrix obtained from (A t i I m ) by replacing its …rst row by the corresponding components of u, then we delete the ( i 1) row-column pairs i 1 ; :::; i j 1 ; i j+1 ; :::; i i . Remark 1 Technically, the i forms i j (u) ; j = 1; :::; i ; can be obtained directly from the expression (16) of i (u) by replacing i with the multiple eigenvalue and deleting for each j; the ( i 1) row-column pairs i 1 ; :::; i j 1 ; i j+1 ; :::; i i : Remark 2 We note that, from the Cauchy interlace Theorem for symmetrizable matrices, all the determinants of the (m i )-order principal submatrices of A (i j )

i I i have the same sign as i given by ( 6) which is independent of j and of course, that given by [START_REF] Kanel | Pointwise a priori bounds for a strongly coupled system of Reaction-Di¤usion Equations with a balance law[END_REF].

In order to clarify the situation for the reader, we give the following example

Example 3 Take m = 4; 0 < 1 < 2 < 3 = 4 =
; i = 1; 2:
The two remained characteristics associated to the double eigenvalue are obtained from the above expression of i (u) for i = 3: one time by deleting its fourth row-column pair to get the …rst characteristic 2 3 (u). Another time we delete its second row-column pair to obtain the second characteristic 4 3 (u):

u 1 u 2 u 3 u 4 a 12 a 22
a 32 a 42 a 13 a 23 a 33 a 43 a 14 a 24 a 34 a 44

! 8 > > > > > > > > < > > > > > > > > : 2 3 (u) = u 1 u 2 u 3 a 12 a 22
a 32 a 13 a 23 a 33

; 4 3 (u) = u 1 u 3 u 4 a 13 a 33
a 43 a 14 a 34 a 44 Now we shall write the reactions [START_REF] Mahara | Three-variable reversible Gray-Scott model[END_REF] and the boundary functions (15) of the diagonal system (10) each one with respect to its corresponding linear form given by ( 17). As each eigenvector x i j ; j = 1; :::; i ; associated to the eigenvalue i of multiplicity i (even if i = 1) is parallel to the …xed eigenvector x i j 0 with components i j ;l k ; k = 1; :::; i ; given by

i j ;l k = ( 1) l k +1 det (a lslr ) s6 =k; r6 =1 i I m i ; (18)
and since the …rst component of any eigenvector in di¤erent from zero, then i j ;1 6 = 0: Consequently the new reaction ( 14) and the boundary functions ( 15) can be written as follows

F i j (t; x; z) = x i j ;1 i j ;1 i j (f (t; x; u)) ; j = 1; :::; i ; (19)
and

G i j (x; z) = x i j ;1 i j ;1 i j (g(x; u)) ; j = 1; :::; i ; (20)
respectively, where u = X 1 z. Recall that all the determinants i j ;1 ; j = 1; :::; i have the same sign depending only on i and not on j. We should also remark that the coe¢ cients of the u l k 's intervening in the expression of the i j (u)'s represent the components non null i j ;l k ; k = 1; :::; i given by ( 18) of the …xed eigenvector x i j 0 .

De…nition 4 A subset (L 1 ( )) m is called a positively invariant region (or more simply an invariant region) for system 1, if all solutions with initial data in remain in for all time in their interval of existence.

If there exists a bounded invariant set of system [START_REF] Kouachi | Invariant regions and global existence of solutions for reaction di¤usion systems with a full matrix of di¤usion coe¢ cients and no homogeneous boundary conditions[END_REF], then using the following well known existence result (see [START_REF] Friedman | Partial Di¤erential Equations of Parabolic Type[END_REF], [START_REF] Henry | Geometric Theory of Semilinear Parabolic Equations[END_REF], [START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Di¤erential Equations[END_REF], [START_REF] Pierre | Global Existence in Reaction-Di¤usion Systems with Control of Mass: a Survey[END_REF], [START_REF] Quittner | Superlinear Parabolic Problems (Blowup, Global Existence and Steady States)[END_REF] and [START_REF] Roth | Global Solutions of Reaction-Di¤usion Systems[END_REF]), the solution of the system (1) is global whenever u 0 is in the bounded subset =: X 1 ( ). 3 Results and proofs

The homogeneous boundary conditions case

In this subsection we consider the case when the Neumann boundary conditions in [START_REF] Chueh | Positively invariant regions for systems of nonlinear di¤usion equations[END_REF] are homogenous (i.e. g i (x; u) 0). The nonhomogeneous case will be deduced easily in the next subsection. We suppose that for ju 1 j su¢ ciently large, the reactions satisfy The dissipativity condition: Each of the functions u 1 f k (t; x; u) ; k = 1; :::; m, does not change sign on R + R m . That is

k :u 1 f k (t; x; u) 0; (21)
where k = 1 is a sign application from the subset of integers f1; 2; :::; mg and takes its values in f 1; +1g :

The balance law condition: The characteristics i j ( jf j) of the system are alternate functions with respect to i

( 1) i+1 i j ( jf j) 0; j = 1; :::; i ; i = 1; :::; l;

where jf j = ( 1 jf 1 j ; :::; m jf m j).

When the eigenvalue i is of multiplicity i > 1; the above condition means that the characteristics i j ( jf j) ; j = 1; :::; i don't change sign. For example, when k = ( 1) k+1 ; the assumptions ( 21) and ( 22) become u 1 f 1 (t; x; u) 0; f k (t; x; u) :f k+1 (t; x; u) 0; k = 1; :::; m 1;

and

( 1) i+1 : i j ( 1) k+1 jf k j 0; i = 1; :::; l;

respectively. Let us de…ne the set

= l \ p=1 p \ q=1 pq ; ( 23
)
where pq represents the rectangle

pq = n z 2 R m : pq z pq pq o ; (24)
p = 1; :::; l; q = 1; :::; p ; with edges

i j i j = z 2 pq : z i j = i j ; p q 6 = i j ; (25)
and

i j i j = n z 2 pq : z i j = i j o ; p q 6 = i j ; (26)
i = 1; :::; l; j = 1; :::; i : Our main result on the invariant regions and global existence of system (1) is the following Theorem 6 Suppose that the di¤usion matrix is symmetrizable, then under conditions ( 21) and ( 22), the region = X 1 ( ) with de…ned by ( 23) is invariant for system [START_REF] Amann | Invariant sets and existence theorems for semilinear parabolic and elliptic systems[END_REF]. Moreover the solution is global and uniformly bounded on for any initial data in L 1 ( ) : Furthermore when the conditions are satis…ed for all u 1 2 R; then

kuk 1 C u 0 1 (27)
where C is a positive constant depending only on the di¤usion matrix A and equal to the unity when A is symmetric.

Proof. Since the matrix A is diagonalizable with positive eigenvalues, the two systems (1) and ( 9) are both parabolic and equivalent. Then if the diagonal system has an invariant rectangle ; the original has an invariant region X 1 ( ). But, thanks to the invariant region's method (see H. j. Kuiper [START_REF] Kuiper | Invariant sets for elliptic and parabolic systems[END_REF] and [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] and J. A. Smoller [START_REF] Smoller | Shock Waves and Reaction-Di¤usion Equations[END_REF]), the region given by ( 23) is invariant for system (9) under the following conditions F i j (t; x; z) 0; , for all z 2 i j i j ; and F i j (t; x; z) 0, for all z 2 i j i j ;

(28) for all i = 1; :::; l; j = 1; :::; i : Since A is symmetrizable, then the system can be transformed to an equivalent system with symmetric diffusion matrix. Therefore, we can suppose that A is symmetric and then the matrix X formed with the rows x i j can be orthogonalized (i.e. chosen such that X 1 = X t ). Consequently, for all i = 1; :::; l; j = 1; :::; i ; the formula

z i j = i k=1 x i j ;l k u l k ;
gives

u i j = l p=1 p q=1 x pq p;i j z pq (29)
Particularly, the …rst component of the unknown u can be written as follows

u 1 = x i j ;1 z i j + l p=1 p P q=1; pq6 =i j x pq;1 z pq ! (30)
where from (8) x i j ;1 6 = 0. Let us verify that the rectangle given by ( 24) is invariant under conditions ( 21) and [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF]. For this purpose, we should verify (28) for all i = 1; :::; l; j = 1; :::; i : To prove this, we need to verify two inequalities along each of the m edges of the rectangle given by ( 24): For z 2 i j i j ; we have z i j = i j : Taking i j negative and su¢ ciently large in absolute value, then the summation

l p=1 p P q=1; pq6 =i j x pq;1 z pq ! , in
the above expression (30) of u 1 is bounded ( pq z pq pq ; p q 6 = i j ). This gives u 1 is also large in absolute value, but its sign depends on that of x i j ;1 and so we have two cases: When x i j ;1 > 0, then u 1 is negative which gives from ( 21) f k = k jf k j ; k = 1; :::m: Consequently, from ( 19) F i j (t; x; z) 0 under the following condition i j ( jf j)

i j ;1 0: (31)
As A is symmetric, then from the well known Cauchy's interlace Theorem, the eigenvalues of A and those of each of its principal submatrices of order m 1 interlace. That is

( 1) i+1 i j ;1 0; i = 1; :::; l: (32) Thus, [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF] gives F i j (t; x; z) 0 for all z 2 i j i j : When x i j ;1 < 0, then u 1 is positive and again from ( 21), we have f k = k jf k j ; k = 1; :::m. This gives

F i j (t; x; z) = x i j ;1 i j ( jf j) i j ;1 : (33)
Since x i j ;1 < 0, then using (32), we can conclude the positivity of F i j (t; x; z) on i j i j under the condition [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF].

Following the same reasoning for z 2 i j i j by choosing i j positive and su¢ ciently large, we can deduce that u 1 is also large in absolute value. Analogously, we have two cases: When x i j ;1 < 0, then u 1 is negative and this gives from ( 21) f k = k jf k j ; k = 1; :::m: Thus F i j (t; x; z) can be written on the form [START_REF] Quittner | Superlinear Parabolic Problems (Blowup, Global Existence and Steady States)[END_REF]. Using ( 31) and (32), we can say that [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF] gives F i j (t; x; z) 0 for all z 2 i j i j : Finally when x i j ;1 > 0, then u 1 is positive, also from ( 21), we have f k = k jf k j ; k = 1; :::m and then (33). Another time with (31) and (32) together, we can conclude that condition [START_REF] Savchik | Application of moments to the general linear multicomponent reaction-di¤usion equations[END_REF] gives F i j (t; x; z) 0 for all z 2 i j i j : This ends the proof of the …rst part of the Theorem when the di¤usion matrix is symmetric. When A is symmetrizable, the di¤usion matrix of the new system becomes B = S 1 2 AS 1 2 which is symmetric and the reaction is S 1 2 f . If we denote by A i j (u) the linear form i j (u) given by ( 16) and associated to the di¤usion matrix A; then B i j S 1 2 f = p s 1 A i j ( jf j) and i j ;1 ; j = 1; :::; i , remains the same, for all i = 1; :::; l. We follow the same steps of the proof by replacing u 1 and f by p s 1 u 1 and S 1 2 f respectively, where S = Diag(s 1 ; :::; s m ): Since is invariant (i.e. kzk kz 0 k), the uniform bounds (27) of the solutions are a trivial consequence of the relations ( 13) which, when written on the matricial form z t = Xu t , give

S 1 2 u X 1 z =: kzk = and kz 0 k XS 1 2 u 0 =: S 1 2 u 0 :
This ends the proof of the Theorem.

Remark 7 When all the eigenvalues of the di¤usion matrix are simple, then P ( i ) = i ( jf j) is a polynomial of degree m 1 with respect to i ; with leading coe¢ cient 1 jf 1 j : But from ( 22) P ( i ) changes sign m times and P ( 1 ) 0: Consequently we should take 1 = 1 (i.e.

u 1 f 1 0 on R + R m ).
Remark 8 When the balance law is strict (i.e. jf k j = jf 1 j ; k = 1; :::; m), the conditions ( 22) are independent of the reactions and become a simple conditions on the di¤usion terms.

The non-homogeneous boundary conditions case

When the Neumann boundary conditions are non-homogeneous, we consider, for ju 1 j is su¢ ciently large, the following hypothesis:

H1) The reaction terms satisfy [START_REF] Roth | Global Solutions of Reaction-Di¤usion Systems[END_REF] and the following conditions analogous to (22):

i j ( jf j)

i j ;1 i (x) ; i = 1; :::; l; j = 1; :::

; i ; (34)
on R + R m ; where (x) is a measurable function.

H2) The functions g i satisfy a conditions analogous to ( 21) and (34):

k :u 1 g k (x; u) 0; k = 1; :::; m; (35)
and i j ( jgj)

i j ;1 (x) ; i = 1; :::; l; j = 1; :::; i ;

on @ R m ; where (x) is a measurable function and is a sign application from the subset of integers f1; 2; :::; mg and takes its values in f 1; +1g , not necessary equal to .

In order to construct invariant regions in this case, we use a given positive and bounded ' (x) on : The edges of the set , given by [START_REF] Smoller | Shock Waves and Reaction-Di¤usion Equations[END_REF] are translated, each one by a vector function of R m which components are equal to x i j ;1 ' (x) ; i = 1; :::; l; j = 1; :::; i : In other words we de…ne the parallelepiped

(') = l \ p=1 p \ q=1 pq (') ; (37)
where

pq (') = n z 2 R m : pq x pq;1 ' z pq pq + x pq;1 ' o ;
p = 1; :::; l; q = 1; :::; i ; with edges the rectangles

i j i j x i j ;1 ' = z 2 pq : z i j = i j
x i j ;1 '; p q 6 = i j ; and i j i j + x i j ;1 ' = n z 2 pq : z i j = i j + x i j ;1 '; p q 6 = i j ; o j = 1; :::; i ; i = 1; :::; l: Following [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] in the case m = 2, we choose the function ' as a strong solution of the following linear elliptic problem

' (x) ; @' @ (x) : (38)
H3) The above elliptic problem has a solution ' two times continuously di¤erentiable on with …rst and second derivatives bounded on : In the case of Neumann boundary conditions on all @ ; it is of course necessary to satisfy the compatibility condition

Z (x) dx Z @ (x) d :
In order to state our results in this case, we apply the following standard proposition which can be found in the wide bibliography cited at the beginning of our introduction. Proposition 9 Suppose that, for some bounded ' on ; we have, for all i = 1; :::; l; j = 1; :::; i ; i z i j F i j (t; x; z); 8z 2 i j i j x i j ;1 ' ;

i z i j F i j (t; x; z); 8z 2 i j i j i + x i j ;1 ' ; (39)
in R + and @z i j @ G i j (x; z) ; 8z 2 i j i j x i j ;1 ' ;

@z i j @ G i j (x; z) ; 8z 2 i j i j i + x i j ;1 ' ; (40)
in R + @ ; then the region (') given by ( 37) is invariant for the diagonal system [START_REF] Kouachi | Invariant regions and global existence of solutions for reaction di¤usion systems with a full matrix of di¤usion coe¢ cients and no homogeneous boundary conditions[END_REF].

Our main result when the boundary conditions are non-homogeneous is the following Theorem 10 Suppose that hypotheses H 1 -H 3 are satis…ed, then the set = X 1 ( (')) ; where (') is given by (37) is invariant for system (1) and for any bounded initial condition on , there exists a unique global solution of problem (1)-(3) uniformly bounded on : Furthermore when the conditions are satis…ed for all u 1 2 R; we have

kuk 1 C u 0 1 + k'k 1 ; (41)
where C is a positive constant depending only on the di¤usion matrix A and equal to the unity when A is symmetric.

Proof. The two systems (1) and ( 9) being equivalent, then by this, if we verify the two inequalities (39) for the reactions and two others (40) for the boundary along each of the m edges; the set (') is invariant and consequently, by a convenient choice of the constant i j and i j i , it can contains any bounded initial data. So, using (38), it remains to verify, for all i = 1; :::; l; j = 1; :::; i , the following inequalities F i j (t; x; z) i x i j ;1 (x) ; for all z 2 i j i j x i j ;1 ' ; and F i j (t; x; z) i x i j ;1 (x) , for all z 2 i j i j + x i j ;1 ' ;

vector x i = (x i1 ; x i1 ) t ; we have

z i = x i1 u 1 + x i2 u 2 = x i1 a 22 i [(a 22 i ) u 1 a 12 u 2 ] ;
and

F i (t; x; z) = x i1 a 22 i [(a 22 i ) f 1 a 12 f 2 ] ; i = 1; 2;
then

F i (t; x; z) = x i1 i (f 1 ; f 2 ) a 22 i ;
where i are given by ( 16). Theorem 10 is applicable with i;1 = a 22 i and we have the same conclusions even the di¤usion matrix is not symmetrizable.

Remark 12 Following [START_REF] Kuiper | Invariant sets for elliptic and parabolic systems[END_REF] in the case of systems of two equations, Theorem 10 can be extended to the mixed boundary conditions or those of the type Dirichlet.

Remark 13

The operator can be replaced by any uniformly elliptic operator L.

Remark 14 Theorem 10 is applicable even when u 1 f k and u 1 g k change sign for u 1 small and this is not the case in the results obtained by authors until now. It is, for example the case of the two component reaction di¤usion system modeling the two reversible chemical reaction with nonlinear terms f 1 (u; v) = f 2 (u; v) = u v u p v q ; where the exponents are positive integers representing the number of molecules of each reactant.

Applications

Strongly coupled reaction di¤usion equations

Let us consider the reaction-di¤usion system (1) for m=2, with boundary conditions (2) and initial data (3). To get positive eigenvalues 0 < 1 < 2 ; we should assume that a 11 a 22 a 12 a 21 > 0;

(44) and a 11 + a 22 > 0:

(45) Suppose that

u 1 f 1 0, u 1 f 2 0, on R + R 2
and u 1 g 1 0, u 1 g 2 0; on @ R 2 ;

(46) with ju 1 j su¢ ciently large. This means that conditions [START_REF] Roth | Global Solutions of Reaction-Di¤usion Systems[END_REF] Remark 16 When jf 1 j = jf 2 j and (x) = (x) = 0, then global existence for system (1)-( 3) occurs under the conditions (46) on the reactions and the following condition on the di¤usion terms analogous to that of [START_REF] Kuiper | Positively Invariant Regions for Strongly Coupled Reaction Di¤usion Systems with a Balance Law[END_REF] but without restriction on the sign of a 12 (a 22 a 11 + a 12 a 21 ) a 21 0:

Three components systems

Let us consider the reaction-di¤usion system (1) for m=3, then the symmetrizer can be taken to be, for example

S =
Constrained by the above inequalities, we take 1 = 1: We see that i ( jf j) is a second degree polynomial with respect to i with positive leading coe¢ cient and then condition (50) is satis…ed by a large class of reactions f without conditions on their growth.

Corollary 17 Suppose that the di¤usion matrix is symmetrizable with positive and simple eigenvalues and that u 1 f 1 (u) 0 with u 1 f 2 (u) and u 1 f 3 (u) don't change sign for all ju 1 j su¢ ciently large, then global existence of solutions for system (1)-( 3) with homogenous Neumann boundary conditions and bounded initial data, occurs under conditions (50).

An eigenvalue of the di¤usion matrix is double

When an eigenvalue of A is double (for example 1 is simple and 1 < 2 = 3 is double), we should be careful when we construct the 0 i s; but the situation becomes less di¢ cult to deduce formulas analogous to (50). We de…ne 1 as above ;

but for the two characteristics corresponding to the double eigenvalue, we use Remark 1: We replace 1 by 2 in the above expression of 1 ( jf j), then by deleting the third row-column pair we obtain the …rst characteristic, denoted by 1 2 ( jf j). The second one, denoted by 2 2 ( jf j), is the second order determinant obtained from the expression of 1 ( jf j) by deleting its second row-column pair with 1 replaced by 2 . That is: The conditions (22) become 1 ( jf j) 0;

1 jf 1 j 2 jf
(a 22 2 ) 1 jf 1 j a 12 2 jf 2 j 0;

(a 33 3 ) 1 jf 1 j a 13 3 jf 3 j 0:

(51)

Since the second order determinants of the matrix A t 2 I are all null, then we can calculate explicitly the eigenvalues of A t : which shows that the coe¢ cients of the jf k j's in the inequalities (51) depend only on the di¤usion matrix, especially in the case when the balance law is strict.

Corollary 18 Suppose that the di¤usion matrix has a double eigenvalue and that u 1 f k (u) ; k = 1; 2; 3 don't change sign for all u 2 R 3 with ju 1 j su¢ ciently large, then under conditions (51) we have global existence of solutions for system (1)-( 3) with homogenous Neumann boundary conditions and bounded initial data.

Proposition 5

 5 The problem (1)-(3) admits a unique classical solution u(t; x) on an interval [0; T max [ and (i) Either ku(t; :)k 1 is bounded on [0; T max [ and the solution is global ( i.e. T max = +1). (ii) Or lim t!Tmax ku(t; :)k 1 = +1 and the solution is not global, we say that it blows up in …nite time T max or that it ceases existing.

2 j 3 jf 3 j a 12 1 2 = 1 jf 1 j 2 jf 2 j a 12 a 22 2 ; 2 2 = 1 jf 1 j 3 jf 3 j a 13 a 33 3 :

 121223

2 = 3 = a 11 a 21 a 13 a 23 and 1 =

 3111 trA 2 2 ;

  and (35) are satis…ed with 1 = 1 = 1 and 2 = 2 = 1: By application of the remark 11, the conditions (34) and (36) become respectively

	We can state the following
	Corollary 15 Suppose that conditions (44), (45) and (46) are satis…ed
	and the problem (38) has a bounded solution, then global existence of
	system (1)-(3) for bounded initial data occurs under conditions (47) when
	a 12 a 21 > 0 and (48) when a 12 a 21 < 0.
	(a 22	i ) jf 1 (t; x; u)j + a 12 jf 2 (t; x; u)j (a 22 i )		i (x) ;	i = 1; 2;
	and
	(a 22	(a 22	i )			(x) ;	i = 1; 2:
	Using (44) and (45), we have have two cases
	When a 12 a 21 > 0 then 1 < a 22 < 2 and the conditions can be
	written as follows
		(a 22	2 ) 2 (x) + (a 22		2 ) jf 1 j
		a 12 jf 2 j (a 22	1 ) jf 1 j + (a 22		1 ) 1 (x) ;
				and		(47)
		(a 22	2 ) (x) + (a 22	2 ) jg 1 j
		a 12 jg 2 j (a 22	1 ) jg 1 j + (a 22	1 ) (x) :
	When a 12 a 21 < 0 then a 11 < 1 < 2 < a 22 : The conditions (34)
	and (36) become respectively
		jf (48)

i ) jg 1 (x; u)j + a 12 jg 2 (x; u)j 1 (t; x; u)j + a 12 (a 22 i ) jf 2 (t; x; u)j i (x) ; and jg 1 (t; x; u)j + a 12 (a 22 i ) jg 2 (t; x; u)j (x) ; i = 1; 2:

  ( jf j) = 1 jf 1 j 2 jf 2 j 3 jf 3 j a 12 a 22

		1	a 32
	a 13	a 23 a 33	1

1

The author gratefully acknowledge Qassim University, represented by the Deanship of Scienti…c Research, on the material support for this research under the number (3388 ) during the academic year 1436 AH / 2015 AD.

for the reaction terms and G i j (x; z)

x i j ;1 (x) ; for all z 2 i j i j x i j ;1 ' ; and G i j (x; z) x i j ;1 (x) , for all z 2 i j i j + x i j ;1 ' ;

for the boundary conditions. We follow the same reasoning as in the proof of Theorem 6 for the homogenous case: Taking z i j = i j x i j ;1 ' with i j negative and su¢ ciently large in absolute value, then from (30) u 1 is also large in absolute value, but its sign depends on that of x i j ;1 and so we have two cases: When x i j ;1 > 0, then z i j = i j x i j ;1 ' and u 1 are both negative which gives from (21) f k = k jf k j ; k = 1; :::m. Then F i j (t; x; z) can be written on the form [START_REF] Quittner | Superlinear Parabolic Problems (Blowup, Global Existence and Steady States)[END_REF]. Thus (34) gives, from our choice of ' the …rst inequality of (42). The needed inequality on the boundary condition follows analogously. Next when x i j ;1 < 0, then z i j = i j + x i j ;1 ' and u 1 are of opposite signs with u 1 positive and again from [START_REF] Roth | Global Solutions of Reaction-Di¤usion Systems[END_REF], we have f k = k jf k j ; k = 1; :::m: Then F i j (t; x; z) can be written on the form (33). Hence from the condition (34) we get another time the …rst inequality of (42). Following the same reasoning for z 2 i j i j + x i j ;1 ' by choosing i j positive and su¢ ciently large, then u 1 is also large in absolute value, but its sign depends on that of x i j ;1 : Analogously, we have two cases: When x i j ;1 < 0, then z i j = i j x i j ;1 ' and u 1 are of opposite signs with u 1 negative which gives from [START_REF] Roth | Global Solutions of Reaction-Di¤usion Systems[END_REF], f k = k jf k j ; k = 1; :::m. Then F i j (t; x; z) has the form [START_REF] Quittner | Superlinear Parabolic Problems (Blowup, Global Existence and Steady States)[END_REF]. Consequently the second inequality of (42) is deduced from the condition (34). Finally when x i j ;1 > 0, then z i j = i j + x i j ;1 ' and u 1 is positive which gives from (21) f k = k ; k = 1; :::m. Then F i j (t; x; z) can be written on the form (33). Thus from the condition (34) we get again the second inequality of (42). We have therefore proved the …rst part of the Theorem when the di¤usion matrix is symmetric. When A is symmetrizable we follow the same reasoning as in the end of the proof for the homogenous case taking in the account that x i j ;1 should be multiplied by p s 1 .

The uniform bounds (41) of the solutions can obtained in the same way as in the homogenous case. This ends the proof of the Theorem.

Remark 11 When m = 2, we take x 0 i = (a 22 i ; a 12 ) t as …xed eigenvector associated to the simple eigenvalues i , then for any other eigen-


Caroline C Colnenne-David

Gilles Grandeau

Véronique Tanneau

M.-H Jeuffroy

Thierry Doré

Origins of the performance gaps in innovative cropping systems under experimental assessment

In the L-GHG ICS:

Very dry conditions in summer 3 years / 6 → Low amount of aerial biomass of cover crops

An unpredicted evolution of the agrosystem occurred

In both the L-EN and the L-GHG ICSs: High weed development → to mow oilseed rape plots in 2014

L-GHG

Origins of the performance gaps in the innovative cropping systems under experimental assessment

Farming System Design, Conference 2015 (September) in Montpellier (France) REFERENCES 1 Reau R. and Doré T. (2008). "Systèmes de culture innovants et durables : quelles méthodes pour les mettre au point et les évaluer ? "Educagri Editions, Dijon, France, pp 175 2 Colnenne-David C. and Doré T., 2014. Designing innovative productive cropping systems with quantified and ambitious environmental goals. "Renewable Agriculture and Food Systems". doi:10.1017/S1742170514000313.

Colnenne-David C. 1,2 , Grandeau G. 1,2 , Tanneau V. 1,2 , Jeuffroy M.H. 1,2 , Doré T.

Results: Classification of the major disparities

To meet new agricultural issues and make agriculture more sustainable, innovative cropping systems (ICSs) targeting a multiplicity of purposes need to be designed. Four ICSs were designed by prototyping (Reau and Doré, 2008 1 ) and assessed in a long-term field system experiment.

Our objective was to analyze the results after the first complete rotation, particularly the bad performances compared to the targeted aims.

Classification

Examples collected in the ICSs

Some agronomical strategies were no suitable to reach the goals

In the L-GHG ICS: No ploughing → No increase of C sequestration as expected C sequestration evolution = -149kgCO 2 ha -1 year -1 (+87kgCO 2 ha -1 year -1 expected )

Some practices were not adapted to satisfy a multiplicity of objectives

In the No-Pest ICS:

No possible to satisfy both the pesticide constraint and the S.O.M. criteria.

Restitution of small organic matter amounts + regular ploughings → Few weeds but adverse effect on C sequestration (C sequestration evolution = -560kgCO 2 ha -1 year -1 ) Some practices were not appropriate in the context of the field-trial conditions

  soil C sequestration → many cereals, continuous soil cover, high yield targets, no ploughing (ii) decrease N 2 O emissions → high number of legume crops in the rotation, improvement of N fertilization management, crops with taproots in order to reduce soil compaction 2. To satisfy multiple environmental criteria: idem PHEP ICS 3. To reach yield targets matching the Ile-de-France yields Crop sequence: catch crop (CC), maize, triticale, CC, spring faba bean, winter oilseed rape, winter wheat, CC, winter barley The L-EN ICS goals: 1. 50% fossil energy consumption compared to the PHEP ICS (i) Low direct energy consumption → no ploughing and using direct sowing machine (ii) Low indirect energy consumption → high number of legume species in the rotation, species with high N efficiency use, decrease N fertilization by reducing yield objectives 2. To satisfy multiple environmental criteria: idem PHEP ICS 3. To reach yield targets: 20% lower than the Ile-de-France yields Crop sequence: winter faba bean, winter wheat, winter flax, winter wheat-trifolium mixture, Trifolium as CC, spring oat The No-Pest ICS goals: 1. No pesticide is allowed → long rotation including a wide diversity of species (e.g. hemp), alternate sowing dates, different dates and densities of sowing, highly resistant varieties or mixtures, ploughing and mechanical weeding 2. To satisfy multiple environmental criteria: idem PHEP ICS 3. To reach yield targets: higher than organic systems in the Ile-de-France Crop sequence: triticale, CC, maize, winter wheat, CC, spring faba bean, winter wheat, CC, hemp PHEP L-GHG North Main characteristics of the field trial located in Grignon (78, France) 6.2 ha (surface plot: ≈ 4000 m²) 3 blocks deep loamy soil Beginning of field assessment: 2008 Materials and Methods  After the first complete rotation the major sources of disparities were classified  Nevertheless, a more complete agronomic diagnosis is necessary to identify and to rank all the causes of bad performances  This knowledge allowed us to improve the innovative cropping systems through a new design step  This experiment contributes to the learning design processes and cropping system management

2,1 1

  INRA, UMR 211, 78850 Thiverval-Grignon, France (caroline.colnenne@grignon.inra.fr), 2 AgroParisTech, UMR 211, 78850 Thiverval-Grignon, France

The PHEP ICS goals: 1. To satisfy multiple environmental criteria: *low pesticide use → high

	crop
	diversity, highly resistant
	varieties
	*low direct energy consumption
	→ only 1 ploughing within the
	rotation
	*low indirect energy
	consumption → legumes in
	the rotation
	*low nitrogen leaching → catch
	crop (CC) before spring crops
	and no N fertilization during
	autumn and winter
	*

stabilization and/or to enrich soil organic matter

→ burying residues of all crops