Fabrice Wendling
email: fabrice.wendling@univ-rennes1.fr

Karim Ansari-Asl

Fabrice Bartolomei

Lotfi Senhadji

From EEG signals to brain connectivity: a model-based evaluation of interdependence measures

Keywords: connectivity, regression analysis, phase synchronization, generalized synchronization, model-based evaluation, EEG

In the past, considerable effort has been devoted to the development of signal processing techniques aimed at characterizing brain connectivity from signals recorded from spatially-distributed regions during normal or pathological conditions. In this paper, three families of methods (linear and nonlinear regression, phase synchronization, and generalized synchronization) are reviewed. Their performances were evaluated according to a model-based methodology in which a priori knowledge about the underlying relationship between systems that generate output signals is available. This approach allowed us to relate the interdependence measures computed by connectivity methods to the actual values of the coupling parameter explicitly represented in various models of signal generation. Results showed that: (i) some of the methods were insensitive to the coupling parameter;

(ii) results were dependent on signal properties (broad band versus narrow band); (iii) there was no "ideal" method, i.e. none of the methods performed better than the other ones in all studied situations. Nevertheless, regression methods showed sensitivity to the coupling parameter in all tested models with average or good performances. Therefore, it is advised to first apply these "robust" methods in order to characterize brain connectivity before using more sophisticated methods that require specific assumptions about the underlying model of relationship. In all cases, it is recommended to compare the results obtained from different connectivity methods to get more reliable interpretation of measured quantities with respect to underlying coupling. In addition, timefrequency methods are also recommended when coupling in specific frequency sub-bands ("frequency locking") is likely to occur as in epilepsy.

Introduction

It is commonly admitted that most of the brain functions are based on interactions between neuronal assemblies distributed within and across different cerebral regions. For instance, it has been shown that distant areas may activate in response to a particular cognitive task, the way coordination is achieved is not resolved yet [START_REF] Uhlhaas | Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology[END_REF]. In neurological disorders like epilepsy, it has also been shown that paroxysmal activity may involve networks extending over rather large regions [START_REF] Bartolomei | Neural networks involving the medial temporal structures in temporal lobe epilepsy[END_REF].

The identification of involved networks under either normal or pathological condition is still considered as a difficult and unsolved problem. Indeed, the characterization of the functional relationship between different brain areas from a measure of the statistical coupling between signals generated by these areas is not trivial, particularly because recording techniques only allow for indirect measurement of the activity in neuronal networks.

Over the past decades, methods aimed at estimating the functional connectivity between spatially-distributed regions from electrophysiological (scalp EEG, MEG or intracerebral EEG) data have received much attention. The consequence of this increasing interest is that now a plethora of methods for estimating functional connectivity is available, all based on different assumptions about the underlying model of relationship between analyzed signals. Therefore, many neuroscientists working in this field are confronted to this crucial question: "Given a particular context of study (cognitive or clinical research), which method should be used in order to best characterize the functional interactions between distant brain areas?"

The work presented in this paper stems from this question. Indeed, considering i) the large number of methods, ii) the possible dependence of results (provided by these methods) with respect to signal properties (stationary, linearity, bandwidth) and iii) the lack of reference regarding the underlying relationship between systems that generate analyzed signals, there is a need for performance comparison.

In this paper, we report a "model-based" evaluation of a number of "well-established" methods that have been used on electrophysiological data to assess brain connectivity. In order to restrict the scope of this study, we focused on bivariate methods in locally-stationary situations (thus providing frequency-independent measures).

These methods were evaluated on the output signals of various models in which a coupling parameter can be tuned. Therefore, we were able to analyze the behavior of the methods with respect to changes of the coupling parameter and, using quantitative criteria, to compare their performances.

Background

Since the middle of the last century, numerous techniques have been proposed for measuring the temporal evolution of the cross-correlation (in a wide sense) between signals recorded from spatially-distributed brain regions.

Pioneer works started with the cross-correlation function [START_REF] Barlow | A note on a correlator for electroencephalographic work[END_REF][START_REF] Brazier | Cross-correlation and autocorrelation studies of electroencephalographic potentials[END_REF] in the time domain and with the coherence function [START_REF] Brazier | Varieties of computer analysis of electrophysiological potentials[END_REF][START_REF] Storm Van Leeuwen | The use of computer analysis for diagnosis in routine electroencephalography[END_REF]) in the frequency domain, just after fast Fourier transform (FFT) algorithms were introduced [START_REF] Cooley | An Algorithm for the Machine Calculation of Complex Fourier Series[END_REF].

These methods were applied to the study of the propagation of inter-ictal events in human intracerebral EEG data [START_REF] Brazier | Spread of seizure discharges in epilepsy: anatomical and electrophysiological considerations[END_REF]. The averaged coherence was also used on signals acquired from both hemispheres in order to study the evolution of inter-hemispheric interactions on the whole duration of partial seizures [START_REF] Gotman | Interhemispheric interactions in seizures of focal onset: data from human intracranial recordings[END_REF]. At the same time, this method revealed the existence of activities that could propagate over short-range or longrange connection fibers [START_REF] Thatcher | Cortico-cortical associations and EEG coherence: a twocompartmental model[END_REF]). Later, the averaged coherence also revealed possible synchronization mechanisms occurring at the onset of seizures [START_REF] Duckrow | Regional coherence and the transfer of ictal activity during seizure onset in the medial temporal lobe[END_REF]. In this context, a frequentlyaddressed question was also the estimation of time delays [START_REF] Gotman | Measurement of small time differences between EEG channels: method and application to epileptic seizure propagation[END_REF][START_REF] Ktonas | Estimation of time delay between EEG signals for epileptic focus localization: statistical error considerations[END_REF] as a measure of the "latency" which could provide insights into propagation phenomena between distant structures. A variant of the classical computation of the coherence function is to use of time-varying linear models (autoregressive models, multiple windowing) in order to estimate cross-and auto-spectra. Some studies reported results from these parametric methods and showed their potential value for measuring the degree of synchronization of interictal and ictal EEG signals and for characterizing the relationship between brain oscillations in the time and/or frequency domain [START_REF] Franaszczuk | An autoregressive method for the measurement of synchronization of interictal and ictal EEG signals[END_REF][START_REF] Haykin | Monitoring neural oscillation and signal transmission between cortical regions using time-frequency analysis of electroencephalographic activity[END_REF]. Later, complementary approaches were developed to estimate the direction of coupling between signals while also taking into account the possible influence of external sources (see review in [START_REF] Gourevitch | Linear and nonlinear causality between signals: methods, examples and neurophysiological applications[END_REF]).

The aforementioned methods are said to be linear. This means that the estimator that is used or the model that is assumed for signals can only capture the linear properties of the relationship between time series. However, most of the mechanisms at the origin of EEG signals are probably nonlinear. Starting from this fact, numerous studies were devoted to the development of nonlinear methods [START_REF] Pikovsky | Synchronization : a universal concept in nonlinear sciences[END_REF]. A first family based on mutual information [START_REF] Mars | Propagation of seizure activity in kindled dogs[END_REF] or on nonlinear regression [START_REF] Pijn | Propagation of electrical activity: nonlinear associations and time delays between EEG signals[END_REF]Wendling et al. 2001b) was introduced in the field of EEG about twenty years ago. A second family developed later on, based on tools already available in the field of nonlinear dynamical systems and chaos [START_REF] Iasemidis | Epileptic seizure prediction and control[END_REF][START_REF] Lehnertz | Non-linear time series analysis of intracranial EEG recordings in patients with epilepsy-an overview[END_REF]). This second family can be divided into two groups: phase synchronization (PS) methods [START_REF] Bhattacharya | Reduced degree of long-range phase synchrony in pathological human brain[END_REF][START_REF] Rosenblum | Synchronization approach to analysis of biological signals[END_REF]) and generalized synchronization (GS) methods [START_REF] Arnhold | A robust method for detecting interdependences: application to intracranially recorded EEG[END_REF][START_REF] Stam | Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets[END_REF]. PS methods estimate the instantaneous phase of each signal and then compute a quantity based on co-variation of extracted phases to determine the degree of relationship. GS methods also consist of two steps: the reconstruction of state space trajectories from time series signals and the computation of a similarity index on reconstructed trajectories.

Evaluated methods

The above brief literature review shows that the number of developed methods and variants is quite large. In this paper, we focus on 10 methods that have been applied to electro-or magneto-encephalographic signals (EEG, depth-EEG or MEG) in numerous studies. These 10 methods can be grouped into three main families: (1) regression methods: Pearson correlation coefficient (R²), coherence function (CF) and nonlinear correlation coefficient (h²); (2) phase synchronization methods: Hilbert phase entropy (HE), Hilbert mean phase coherence (HR), wavelet phase entropy (WE) and wavelet mean phase coherence (WR); (3) generalized synchronization methods: two similarity indexes (S, N) and a synchronization likelihood (SL). Main theoretical aspects are summarized hereafter.

Regression methods: R², CF and h²

For two time series ( )

x t and , the Pearson correlation R² coefficient is defined as

( ) y t ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 cov , max var var x t y t R x t y t τ τ τ + = ⋅ +
where var, cov, and τ denote respectively variance, covariance, and time shift.

The magnitude-squared coherence function (CF) can be formulated as [START_REF] Bendat | Random data: analysis and measurement procedures[END_REF]: ) f can be interpreted as the squared modulus of a frequency-dependent complex correlation coefficient. Regarding the nonlinear regression analysis, we implemented the so-called h² method. This method was introduced in the field of EEG analysis by Lopes da Silva and colleagues (Lopes da [START_REF] Lopes Da Silva | Interdependence of EEG signals: linear vs. nonlinear associations and the significance of time delays and phase shifts[END_REF]. It was also evaluated in several studies [START_REF] Ansari-Asl | Quantitative evaluation of linear and nonlinear methods characterizing interdependencies between brain signals[END_REF]Wendling et al. 2001a) using coupled oscillators. Main theoretical aspects regarding this approach were revisited in [START_REF] Kalitzin | Quantification of unidirectional nonlinear associations between multidimensional signals[END_REF]. In brief, in this method, a nonlinear correlation coefficient referred to as h² is computed based on the fitting of a nonlinear curve which approximates the statistical relationship between ( ) g ⋅ ( ) x t and ( )

( ) ( ) ( ) ( )
y t : ( ) ( ( ) ( ) (
)

) 2 var / max 1 var xy y t x t h y t τ τ τ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ + = - + where ( ) ( ) ( ) ( ) ( )) ( ( ) 2 / arg min g x t E y t t τ τ var y t g x ⎡ ⎤ + + - ⎣ ⎦
In practice, function g(.) can be obtained from the piece-wise linear approximation between the samples of the two time series ( ) x t and ( ) y t [START_REF] Pijn | Quantitative evaluation of EEG signals in epilepsy, nonlinear associations, time delays and nonlinear dynamics[END_REF].

Phase synchronization methods: HE, HR, WE, WR

The first step for estimating the phase synchronization is to extract the instantaneous phase of each signal [START_REF] Rosenblum | Synchronization approach to analysis of biological signals[END_REF]. Two different techniques are considered in this paper: the Hilbert transform and the wavelet transform. The second step is the definition of an appropriate index to measure the degree of synchronization between estimated instantaneous phases.

The Hilbert transform is used to determine the analytical signal associated to a real time series ( )

x t : ( ) ( ) [ ] ( )
( ) H x x Z t x t iH x t A t = + = ( ) H x i t e φ
where H, H

x φ , and

( )

H x
A t are respectively the Hilbert transform, the phase, and the amplitude of x(t).

The complex continuous wavelet transform can also be used for estimating the phase of a signal [START_REF] Delprat | Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies[END_REF][START_REF] Van Quyen | Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF][START_REF] Senhadji | Continuous wavelet transform: ECG recognition based on phase and modulus representations and hidden Markov models[END_REF]:

( ) ( )( ) ( ) ( ( ) * i t W x x W t x t t x t A ψ ψ ′ ′ ′ = = ⋅ ∫ ) t dt - ( ) W x t e φ =
where ψ , W x φ , and

( ) W x
A t are respectively a wavelet function (e.g., Morlet used here), the phase, and the amplitude of ( ) x t . Once the phase estimation has been performed on two considered signals, a synchronization index can be defined to quantify the phase relationship. In this paper, we present two different indexes, both based on the shape of the probability density function (pdf) of the modulo 2π phase difference ( ( ) mod 2

x y φ φ φ π

= -

).

The first index is devised from the Shannon entropy and defined by

( ) 1 1 1 ln i i p p ln M i M = = + ∑ ρ
where M is the number of bins used to obtain the pdf, is the probability of finding the phase difference i p φ within the i th bin [START_REF] Munari | Analysis of clinical symptomatology during stereo-EEG recorded mesiotemporal lobe seizures[END_REF]). The second index, refered to as the "mean phase coherence", is given by

i e φ ⎡ ⎤ ⎣ ⎦ E .
As described in [START_REF] Mormann | Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients[END_REF], it can be estimated by:

( )

1 0 N t φ - = ∑ 1 i t R e N =
where N is the length of the two time series.

Combining the two phase estimators (H, W) and the two synchronization indexes (E, R), we considered four different measures of interdependencies, denoted as Hilbert entropy (HE), Hilbert mean phase coherence (HR), wavelet entropy (WE) and wavelet mean phase coherence (WR).

Generalized synchronization based methods: S, N, SL

Generalized synchronization approaches were introduced to investigate the interactions between nonlinear dynamical systems without any knowledge about the governing equations. They generally proceed according to two steps. First, a state space trajectory is reconstructed from each scalar time series using a time delay embedding method [START_REF] Takens | Lecture Nontes in Mathematics[END_REF]. For each discrete time n a delay vector corresponding to a point in the state space reconstructed from x is defined as:

( ) ( ) 1 , , , ; 1, , n n n n m X x x x n N τ τ + + - = = … …
where m is the embedding dimension and τ denotes time lag. The state space trajectory of y is reconstructed in the same way. Second, a synchronization degree is determined using a suitable measure. Three measures, all based on conditional neighborhood, are presented in this section. The general principle is to quantify the proximity, in the second state space, of the points whose temporal indices correspond to neighbor points in the first state space. Two of these measures S and N [START_REF] Arnhold | A robust method for detecting interdependences: application to intracranially recorded EEG[END_REF], which are also sensitive to the direction of interaction, originate from this principle. They are based on an Euclidean distance:

( ) ( ) ( ) ( ) ( ) (
)

1 1 | | k N k n k n n R X S X Y N R X Y = = ∑ ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 | 1 | N k N k n n N n n R X R X Y N X Y N R X - - = - = ∑ with ( ) ( ) ( ) , 2 1 1 n j k k n n r ( ) ( ) j R X X X k = = - ∑ and ( ) , 2 1 1 | n j k k n n j R X Y X X k = = - ∑ s
where and , , 1, ,

n j r j k = … , , 1, , n j
s j= … k respectively stand for the time indices of the k nearest neighbors of n X and .

n Y

It is noteworthy that the third measure, referred to as the synchronization likelihood (SL) [START_REF] Stam | Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets[END_REF], is a measure of multivariate synchronization. Here, for simplicity, we only consider the bivariate case.

The estimated probability that embedded vectors n X are closer to each other than a distance ε is:

( ) 2 1 1 2 1 , 2( ) 1 N x n n w w j w n j w P X ε θ ε - = < -< = - ∑ j X -
where ⋅ is the Euclidean distance; θ stands for the Heaviside step function, is the Theiler correction and determines the length of sliding window. Letting

ε from which we can determine if simultaneously n X is close to j X and is close to n Y j Y , n j H , i.e., 2 = in the equation below ( ) ( ) , , , n j x n n j y n n j H X X Y Y θ ε θ ε = - - + --
Synchronization likelihood at time n can be obtained by averaging over all values of j ( ) ( )

1 2 , 1 2 1 1 1 2 N n j ref w n j w SL H P w w = < -< = - ∑ n j -
All aforementioned measures are normalized between 0 and 1; the 0 value means that the two signals are completely independent. On the opposite, the 1 value means that the two signals are completely synchronized.

Finally, in order to deal with the evolution, in time, of brain connectivity the three measures described above can be estimated over a sliding window.

Model-based evaluation methodology

In order to perform an objective comparison of methods, we propose a model-based approach in which a priori knowledge about the underlying relationship between systems that generate output signals is available. As illustrated in figure 1, connectivity methods are applied on the output signals X and Y produced by two models of signal generation S 1 and S 2 . The coupling between these two models can be adjusted using a coupling parameter C. Some noise (N 1 and N 2 ) can be added to output signals. Evaluated methods provide a quantity Q (normalized between 0 and 1) measured from X and Y and supposed to characterize the connectivity between the two models.

Therefore, this framework offers the possibility to study the behavior of candidate connectivity methods with respect i) to the coupling parameter C, ii) to the type of model that is used to generate X and Y and iii) to the features of noise present on output signals. In this context, the "ideal" method would provide a quantity Q which would reflect any variation of coupling parameter C for any model of signal generation, with low bias and variance and high sensitivity. This consideration led us to propose different models for generating outputs signals (section 4.1) and three criteria (section 4.2) to evaluate the behavior of connectivity methods with respect to the variation of parameter C.

Models of signal generation

Table 1 gives an overview of signal generation models used in this study. As depicted, these models can be divided into three main families: i) models of coupled stochastic signals (M1 and M2), ii) models of coupled nonlinear dynamical systems (M3 and M4) and iii) models of coupled neuronal populations (M5). In each family, various situations were considered in order to check for the influence of signal properties (narrow versus broad band activity) or the influence of noise on connectivity measures. The main characteristics of each model are described hereafter.

Coupled broad band signals: M1

Model M 1 generates two broad band signals ( 1 , 2 x x ) built using two independent and a common white noises respectively (N 1 , N 2 ) and (N3):

( ) ( )

1 1 2 2 1 1 x C N C N 3 3 x C N CN = - + = - +
where is the coupling degree; for 0 C ≤ ≤1 0 C = the signals are independent and for they are identical.

1 C =

Coupled narrow band signals: M2(PR) and M2(AR)

In M 2 , two narrowband signals around a frequency 0 f are generated from four lowpass filtered white noises (NF 1 , NF 2 , NF 3 , and NF 4 ). They are combined in two ways in order to share either a phase relationship (PR) or an amplitude relationship (AR), only: We also used some nonlinear deterministic systems for modeling nonlinear synchronized coupled oscillators.

( ) ( ) ( ) ( ) ( ) ( ) ( 1 1 0 1 2 2 0 1 2 1 1 0 1 2 1 2 0 cos 2 : cos 2 1 cos 2 : 1 c o s 2 x A f t PR x A f t C C x A f t AR x CA C A f t π φ π φ φ π φ ) 2 π φ = + = + + - = + = + - + ⎧ ⎨ ⎩ ⎧ ⎨ ⎩
Here we report results obtained for two of them: Rössler-Rössler coupled system (M 3 ) [START_REF] Pikovsky | Synchronization in a population of globally coupled chaotic oscillators[END_REF] and Hénon-Hénon coupled system (M 4 ) [START_REF] Bhattacharya | Reduced degree of long-range phase synchrony in pathological human brain[END_REF])

Coupled Rössler systems: M3

In M 3 two Rössler systems [START_REF] Rossler | An equation for continuous chaos[END_REF]) are coupled and the driver system is given by: ( ) Hénon map [START_REF] Hénon | A two-dimensional mapping with a strange attractor[END_REF]) is a nonlinear deterministic which is discrete by construction. In model M 4 , we make use of two of them to simulate a unidirectional coupled system. The driver system is given by:

[ ] [ ] [ ] 2 1 1.4 1 x x n x n b x n + = - + - and the response system is [ ] [ ] [ ] ( ) [ ] ( ) [ ]

Comparison criteria

Three criteria were defined in order to quantitatively evaluate the behavior of the methods presented in section 3 on time-series simulated from the coupled models presented in section 4.1 for different values of the coupling parameter C. The two first criteria are classical:

-the mean square error (MSE) under null hypothesis (i.e., independence between two signals) can be interpreted as a quadratic bias, defined by where E is the mathematical expectation, (no coupling) and is the estimation of Q , ( ) (

2 Q Q - ⎡ ⎣ E ⎤ ⎦ I 0 Q = Q -the
2 1 1 Î i i i Q Q I = ⎡ - ⎢ ⎣ ⎦ ∑ E E ) ⎤ ⎥
where I is number of coupling degree points and is the estimated relationship for the coupling degree .

ˆi

Q i c
The third performance criterion, referred to as the median of local relative sensitivity (MLRS), is novel. It was introduced to quantify the sensitivity of the method with respect to changes in the coupling degree. The MLRS is given by: ( )

2 2 1 1 1 ˆˆˆ, , 2
i i i i i i i i i i Q Q MLRS Median S S c c σ σ σ σ + + + - + = = = -
where is the increase rate of the estimated relationship and i S i σ is the square root of the average of estimated variances associated to two adjacent values of the coupling degree. One can also use the median of the distribution of local relative sensitivity instead of its mean because the fluctuation in its estimation may make this distribution very skewed.

Methods with lower values of MSE and MV can be considered to have better performances. Conversely, higher MLRS values indicate better performances.

Results

For all models of signal generation (M1-M5) and for all values of parameter C (from 0 to 1 by steps of 0.1), Monte Carlo simulations were performed in order to assess statistical properties of connectivity measures provided by methods described in section 3 and to comparatively evaluate their performances according to criteria presented in section 4.2. To proceed, long duration time-series signals (200000 samples) sampled at 256 Hz were generated using each model, for discrete values of the coupling parameter C, ranging from 0 to 1, by steps of 0.1. All quantities were estimated over a sliding window of 512 samples (2 seconds of activity) moving by steps of 64 samples (75% overlapping) and then averaged over time. In GS methods, parameter τ was estimated as follows. First the mutual information as a function of positive time lag was plotted and then, as described in [START_REF] Sills | Concentration-effect studies with topiramate on selected enzymes and intermediates of the GABA shunt[END_REF], time lag τ was chosen as the abscissa value corresponding to the first minimum this curve.

The embedded dimension m, in this family of methods, was determined from the Cao method (Cao 1997).

To limit the length of this section, we will only detail the results obtained in two cases (models M2(PR) and M5

(SPK)) and provide a global synthesis that summarizes the performances of the three families of methods in all considered situations.

Results obtained in model M 2 for phase relationship only (PR), are shown in figure 3. They show that phase synchronization methods (HE, HR, WE, WR) exhibited higher performances than other methods, as it could be expected. Similarly, results obtained with R² and h² methods were correct. On the opposite, generalized synchronization methods and coherence (CF) had lower performances (S, SL: low sensitivity to C variation, S: strong bias, SL: strong variance, CF: low sensitivity). Regarding the variance of estimation, it can be noticed that some methods (R², WE) perform better than others. For higher coupling values, the variance decreases for all methods, except SL. As depicted, most of the methods could detect the increase of the coupling parameter in the model. Some methods exhibited strong bias (HE). As an interesting result, it could be observed that WE and CF were almost blind to the increase of coupling. Similarly, HE and WR only displayed small increase with increasing of degree of coupling but their variance was low. R², h², S and HR methods exhibited good sensitivity. However, MSE under null hypothesis was found to be high for HR. Globally, the variance was found to decrease with the coupling parameter all methods, except CF.

The two examples described above show that the behavior can strongly vary from one method to another within the same situation. In order to make a global comparison of the 10 studied methods, we computed the MSE under null hypothesis, the MV and the MRLS (both averaged over the coupling values) for all situations. Results are summarized in table 2. From this table, it can be observed that the R² method exhibited good results in model M1 and M5(BKG) and that the h² method showed good sensitivity to coupling changes in model M5 (SPK).

Among regression methods, the average coherence CF showed poorer performances. Among phase synchronization methods, the wavelet entropy (WE) performed better than other methods of the same family, in particular in model M2(PR). It is noteworthy that the four phase synchronization methods were not able to detect coupling changes in model M2(AR). Regarding generalized synchronization methods, results showed that best performances are obtained in models of coupled nonlinear systems. In particular, the best results were obtained with the similarity index S in model M4b.

We also averaged the results obtained in each family of methods (regression, PS or GS). Results are synthesized in figure 5. On one hand, a criterion-by-criterion analysis leads to the following remarks. Regarding the sensitivity to the coupling parameter (MRLS, figure 5-a for the three families. In models M2 (AR), lower variance was obtained for GS methods still keeping some sen-sitivity to the coupling parameter. PS methods showed very low variance in model M3 (Rössler). In model M4 (Hénon), GS methods exhibited the highest variance but were, in turn, the most sensitive to the coupling parameter. For model M5, PS and regression methods performed better than GS methods. Regarding the mean square error under null hypothesis (MSE, figure 5-c), results showed that some methods can be characterized by a high bias value like GS methods in model M2, PS methods in model M3 and GS methods in model M5.

On the other hand, a model-by-model analysis of figure 5 shows that for model M 1 regression methods perform better than others as the MV is the lowest while the MRLS is the highest. For model M 2 (in the case of PR), it is evident that PS methods are the most appropriate. For model M 2 (in the case of AR), there is no consensus for the best method. For model M 3 , PS methods outperform others although they are characterized by higher MSE values. For model M 4 (considering the four situations), GS methods have the lowest MSE and PS methods have the lowest MV. As far as the MRLS is concerned, these two groups of methods perform equally. Finally, for the neuronal population model M 5 , regression methods outperform others in the case of normal background EEG activity. For spiking epileptic-like activity, these methods, in addition to PS methods have also higher performances than GS methods.

Discussion

Relevance in the field of brain research and progress with respect to already published studies

Numerous methods have been introduced to tackle the difficult problem of characterizing the statistical relationship between signals without a priori knowledge about the nature of this relationship. In brain research, a number of methods have been proposed and/or used to study normal or pathological processes. Such methods play a key role as they are supposed to provide important information regarding brain connectivity from electrophysiological recordings.

During the ten past years, some efforts have been made for comparing methods but mainly qualitatively [START_REF] David | Evaluation of different measures of functional connectivity using a neural mass model[END_REF][START_REF] Quian Quiroga | Performance of different synchronization measures in real data: a case study on electroencephalographic signals[END_REF] and for particular applications [START_REF] Mormann | On the predictability of epileptic seizures[END_REF][START_REF] Pereda | Comparing Generalized and Phase Synchronization in Cardiovascular and Cardiorespiratory Signals[END_REF].

In this paper, we presented three families of methods (regression, phase synchrony and general synchronization)

and we compared their performances on the basis of simulations produced by various models (including a neurophysiologically-relevant one) in which a parameter can be tuned in order to adjust the coupling between two "systems" from which output signals are generated. Through this model, we were able to study the relationship between the coupling parameter in the models and the quantity actually measured on output signals using "connectivity" methods.

In this regard, this approach differs from that of [START_REF] Schiff | Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble[END_REF] who evaluated one method to characterize dynamical interdependence (based on mutual nonlinear prediction) on both simulated (coupled identical and non identical chaotic systems) and real (activity of motoneurons within a spinal cord motoneuron pool) data. It also differs from other evaluation studies which mainly focused on qualitative comparisons [START_REF] David | Evaluation of different measures of functional connectivity using a neural mass model[END_REF][START_REF] Quian Quiroga | Performance of different synchronization measures in real data: a case study on electroencephalographic signals[END_REF] and for specific applications [START_REF] Mormann | On the predictability of epileptic seizures[END_REF][START_REF] Pereda | Comparing Generalized and Phase Synchronization in Cardiovascular and Cardiorespiratory Signals[END_REF]). In the particular field of EEG analysis, the model of coupled neuronal populations is of particular relevance since it generates realistic temporal dynamics. In this model, for background activity (that can be considered as a broadband random signal), we found that coherence and phase synchrony methods (except HR) were not sensitive to the increase of the coupling parameter whereas regression methods (linear and nonlinear) exhibited better sensitivity. This result may be explained by the fact that the interdependence between simulated signals is not entirely determined by a phase relationship. This point is crucial since it illustrates the fact that the choice of the method used to characterize the relationship between signals is critical and may lead to possible incorrect interpretation of results obtained on EEG data.

In addition, as background activity can be recorded in epileptic patients during interictal periods, our results also relate to those recently published by Morman et al. [START_REF] Mormann | On the predictability of epileptic seizures[END_REF] in the context of seizure prediction.

For thirty different measures obtained from univariate and bivariate approaches, authors evaluated their ability to distinguish between the interictal period and the pre-seizure period (sensitivity and specificity of all measures were compared using receiver-operating-characteristics). In both types of approach (and consequently for bivariate methods similar to those implemented in the present study) they also found that linear methods performed equally good or even better than nonlinear methods.

Limitations of this study

In this report, results about the characterization of the direction of coupling were not dealt with. This difficult issue has already been addressed in various reports. For instance, Quian Quiroga et al. [START_REF] Quian Quiroga | Performance of different synchronization measures in real data: a case study on electroencephalographic signals[END_REF] quantitatively tested two interdependence measures on coupled nonlinear oscillators for their ability to determine whether one the two systems drives the other. Other families of approaches where also developed to estimate directional properties of the relationship between signals while taking into account possible influences of external sources [START_REF] Gourevitch | Linear and nonlinear causality between signals: methods, examples and neurophysiological applications[END_REF]). According to these approaches, causality between two signals (in Granger's sense) is estimated based on the predictability of one signal from the immediate past of the other signal.

Besides, the methods evaluated in this study are independent from frequency. However, in some situation, it might be crucial to take frequency into account. Indeed, frequency-independent methods may not be able to reveal some phenomena like a hypersynchronization in a narrow frequency band as sometimes observed at the onset of partial seizures. An example is illustrated in figure 6 More importantly, it should be noticed that time-independent methods (R² and h²) that showed relatively robust in simulations do not show significant changes at seizure onset (figure 6-d and 6-e). As these methods apply on broadband signals, they are unable to detect the increase of correlation occurring at a specific frequency.

Summary of the main findings

The main findings of this study can be summarized as follows: (i) some of the compared methods are insensitive to the coupling parameter in the model; (ii) results are dependent on signal properties (broad band versus narrow band); (iii) generally speaking, there is no universal method to deal with signal coupling, i.e., none of the studied methods performed better than the other ones in all studied situations.

Nevertheless, we notice that simple methods like R² and h² methods showed to be sensitive to the coupling parameter in the model with average or good performances. Therefore, it might be reasonable to first apply these "robust" regression methods in order to characterize brain connectivity before using more sophisticated methods that require specific assumptions about the underlying model of relationship. In all cases, it is highly recommended to compare the results obtained from different connectivity methods to get more reliable interpretation of measured quantities with respect to underlying coupling. Following this idea, results suggest that some information could be inferred about the nature of the relationship based on the comparison of indexes provided by the different methods although the question of the physiological plausibility of the nature of the relationship will always remain.

In addition, some coupling phenomena might correspond to increased correlation between output signals, but in one or several specific sub-band(s). In such cases, time-frequency methods may also be of help. In this category, it should be reminded that a characterization of the relationship based on a "same frequency to same frequency" approach can only reveal the linear component of this relationship. A more general approach would allow for the computation of the correlation (in a wide sense) between signals filtered in different frequency sub-bands but might also lead to additional difficulties in the interpretation of results. In some cases, methods were found to be insensitive to the variations of the coupling parameter C in the models. In these cases, the three criteria are not applicable (N/A). Averaged MRLS Averaged MV

Captions

Averaged MSE

Figure 5: A summary of the performances (according to quantitative criteria introduced in this study) of the three families of methods for all considered models of signal generation. For M2(AR), results obtained for phase synchronization methods are not represented since these methods were found to be insensitive w.r.t changes of the coupling parameter in this case (see table 2).

  t and ( ) y t , and where denotes their cross-spectral density. It is the counterpart of the R² coefficient in the frequency domain. It ( xy S

  two signals have independent phase and amplitude and for 1 C = they have identical phase or amplitude.

  Hénon systems: M4a and M4b

2 N

 2 to 0.3 to have two identical systems (M 4a ) and once y b is set to 0.1 to have two non-identical systems (M 4b ).4.1.5 Coupled neuronal populations: M5In order to simulate realistic temporal dynamics encountered in real depth-EEG signals, as recorded in patients with drug-resistant epilepsy during pre-surgical evaluation, a physiologically-relevant computational model of EEG generation was used. As illustrated in figure2, this model represents the activity of two distant -and possibly coupled -neuronal populations. Each population generates a local field potential (X and Y) that can be seen as a depth-EEG signal if one does not consider the source-electrode transfer function. Readers may refer to[START_REF] Wendling | Relevance of nonlinear lumpedparameter models in the analysis of depth-EEG epileptic signals[END_REF] for detailed description and generalization to more than 2 populations. In this model, each population contains two subpopulations of neurons that mutually interact via excitatory or inhibitory feedback: main pyramidal cells and local interneurons. The excitatory influence from non-specific afferents is modeled by a Gaussian input noise ( or ) that globally represents the average density of afferent action potentials on population 1 and population 2. Pyramidal cells are excitatory neurons that project their axons to distant areas of 1 N the brain. The model accounts for this type of connection by using the average pulse density of action potentials from the main cells of population 1 as an excitatory input to population 2. In addition, this uni-directional connection from population 1 to population 2 is characterized by parameter C which represents the "synaptic gain", allowing for tuning the degree of coupling between the two populations. Other parameters include excitatory and inhibitory gains in feedback loops as well as average number of synaptic contacts between subpopulations. The model was used to generate two kinds of signals: background (BKG) and spiking (SPK) activity respectively obtained for normal and increased excitation/inhibition ratio in each population. For both cases, bivariate data were simulated. The normalized coupling parameter C was varied from 0 (independence between signals) to 1 (dependence between signal, similar temporal dynamics).

  mean variance (MV) computed over all values C = ,

Figure 4

 4 Figure4shows the signals generated by the neuronal population model (M5) in the case of spiking (SPK) activity. The temporal dynamics of these signals reproduce those often observed during epileptic seizures. It is noteworthy that similar signals were used in a previous attempt for qualitatively comparing different correlation estimators[START_REF] Quian Quiroga | Performance of different synchronization measures in real data: a case study on electroencephalographic signals[END_REF]. As depicted, most of the methods could detect the increase of the coupling

  ), regression methods (R², h² and CF) and PS methods (HE, HR, WE, WR) perform better than GS (S, N, SL) methods in model M1. As it could be expected, PS methods showed no sensitivity to coupling in model M2(AR) and highest sensitivity in model M2(PR). Globally GS methods showed good sensitivity to coupling in model M4(Hénon) but were outperformed by PS methods in model M3 (Rössler). Regression methods (R², h² but not CF), showed good sensitivity to the coupling parameter in the model of coupled neuronal populations. Regarding the mean variance (MV, figure5-b), regression methods (R², h²) displayed low values in model M1. In models M2 (PR), quasi-similar variance values were obtained

  which displays the depth-EEG signals recorded from mesial temporal lobe structures in a patient candidate to surgery. As observed, the seizure onset is marked by the appearance of a fast activity clearly localized in the time-frequency representations (figure 6-b) of the two depth-EEG signals (figure 6-a). A "frequency-locking" phenomenon is revealed by the computation of the linear correlation as the function of time and frequency (see methods in (Ansari-Asl et al. 2005)). As shown in figure 6-c, at seizure onset, the two depth-EEG signals are correlated in a well-defined frequency sub-band (around 30 Hz).

Figure 1 :

 1 Figure 1: general model proposed to study the behavior of bivariate methods aimed at characterizing the connectivity between coupled systems (unidirectional coupling is tuned using parameter C) from output signals X and Y.

Figure 2 :

 2 Figure 2: Model of coupled neuronal populations used to generate bivariate time-series data that mimic the temporal dynamics of actual depth-EEG signals (local field potentials).

Figure 3 :

 3 Figure 3: Results obtained from model M2 in the case coupling parameter C induces a phase relationship only between signals. (a) Example of simulated signals. (b) Estimated relationship (mean value of the estimated quantity Q(X,Y) over all realizations as a function of the coupling degree in the model). (c) Variance of estimation.

Figure 4 :Figure 5 :

 45 Figure 4: Results obtained from the neuronal population model in the case of epileptic EEG activity (sustained spiking activity as observed during ictal periods in epilepsy). (a) Example of simulated signals. (b) Estimated relationship (mean value of the estimated quantity Q(X,Y) over all realizations as a function of the coupling degree in the model). (c) Variance of estimationFigure5: A summary of the performances (according to quantitative criteria introduced in this study) of the three families of methods for all considered models of signal generation. For M2(AR), results obtained for phase syn-chronization methods are not represented since these methods were found to be insensitive w.r.t changes of the cou-pling parameter in this case (see table2).

Figure 6 :

 6 Figure6: A "frequency-locking" phenomenon observed at the onset of seizure in a patient with medial temporal lobe epilepsy. a) depth-EEG recordings from medial structures (X(t) : hippocampus, Y(t) : entorhinal cortex) performed in a patient with mesial temporal lobe epilepsy. b) Time-frequency representation of the signal energy (spectrogram method computed using the short-term Fourier Transform -STFT-of signals X and Y). White arrows show the occurrence of a fast discharge (around 30 Hz) at seizure onset. c) Time-frequency representation of the linear relationship between signals X and Y. The method consists in the computation of the linear correlation coefficient between signals filtered in narrow sub-bands. Grey arrow shows that the linear correlation increases within the specific frequency sub-band corresponding to the fast activity. d, e) Results obtained from "standard" (i.e. frequency-independent) linear (R²(t)) and nonlinear (h²(t)) methods. Black arrows show that the "localized-in-frequency" correlation increase is not detected.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: general model proposed to study the behavior of bivariate methods aimed at characterizing the connectivity between coupled systems (unidirectional coupling is tuned using parameter C) from output signals X and Y.

Table 1 :

 1 Models of signal generation. Abbreviation SNR denotes "Signal-to-Noise Ratio".

Table 2 :

 2 Summary of results. For each model (horizontal: M1-M5) and for each studied method (vertical: R²-SL), the value of the mean square error (MSE) under null hypothesis (H0, i.e. no coupling), the mean variance (MV) and the median of local relative sensitivity (MLRS) were computed from simulated signals. Values highlighted in grey color correspond to the "best performances', i.e. lowest MSE, lowest MV and highest MLRS.

Table 1 :

 1 Models of signal generation. Abbreviation SNR denotes "Signal-to-Noise Ratio".

Table 2 :

 2 see Captions section for legend

					Coupled stochastic			Coupled nonlinear		Coupled neuronal
						Signals				dynamical systems		populations
					M 1	M 2 (PR)	M 2 (AR)	M 3	M 4a (SNR=inf)	M 4a (SNR=2)	M 4b (SNR=inf)	M 4b (SNR=2)	M 5 (SPK)	M 5 (BKG)
				MSE (H 0 )	0.12	76.42	0.2	109.55	0.28	0.22	0.26	0.22	63.17	1.54
	Regression		R² CF	MV MRLS MSE (H 0 ) MV MRLS	3.6 57.6 104.48 5.0 56.4	200.1 3.94 108.22 203.7 1.3	366.6 0.41 N/A N/A N/A	65.5 1.38 91.14 199.4 2.20	51.7 6.99 107.14 26.6 5.68	23.9 2.93 107.83 13.1 1.94	17.4 21.31 102.53 11.6 17.17	10.9 16.85 104.17 9.2 9.41	215.8 1.3 N/A N/A N/A	21.2 1.2 N/A N/A N/A
				MSE (H 0 )	1.10	117.45	0.31	151.14	3.43	1.33	2.73	1.12	103.79	5.99
			h²	MV	3.7	161.10	274.8	60.4	42.0	23.3	21.5	11.8	205.0	22.6
				MRLS	35.6	4.06	0.36	0.98	6.91	0.56	20.62	16.42	1.0	1.1
				MSE (H 0 )	4.98	23.11	N/A	63.82	8.09	6.87	6.00	5.69	28.75	N/A
			HE	MV	2.4	57.2	N/A	19.0	25.5	4.2	4.2	2.3	45.3	N/A
				MRLS	40.9	6.58	N/A	15.5	7.16	3.98	20.54	15.58	1.2	N/A
	Phase	Synchronization	HR WE	MSE (H 0 ) MV MRLS MSE (H 0 ) MV MRLS	3.00 8.8 42.5 10.54 1.5 47.0	175.56 206.7 6.5 20.48 29.4 6.69	N/A N/A N/A N/A N/A N/A	473.07 10.1 8.87 76.47 6.3 13.8	6.26 49.6 6.79 13.01 20.1 11.8	4.99 28.0 4.77 13.27 3.3 5.07	2.74 17.9 21.45 12.98 3.1 19.20	3.04 14.8 15.38 12.76 1.7 11.75	249.31 217.5 1.2 N/A N/A N/A	18.99 65.5 0.7 N/A N/A N/A
				MSE (H 0 )	8.78	113.65	N/A	161.50	68.97	65.79	62.26	64.22	53.39	N/A
			WR	MV	6.0	118.6	N/A	17.6	23.2	14.6	11.4	10.1	38.4	N/A
				MRLS	46.6	6.76	N/A	8.83	10.27	4.68	20.02	14.80	1.2	N/A
	Generalized Syn-	chronization	S N	MSE (H 0 ) MV MRLS MSE (H 0 ) MV MRLS MSE (H 0 )	75.51 2.8 31.1 0.41 5.0 29.0 4.28	28.45 58.1 2.23 378.44 120.5 3.02 115.25	0.17 74.6 0.84 0.6 142.3 0.60 0.31	26.59 48.8 6.91 116.76 68.4 3.46 8.50	0.03 20.8 19.40 0.63 12.5 25.42 4.18	27.47 6.6 10.26 0.30 1.36 15.70 4.47	0.03 1.2 31.51 0.55 8.4 12.32 3.83	27.41 4.7 18.03 0.33 13.4 17.04 3.74	107.53 183.7 0.9 201.46 501.1 0.4 41.32	120.04 44.1 0.05 4.90 60.1 0.9 6.16
			SL	MV	44.8	253.4	209.1	253.5	104.2	138.2	163.5	139.2	383.8	59.2
				MRLS	8.32	0.772	0.41	3.52	2.41	2.61	4.90	3.84	1.3	0.007

"standard" (i.e. frequency-independent) linear (R²(t)) and nonlinear (h²(t)) methods. Black arrows show that the "localized-in-frequency" correlation increase is not detected.


Nader El Khatib

Mamdouh Nicolas Forcadel

M Zaydan

Regularizing effect for unbounded flux-limited viscosity solutions of a discontinuous Hamilton-Jacobi equation on junction

Keywords: AMS Classification: 35F21, 35D35, 35D40 Regularizing effect, viscosity solutions, discontinuous Hamilton-Jacobi on a junction

Regularizing effect for unbounded flux-limited viscosity solutions of a discontinuous Hamilton-Jacobi equation on junction

Introduction

The term regularizing effect indicates that the solution of a non-linear PDE becomes more regular for t > 0 than it is at the initial time t = 0. Such effects were studied in the framework of first order continuous Hamilton-Jacobi equations in many papers like [START_REF] Barles | Regularity results for first order hamilton-jacobi equations[END_REF][START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF][START_REF] Lions | Regularizing effects for first-order hamilton-jacobi equations[END_REF][START_REF]New regularity results for hamilton-jacobi equations and long time behavior of pathwise (stochastic) viscosity solutions[END_REF]. In these works, the regularizing effect is obtained for Hamiltonian H satisfying that (H p • p -H)(p) is large when |p| is large or when H(p) is large. This regularizing effect can be used to study the large time behavior of viscosity solutions. For example, in [START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF] the authors proved that the viscosity solution u satisfy u t ≥ -η(t)e u . This result implies first more regularity in space: local Lipschitz continuity if the Hamiltonian is coercive and Holder regularity for hypo-elliptic Hamiltonian. As a second application, this result can be used in the study of the large time behavior of the viscosity solution if η(+∞) = 0. For more details, the reader can refer to [START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF][START_REF] Barles | On the large time behavior of solutions of hamiltonjacobi equations[END_REF][START_REF]An introduction to the theory of viscosity solutions for first-order hamilton-jacobi equations and applications[END_REF]. Let us mention also the work [START_REF]New regularity results for hamilton-jacobi equations and long time behavior of pathwise (stochastic) viscosity solutions[END_REF] where authors exploited the regularizing effect for a Hamilton-Jacobi equation in order to study the large time behavior of the stochastic viscosity solution.

In this paper, we want to develop regularizing results for discontinuous Hamilton-Jacobi equations. The study of discontinuous (in space) Hamilton-Jacobi equations is recent (see [START_REF] Achdou | Hamilton-jacobi equations constrained on networks[END_REF][START_REF] Imbert | A hamilton-jacobi approach to junction problems and application to traffic flows[END_REF]) but the basic results, like comparison principle, existence, stability, homogenization are now well understood (see [START_REF] Imbert | Flux-limited solutions for quasi-convex hamilton-jacobi equations on networks[END_REF][START_REF] Lions | Well-posedness for multi-dimensional junction problems with kirchoff-type conditions[END_REF][START_REF] Barles | Flux-limited and classical viscosity solutions for regional control problems[END_REF][START_REF] Imbert | Quasi-convex hamilton-jacobi equations posed on junctions: the multi-dimensional case[END_REF][START_REF] Siconolfi | Time-dependent hamilton-jacobi equations on networks[END_REF][START_REF] Achdou | Effective transmission conditions for hamiltonjacobi equations defined on two domains separated by an oscillatory interface[END_REF][START_REF]Effective transmission conditions for second-order elliptic equations on networks in the limit of thin domains[END_REF][START_REF] Khatib | Homogenization of a microscopic pedestrians model on a convergent junction[END_REF]) and we are now able to attack more difficult problems like the regularity or the regularizing effect.

More precisely, in this paper, we consider regularizing effect for a Hamilton-Jacobi posed on a junction, which is a metric space formed by the union of a N branches J 1 , J 2 , .., J N glued to one point called junction point. For i ∈ {1, ..., N }, the branch J i is isometric to R + and the junction point is x = 0. We denote by J the junction and it is defined by J = i=1,...,N J i with J i ̸ = J j if i ̸ = j.

To be more precise, the definition of elements in J is given by the following definition: if x ∈ J, then x = x i • e i , x i > 0 if x ∈ J * i = J i \{0} or x = (0, 0) where e i is a unit vector in R 2 and e i ̸ = e j if i ̸ = j. We define on the junction J the geodesic distance by

d(x, y) = |x i -y i | if x, y ∈ J * i , |x i | + |y j | if x ∈ J i and y ∈ J j with i ̸ = j.
For T > 0, let u : [0, T ]×J → R be a real function. The space gradient of u at (t, x) ∈ [0, T ]×J is defined by

u x (t, x) = ∂ i u(t, x) if x ∈ J * i , (∂ 1 u(t, 0), ..., ∂ N u(t, 0)) if x = 0
where ∂ i u(t, x) is the derivative of u with respect to x ∈ J * i . We consider a continuous bounded by below viscosity solution of the following Hamilton-Jacobi with flux limited condition at the point 0 introduced in [START_REF] Imbert | Flux-limited solutions for quasi-convex hamilton-jacobi equations on networks[END_REF] and given by

u t + H i (t, x, u x ) = 0 if (t, x) ∈ [0, T ] × J * i , u t + F A (t, 0, u x ) = 0 if t ∈ [0, T ], (1.1)
where for each i ∈ {1, ..., N }, the Hamiltonian H i satisfies the following assumptions:

     H i ∈ C([0, T ] × J i × R), for all (t, x) ∈ [0, T ] × J i , the Hamiltonian H i (t, x, •) is quasi-convex, for all (t, x) ∈ [0, T ] × J i , the Hamiltonian H i (t, x, •) is coercive.
For t ∈ [0, T ] and p = (p 1 , ..., p N ) ∈ R N , the flux limiter function F A is defined by

F A (t, 0, p) = max A(t), max i=1,...,N H - i (t, 0, p i ) (1.2)
where

H - i (t, 0, •) is the non-increasing part of H i (t, 0, •).
The function A is called the "flux limiter" and we assume that A ∈ C([0, T ]). Moreover, we suppose that for t ∈ [0, T ],

A(t) ≥ max i∈{1,...,N } min p∈R H i (t, 0, p).
Additional assumptions (H) (see below) are required to obtain a regularizing effect. These assumptions are satisfied for example by the following Hamiltonians,

H i (t, x, p) = |A i (t, x)| mi |p| mi -f i (t, x)
where m i > 1, f i is non negative (or bounded by below), A i and f i are space-time Lipschitz continuous functions and there exists 0

< δ ≤ 1 such that δ ≤ |A i (t, x)| ≤ 1 δ for (t, x) ∈ [0, T ] × J i .
More precisely, we will provide a lower bound for the time derivative of the viscosity solution u.

Up to constant addition, we can assume that u ≥ 0. We will prove that there exists a time t * > 0 and function η defined on (0, t * ) such that for all t ∈ (0, t * ) and x ∈ J, we have ,x) .

u t ≥ -η(t)e u(t
(1.3)

Using inequality (1.3), we can prove that u is locally Lipschtiz continuous in J × (0, t * ) (see Proposition 2.2 ). To obtain (1.3), we will first provide a lower gradient bound on each branch of the junction using the assumption ((H i ) p • p -H i )(t, x, p) is positive and large when p is negative and small (Lemma 4.3). This bound, joint to the idea of studying the sign of the derivative of d 2 (x, y) used in the proof of the comparison principle in [START_REF] Barles | Flux-limited and classical viscosity solutions for regional control problems[END_REF] will allow us to obtain our main result by exploiting the assumption ((H i ) p •p-H i )(t, x, p) is positive and large when H i (t, x, p) is positive and large.

Organization of the paper

In section 2, we state the assumptions imposed on the Hamiltonians and then we give the main result (Theorem 2.1) of this paper. In section 3, we derive consequences from assumptions of section 2 and then we give examples of Hamiltonians. In section 4, first, using the variable change v = -e -u , we obtain a new Hamilton-Jacobi formulation and we derive properties of the new Hamiltonians G i . Then, we prove the main result of this paper using this new equation (Theorem 4.2). In the last section, as application of Theorem 2.1, we show that our solution is locally Lipschitz.

Assumptions on the Hamiltonians and main results

Assumptions (H)

Our main goal is to process the proof in the vicinity of the junction point x = 0. Far from it, we work with one type of Hamiltonian and we can impose weaker assumptions (like non-coercivity) on the Hamiltonians. To be more precise, far from the junction point, we can use the assumptions on Hamiltonians provided in [START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF]. Therefore, all the assumptions we will see below can be imposed on a neighborhood of the junction point. For simplicity, we assume that they are satisfied for all x ∈ J. Let u : [0, T ] × J → R be a continuous, non-negative viscosity solution of (1.1). We set

C u = min t∈[0,T ]
(e -u(t,0) ) > 0.

(2.1)

We state now the assumptions concerning the Hamiltonians.

(H0) For all i ∈ {1, ..., N } and for all t ∈ [0, T ], x ∈ J i , the Hamiltonian H i (t, x, •) is quasi-convex and coercive on R and

H i is continuous in [0, T ] × J i × R.
1) We assume that there exists an interval [p - i (t, x), p + i (t, x)] such that

     H i (t, x, •) is non increasing on (-∞, p - i (t, x)), H i (t, x, •) is constant on [p - i (t, x), p + i (t, x)], H i (t, x, •) is non decreasing on (p + i (t, x), +∞).
Moreover, we assume that p - i is bounded in [0, T ]×J i and that p + i is continuous in [0, T ]×J i . 2) We assume that for all i ∈ {1, ..., N }, there exists a continuous, quasi-convex and coercive function h i such that for all t ∈ [0, T ] and p ∈ R,

H i (t, 0, p) ≤ h i (p).
For i ∈ {1, ..., N }, there exists an interval [p - hi , p + hi ] such that

     h i is non increasing on (-∞, p - hi ), h i is constant on [p - hi , p + hi ], h i is non decreasing on (p + hi , +∞).
3) For all i ∈ {1, ..., N }, there exists a constant c i such that

c i ≥ max max t∈[0,T ],x∈Ji |p - i (t, x)|, |p - hi |C u
such that:

the function H i is locally Lipschitz continuous in a neighborhood of the set

{(t, x, p) ∈ [0, T ] × J i × R; p ≤ -c i }.
-There exists a continuous, positive and decreasing function ϕ i defined on (-∞, -c i ) such that for some

B i > c i , +∞ Bi 1 sϕ i (-s) ds < +∞ (2.2)
and for almost all t ∈ [0, T ], x ∈ J i and p ∈ R,

((H i ) p • p -H i )(t, x, p) ≥ ϕ i (p) a.e. in {(t, x, p) ∈ [0, T ] × J i × R; p ≤ -c i }. (2.3)
-There exists κ i > 0 such that

|(H i ) p (t, x, p)| ≤ κ i • ((H i ) p • p -H i )(t, x, p) a.e. in {(t, x, p) ∈ [0, T ] × J i × R; p ≤ -c i }.
(H1) Let C ≥ 0 satisfying the following property: for all i ∈ {1, ..., N }, Regularizing effect (H2) For all i ∈ {1, ..., N }, we define the following set

h i (p) ≥ C ≥ min p∈R h i (p) ⇔ p ≥ max t∈[0,T ] p + i (t,
I i = {p ∈ R; p ≤ p - hi and h i (p) ≥ c}. Then, we have that ϕ (h i (p)) h i (p) < ϕ i (C u p)h ′ i (p)p a.
e. in I i . The next two assumptions concern the dependence of the Hamiltonians on the space and the time.

(H3) For all i ∈ {1, ..., N }, there exists

t i > 0 such that H i is locally Lipschitz in x in a neighborhood of the set {(t, x, p); t ∈ [0, t i ], x ∈ J i , p ≤ -η i (t)} and for all t ∈ [0, T ], x ∈ J i and p ∈ R, |(H i ) x (t, x, p)| ≤ 1 2 η i (t) (p(H i ) p -H i ) (t, x, p) a.e in {(t, x, p); t ∈ [0, t i ], x ∈ J i , p ≤ -η i (t)} (H4) For all i ∈ {1, ..., N }, there exists τ i > 0 such that H i is locally Lipschitz in t in a neighborhood of the set {(t, x, p) ∈ [0, τ i ] × J i × R; H i (t, x, p) ≥ η(t)} and for all t ∈ [0, T ], x ∈ J i and p ∈ R, |(H i ) t (t, x, p)| ≤ 1 2 ψ (η(t)) (p(H i ) p -H i ) (t, x, p) a.e in {(t, x, p) ∈ [0, τ i ] × J i × R; H i (t, x, p) ≥ η(t)}
where ψ is defined in (3.3).

Main results

The main result of this paper is given by the following theorem.

Theorem 2.1. Assume (H). Let u be a continuous, non-negative viscosity solution of (1.1). Then, there exists t * and a continuous function η such that for all t ∈ (0, t * ) and x ∈ J, we have ,x) .

u t (t, x) ≥ -η(t)e u(t
A consequence of this result is the following local Lipschitz coninuity of u:

Proposition 2.2. Assume that (H) holds. If u is a non-negative continuous viscosity solution of (1.1), then u is locally Lipschitz continuous on (0, t * ) × J.

Remarks on the assumptions and examples

In this section, we first exhibit some consequences of the assumptions and we construct in particular the function η. We then provide typical examples of Hamiltonians satisfying the assumptions.

consequences of the assumptions

In this subsection, we derive consequences from assumptions (H0)-(H2).

Consequence of (H0). The function

F i : (c i , +∞) → (0, F i (c i )) defined by F i (τ ) = 2 ∞ τ dσ σϕ i (-σ)
is decreasing and invertible. We define for s ∈ (0, F i (c i )) the function η i (s) = F -1 i (s). We can easily check that the following hold:

   η ′ i (s) = - η i (s)ϕ i (-η i (s)) 2 η i is decreasing, positive and η i (0 + ) = +∞. (3.1)

Consequence of (H1):

The function F : (c, +∞) → (0, F (c)) defined by

F (τ ) = 2 ∞ τ dσ σϕ(σ)
is decreasing and invertible. We define for s ∈ (0, F (c)) the function η(s) = F -1 (s). We can easily check that the following hold:

   η ′ (s) = - η(s)ϕ(η(s)) 2 η is decreasing, positive and η(0 + ) = +∞. (3.2) Moreover, x → x + ϕ(x) is increasing from (c, +∞) to (c + ϕ(c), +∞).
We consider its inverse function denoted by ψ. We have that

ψ(x + ϕ(x)) = x for x ∈ (c, +∞). (3.3)
Consequence of (H2). For t ∈ (0, F (c)) and x ∈ J i , we define the continuous function q - i (t, x) by

q - i (t, x) = min{p < p - i (t, x); H i (t, x, p) = η(t)}.
Then, for t ∈ 0, min

F (c), min i∈{1,...,N } F i (c i ) , we have for any i ∈ {1, ..., N }, -η i (t) > q - i (t, 0)C u (3.4)
where C u is defined in (2.1).

Proof of (3.4). Let q - hi (t) ≤ p - hi such that h i (q - hi (t)) = η(t). We remark that q - hi (t) exists because η(t) > c ≥ min p∈R h i (p). In addition, we have

q - hi (t) ≥ q - i (t, 0). (3.5)
In fact, using that η(t) > c ≥ h i (p - i (t, 0)), we have

q - hi (t) = h -1 i (η(t)) < p - i (t, 0) (3.6) where h -1 i is the inverse of the function h i ↾ (-∞,p - h i
) . Moreover, we have

H i (t, 0, q - hi (t)) ≤ h i (q - hi (t)) = H i (t, 0, q - i (t, 0)) = η(t). (3.7)
Using (3.6) and (3.7), we deduce (3.5). We claim that

-η i (t) > q - hi (t)C u . (3.8) Using that η(t) = h i (q - hi (t)) > c ≥ h i (- c i C u ) with - c i C u ≤ p - hi , we get q - hi (t) ≤ - c i C u .
Therefore, using the inverse of η i , (3.8) is equivalent to

t > F i (-q - hi (t)C u ). (3.9)
We have

F i (-q - hi (t)C u ) = 2 +∞ -q - h i (t)Cu dσ σϕ i (-σ)
.

Setting the variable change σ = -q - hi (s)C u , we obtain

F i (-q - hi (t)C u ) = 2 0 t (q - hi ) ′ (s) q - hi (s)ϕ i (q - hi (s)C u )
ds.

Deriving the equality h i (q - hi (s)) = η(s), we get

(q - hi ) ′ (s) = - ϕ(h i (q - hi (s)))h i (q - hi (s)) 2h ′ i (q - hi (s))
.

Replacing the last equality in the integral, we get

F i (-q - hi (t)C u ) = t 0 ϕ(h i (q - hi (s)))h i (q - hi (s)) ϕ i (q - hi (s)C u )h ′ i (q - hi (s))q - hi (s)
ds.

Using assumption (H2), we get (3.9) and (3.8). Finally, (3.5) and (3.8) implies (3.4).

Examples of Hamiltonians

We give now examples of Hamiltonians satisfying assumptions (H). Let us first consider the model Hamiltonians defined for i ∈ {1, ..., N } by

H i (t, x, p) = |A i (t, x)| mi |p| mi -f i (t, x) (3.10) with            m i > 1, f i is non negative, A i and f i are space-time Lipschitz functions, there exists 0 < δ ≤ 1 such that δ < |A i (t, x)| < 1 δ for (t, x) ∈ [0, T ] × J i .
For such Hamiltonians, our assumptions are satisfied for c i = 0 and c = max

t∈[0,T ] |A(t)| where A(t) is the flux limiter. We take ϕ i (s) = (m i -1)δ mi |s| mi , ϕ(s) = (m -1) Cs and h i (p) = |p| mi δ mi and C < (δ 2 C u ) m where m = min i∈{1,...,N } m i and m = max i∈{1,...,N } m i . Moreover, the function η is defined for s ∈ (0, F (c)) by η(s) = 2 C(m -1)s .
For i ∈ {1, ..., N }, the function η i is defined for s > 0 by

η i (s) = δ 2 m i (m i -1)s 1 m i .
Concerning (H3) and (H4), we have the following

((H i ) p • p -H i )(t, x, p) = (m i -1)|A i (t, x)| mi |p| mi + f i (t, x), |(H i ) x (t, x, p)| ≤ m i |(A i ) x (t, x)||A i (t, x)| mi-1 |p| mi + |(f i ) x (t, x)|, |(H i ) t (t, x, p)| ≤ m i |(A i ) t (t, x)||A i (t, x)| mi-1 |p| mi + |(f i ) t (t, x)|.
We can easily verify that (H3) is satisfied if we take t < t i with

η i (t i ) ≥ max 4K 1 m i (m i -1)δ 2mi-1 , 4K 2 (m i -1)δ mi 1 m i +1
where K 1 and K 2 are respectively bounds of

|(A i ) x | and |(f i ) x |. Assumption (H4) is satisfied if we take t < τ i with η(τ i ) ≥ max 4K 3 m i ((m -1) C + 1) (m i -1)δ 2mi-1 , 4K 4 ((m -1) C + 1) (m i -1) 1 2
where K 3 and K 4 are respectively bounds of

|(A i ) t | and |(f i ) t |.
Another example is the following Hamiltonians:

H 1 (p) = e |p| , H 2 (p) = |p| m , H 3 (p) = |p| n -p, H 4 (p) = |p| log(1 + |p|)
with m, n > 1. In this case, H i = h i . For s negative and small enough, we define for i = 1, 2, 3, 4,

ϕ i (s) = φ(h i (s))
with φ(s) = (log(s)) 2 . We remark that for s negative and small enough, we have that

h i (s) > |s| > 1.
Therefore,

+∞ τ dσ σϕ i (-σ) < +∞ τ dσ σ(log(σ)) 2 < +∞.
Moreover, we can easily verify that

             e |p| (|p| -1) -ϕ 1 (p) → +∞ as p → -∞ (m -1)|p| m -ϕ 2 (p) → +∞ as p → -∞ (n -1)|p| n -ϕ 3 (p) → +∞ as p → -∞ |p| 2 1 + |p| -ϕ 4 (p) → +∞ as p → -∞.
Hence, we deduce that for all i ∈ {1, 2, 3, 4}, there exists c i positive and big enough such that (H0) is satisfied for ϕ i (s) = φ(h i (s)). Concerning (H1)-(H2), we can take ϕ(s) = C φ(s) = C(log(s)) 2 with C < C 2 u . There exists c big enough such that (H1) and (H2) are satisfied. In fact, we have

             e |p| (|p| -1) -ϕ(e |p| ) → +∞ as |p| → +∞ (m -1)|p| m -ϕ(|p| m ) → +∞ as |p| → +∞ (n -1)|p| n -ϕ(|p| n -p) → +∞ as |p| → +∞ |p| 2 1 + |p| -ϕ(|p| log(1 + |p|)) → +∞ as |p| → +∞.
We recall that the above inequalities imply that h ′ i (p)p -h i (p) > ϕ(h i (p)) > 0 for |p| big enough. Thus, to prove (H2), it's sufficient to prove that

C φ(h i (p)) < φ(h i (C u p)) in I i .
For c big enough, the last inequality is true for any i ∈ {1, 2, 3, 4} since we have

lim p→-∞ φ(h 1 (C u p)) φ(h 1 (p)) = C 2 u > C and for i = 2, 3, 4, lim p→-∞ φ(h i (C u p)) φ(h i (p)) = 1 > C.

Proof of Theorem 2.1

Adapting the same strategy in [START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF], we consider the function v(t, x) = -e -u(t,x) . To prove Theorem 2.1, we will prove that for all t ∈ (0, t * ) and x ∈ J, we have

v t ≥ -η(t).
This is done in Theorem 4.2. Before to do that, we give some properties of the equation satisfied by v

Viscosity solutions

The function v is a continuous viscosity solution of

v t + G i (t, x, v, v x ) = 0 if (t, x) ∈ [0, T ] × J * i , v t + G A (t, 0, v, v x ) = 0 if t ∈ [0, T ] (4.1)
where

     G i (t, x, v, p) = (-v)H i t, x, - p v , p ∈ R, G A (t, 0, v, p) = (-v) max A(t), max i=1,...,N H - i t, 0, - p i v , p = (p 1 , p 2 , .., p N ) ∈ R N .
The definition of viscosity solutions of (4.1) can be derived from the one of (1.1). For the reader convenience, we state it: for T > 0, set J T = (0, T ) × J. We define the class of test functions on J × (0, T ) by

C 1 (J T ) = {φ ∈ C(J T ) such that the restriction of φ to (0, T ) × J i is C 1 for all i ∈ {1, .., N }}.
Definition 4.1. Let v : J T → R be a continuous function. We say that v is a viscosity sub-solution (resp. super-solution) of (4.1) in J T if for any test function φ ∈ C 1 (J T ) touching v from above (resp. from below) at some point (t 0 , x 0 ) ∈ J T , we have

φ t (t 0 , x 0 ) + G i (t 0 , x 0 , v(t 0 , x 0 ), φ x (t 0 , x 0 )) ≤ 0 (resp. ≥ 0) if x 0 ∈ J * i φ t (t 0 , 0) + G A (t 0 , 0, v(t 0 , 0), φ x (t 0 , 0)) ≤ 0 (resp. ≥ 0) if x 0 = 0.
If v is viscosity sub-solution and viscosity super-solution, we say that v is a viscosity solution.

Properties of G i

We derive properties of the Hamiltonian G i (t, x, v, p) = (-v)H i t, x, p -v from the assumptions (H). The derivative of this function are involved and we will state new conditions (G0)-(G1)-(G2) and (G3). We define first the time t * for which assumptions (H) are satisfied. We set

t * = min T, min i=1,...,N η -1 i (c i ), η -1 (c), min i=1,...,N t i , min i=1,...,N τ i . (4.2)
By simple computations, we can check that the function G i has the following properties.

(G0) Let i ∈ {1, ..., N }. For all t ∈ (0, t * ), if p ≤ -η i (t), then by using (H0), we have that

(G i ) v (t, x, v, p) ≥ ϕ i (p) and |(G i ) p (t, x, v, p)| ≤ κ i • (G i ) v (t, x, v, p) and (G i ) p (t, x, v, p) ≤ 0. (G1) Let i ∈ {1, ..., N }. For all t ∈ (0, t * ), if G i (t, x, v, p) ≥ η(t)
, then by using (H1) we have that

(G i ) v (t, x, v, p) ≥ ϕ(G i (t, x, v, p)) and |(G i ) p (t, x, v, p)| ≤ κ • (G i ) v (t, x, v, p) .
(G2) Let i ∈ {1, ..., N }. For all t ∈ (0, t * ), if p ≤ -η i (t), then by using (H3) we have that

|(G i ) x (t, x, v, p)| ≤ 1 2 η i (t)(G i ) v ((t, x, v, p). (G3) Let i ∈ {1, ..., N }. For all t ∈ (0, t * ), if G i (t, x, v, p) ≥ η(t)
, then by using (H4) we have that

|(G i ) t (t, x, v, p)| ≤ 1 2 ψ(G i (t, x, v, p))(G i ) v (t, x, v, p).

The main result wrote on v

Theorem 2.1 is a direct consequence of the following result on v.

Theorem 4.2. Assume (H).

Let u be a continuous, non-negative viscosity solution of (1.1). Then, for all t ∈ (0, t * ) and x ∈ J, the function v = -e -u is a viscosity solution of (4.1) and it satisfies:

v(t, x) -v(s, x) ≥ -η(s)(t -s) if 0 < s < t < t * and x ∈ J.
To prove Theorem 4.2, we need the following lemma in which we obtain a gradient lower bound on each branch.

Lemma 4.3. Assume (H). Let u be a continuous, non-negative viscosity solution of (1.1).

The function v = -e -u is a viscosity solution of (4.1) and it satisfies the following: for all i ∈ {1, ..., N } and t ∈ (0, t * ), if x, y ∈ J * i such that x = x i • e i and y = y i • e i with x i > y i , then

v(t, x) -v(t, y) ≥ -η i (t)d(x, y).
Proof of Lemma 4.3. Let i ∈ {1, ..., N }. We define

M = sup t∈(0,t * );x,y∈J * i ,xi>yi>0 {v(t, y) -v(t, x) -η i (t)d(x, y)}.
By contradiction, we assume that M > 0. We introduce

M ν,γ = sup t,s∈(0,t * );x,y∈J * i ,xi>yi>0 v(t, y) -v(s, x) -η i (t)d(x, y) - (t -s) 2 2ν - γ t * -t -αd(0, x) - γ d(0, y) .
Since M > 0, we deduce that M ν,γ > 0 for γ, α small enough. Using that -1 ≤ v < 0, we have 0 < M ν,γ ≤ 1 -αd(0, x). This implies that d(0, x) ≤ 1 α < +∞. Moreover, using that M ν,γ is a supremum of a continuous function, we deduce that it's reached at some point (t, s, x, y).

Step 1: t, s ̸ = 0 and t, s ̸ = t * : Using that η i (t) → +∞ as t → 0 and γ t * -t → +∞ if t → t * , we have that t ̸ = 0 and t ̸ = t * for all ν > 0. Moreover, we have

(t -s) 2 2ν ≤ 1. (4.3)
Inequality (4.3) implies that |t -s| → 0 as ν → 0. Thus, for ν small enough, we have that s ̸ = 0 and s ̸ = t * .

Step 2: y i ̸ = 0 and x i > y i for ν small enough. Using that γ d(0, y)

→ +∞ if d(0, y) → 0, we get that d(0, y) > 0 and thus y i > 0. Assume that there exista ν → 0 such that we have

x i = y i . This implies that v(t, x) -v(s, x) > γ t * . ( 4.4)
Using that d(0, x) = x i ≤ 1 α , we deduce that x converges (up to a sub-sequence) as ν → 0.

Denoting by t (resp. x) the common limit of t and s (resp. x and y), and taking ν to zero in (4.4), we obtain

v( t, x) -v( t, x) ≥ γ t * > 0
which gives a contradiction. Therefore, x i > y i for ν small enough.

Step 3: Use of the viscosity inequalities. Using the viscosity inequalities, we obtain

     η ′ i (t)(x i -y i ) + t -s ν + γ (t * -t) 2 + G i t, y, v(t, y), -η i (t) - γ (y i ) 2 ≤ 0 t -s ν + G i (s, x, v(s, x), -η i (t) -α) ≥ 0.
This implies

η ′ i (t)(x i -y i ) + γ (t * -t) 2 + G i t, y, v(t, y), -η i (t) - γ (y i ) 2 -G i (s, x, v(s, x), -η i (t) -α) ≤ 0.
Denoting by t the common limit of t and s and by x (resp. ȳ) the limit of x (resp. y) as ν goes to zero and using the continuity of G i and v, we obtain after taking ν to zero,

η ′ i ( t)(x i -ȳi ) + γ (t * -t) 2 + G i t, ȳ, v 1 , -η i ( t) - γ (ȳ i ) 2 -G i t, x, v 2 , -η i ( t) -α ≤ 0 (4.5)
where

v 1 = v( t, ȳ) and v 2 = v( t, x).
For σ ∈ [0, 1], we define

ξ σ = ( t, x, v 2 , -η i ( t) -α) + σ 0, ȳ -x, v 1 -v 2 , - γ ȳ2 i + α
We argue as if G i is C 1 (otherwise a standard mollification argument allows to reduce to this case). Using (4.5), we have

η ′ i (t)(x i -ȳi ) + γ (t * -t) 2 + 1 0 d dσ G i (ξ σ )dσ ≤ 0. (4.6)
Step 4: Getting the contradiction. We have that

d dσ G i (ξ σ ) = (G i ) x (ξ σ )(ȳ -x) + (G i ) v (ξ σ )(v 1 -v 2 ) + (G i ) p (ξ σ ) - γ ȳ2 i + α = I 1 + I 2 + I 3 with        I 1 = (G i ) x (ξ σ )(ȳ -x) + 1 2 (G i ) v (ξ σ )(v 1 -v 2 ), I 2 = 1 2 (G i ) v (ξ σ )(v 1 -v 2 ) + (G i ) p (ξ σ )α, I 3 = -(G i ) p (ξ σ ) γ ȳ2 i .
Bound for I 1 : We recall that

ξ σ = t, x + σ(ȳ -x), v 2 + σ(v 1 -v 2 ), p with p = -η i ( t) -α + σ - γ ȳ2 i + α .
Using that p ≤ -η i ( t), we can use (G0) and (G2) and we have

   (G i ) v (ξ σ ) ≥ ϕ i (p) ≥ ϕ i (-η i ( t)), |(G i ) x (ξ σ )| ≤ 1 2 (G i ) v (ξ σ )η i ( t)
where we use that ϕ i is decreasing. Using that

v 1 -v 2 > η i ( t)(x i -ȳi ) + γ t * - t + αx i
we obtain that

I 1 ≥ 0.
Bound for I 2 : Using again (G0), we have

(G i ) p (ξ σ ) ≥ -κ i (G i ) v (ξ σ ).
We obtain that

I 2 ≥ 1 2 (G i ) v (ξ σ )(v 1 -v 2 ) -κ i (G i ) v (ξ σ )α ≥ 1 2 ϕ i (-η i ( t))η i ( t)(x i -ȳi ) + (G i ) v (ξ σ ) 1 2 αx i + γ 2t * -κα > 1 2 ϕ i (-η i ( t))η i ( t)(x i -ȳi )
where we use that (G i ) v (ξ σ ) > 0 and γ 2t * -κ i α > 0 for α small enough.

Bound for I 3 : Using (G0), we have (G i ) p (ξ σ ) ≤ 0 and thus I 3 ≥ 0. Using bounds of I 1 , I 2 and I 3 , we obtain in (4.6)

γ (t * -t) 2 + (x i -ȳi ) η ′ i ( t) + 1 2 ϕ i (-η i ( t))η i ( t) ≤ 0
which gives a contradiction.

We turn now to the proof of Theorem 4.2.

Proof of Theorem 4.2. We introduce

M = sup 0<s<t<t * ,x∈J {v(s, x) -v(t, x) -η(s)(t -s)}.
By contradiction, assume that M > 0. We define

M α = sup 0<s<t<t * ,x∈J v(s, x) -v(t, x) -η(s)(t -s) -βη(s) - δ t * -t -αd(0, x) 2 .
We have that M α ≥ M 2 > 0 for β, δ and α small enough. We classically have that M α is reached at some point (s α , t α , x α ). Using that η(s) → +∞ as s → 0 and δ t * -t → +∞ as t → t * , we have s α ̸ = 0 and t α ̸ = t * . If t α = s α , we get a contradiction using that M α > 0. Till the end of this proof, we will denote

v(s α , x α ) = a and v(t α , x α ) = b.
Case 1: x α = 0. For i ∈ {1, ..., N }, we define the constants λ i and γ i that satisfy the following:

λ i ≥ γ i and γ i -b ≥ p + i (t α , 0), G + i (s α , 0, a, λ i ) > η(s α ) -η ′ (s α )(t α -s α + β) and G + i (t α , 0, b, γ i ) < η(t α ) (4.7)
where

G + i (t, 0, v, p) = (-v)H + i t, 0, p -v with -1 ≤ v < 0 and H + i (t, 0, •) is non-decreasing part of H i (t, 0, •). For ε > 0 small, we define M α,ε = sup 0<s<t<t * ,x,y∈ BR (0) {v(s, x) -v(t, y) -η(s)(t -s) -βη(s) - δ t * -t -αd(0, x) 2 - (d(x, y)) 2 2ε -L(s, t, x) -φ(x) + Ψ(y)}
where

     L(s, t, x) = (s -s α ) 2 + (t -t α ) 2 + d(x α , x) 2 , φ(x) = λ i x i if x = x i e i ∈ J i , Ψ(x) = γ i x i if x = x i e i ∈ J i .
We have M α,ε ≥ M α > 0 and M α,ε is reached at some point (s, t, x, y). Classically, we have

d(x, y) → 0 as ε → 0.
Denoting by s (resp. t, x) the limit of s (resp. t, x and y) as ε → 0, we have

M α ≤ M α,ε ≤ M α -L(s, t, x) -φ(x) + ψ(x).
Using that ψ(x) -φ(x) ≤ 0, we obtain that L(s, t, x) = 0 and hence (s, t, x) = (s α , t α , 0).

Case 1.1: x = x i • e i ∈ J * i , y = y j • e j ∈ J j with i ̸ = j. Writing the viscosity sub-solution's inequality, we obtain

η ′ (s)(t -s + β) -η(s) + 2(s -s α ) + G i (s, x, v(s, x), p ε + 2αx i + 2x i + λ i ) ≤ 0 with p ε = x i + y j ε . Using that G i (s, x, v, p) ≥ G + i (s, x, v, p), p ε + 2αx i + 2x i + λ i > λ i we get η ′ (s)(t -s + β) -η(s) + 2(s -s α ) + G + i (s, x, v(s, x), λ i ) ≤ 0.
Taking ε to zero, we obtain a contradiction using (4.7).

Case 1.2: x = 0 and ∃j such that y ∈ J * j .

Writing the viscosity super-solution's inequality, we obtain

G j (t, y, v(t, y), p ε + γ j ) ≥ η(s) + δ (t * -t) 2 + 2(t -t α ) with p ε = - y j ε
. This implies that for ε small enough,

H j t, y, p ε + γ j -v(t, y) > η(s) > η(t) (4.8)
where we use that η is decreasing. We deduce that

     p ε + γ j -v(t, y) > q + j (t, y) or p ε + γ j -v(t, y) < q - j (t, y) (4.9)
where q - j (t, y) and q + j (t, y) are defined by q + j (t, y) = max{p > p + j (t, y); H j (t, y, p) = η(t)} q - j (t, y) = min{p < p - j (t, y); H j (t, y, p) = η(t)}. Using successively Lemma 4.3 and (3.4) , we have

p ε + γ j ≥ -η j (t) > q - j (t, 0)C u which implies for ε small enough p ε + γ j > q - j (t, y)(-v(t, y)).
Therefore, we obtain p ε + γ j -v(t, y) > q - j (t, y). Using (4.8) and (4.9) , we deduce that

p ε + γ j -v(t, y) > q + j (t, y) (4.10)
and G + j (t, y, v(t, y), p ε + γ j ) > η(t). But using that p ε < 0, we have for ε small enough G + j (t, y, v(t, y), p ε + γ j ) ≤ G + j (t, y, v(t, y), γ j ) < η(t) where we use (4.7). This gives a contradiction. Case 1.3: x = 0 and y = 0. Writing the viscosity super-solution's inequality, we obtain

G A (t, 0, v(t, 0), p) ≥ η(s) + δ t 2 * + 2(t -t α ) ≥ max [0,t * ] A(t) + δ t 2 * + 2(t -t α )
with p = (γ 1 , γ 2 , ..., γ N ). This implies max A(t), max i=1,...,N

H - i t, 0, γ i -v(t, 0) ≥ max [0,t * ] A(t) + δ t 2 * + 2(t -t α ). (4.11)
Taking ε to zero, we obtain a contradiction using (4.7). In fact, using that γ i -b ≥ p + i (t α , 0), we have

H - i t α , 0, γ i -b = min p∈R H i (t α , 0, p) ≤ A(t α ).
Therefore (4.11) implies (after taking ε to zero) that

A(t α ) ≥ A(t α ) + δ t 2 *
.

Case 1.4: ∃i such that x, y ∈ J * i . Writing the viscosity sub and super-solution's inequalities, we have

η ′ (s)(t -s + β) -η(s) + 2(s -s α ) + G i (s, x, v(s, x), p ε + 2αx i + 2x i + λ i ) ≤ 0 (4.
h i p ε + γ i -v(t, y) > C.
Using the definition of C in (H1), we obtain that

p ε + γ i -v(t, y) ≥ max τ ∈[0,t * ]
p + i (τ, 0) or

p ε + γ i -v(t, y) ≤ min τ ∈[0,t * ] p - i (τ, 0).
But, from (4.15), we know that p ε + γ i -v(t, y) > p + i (t, 0). Therefore, we deduce that

p ε + γ i -v(t, y) ≥ max τ ∈[0,t * ] p + i (τ, 0).
The coercivity of H i and inequality (4.12) implies that

|p ε + 2αx i + 2x i + λ i | ≤ C tα,sα,o(α)
and in particular

|p ε | ≤ C tα,sα,o(α) + o(α) + 2r + λ i .
We deduce that up to a sub-sequence, lim ε→0 p ε exists and we denote by p = lim ε→0 p ε . Inequalities (4.12) and (4.13) gives after taking ε to zero

δ (t * -t α ) 2 + η ′ (s α )(t α -s α + β)) + G i (s α , 0, a, p + λ i ) -G i (t α , 0, b, p + γ i ) = δ (t * -t α ) 2 + η ′ (s α )(t α -s α + β)) + 1 0 d dσ G i (ξ σ )dσ ≤ 0 (4.17)
where

ξ σ = (t α + σ(s α -t α ), 0, b + σ(a -b), p + γ i + σ(λ i -γ i )) and σ ∈ [0, 1]. We have that d dσ G i (ξ σ ) = (s α -t α )(G i ) t (ξ σ ) + (a -b)(G i ) v (ξ σ ) + (λ i -γ i )(G i ) p (ξ σ ) = I 1 + I 2 + I 3 with          I 1 = (s α -t α )(G i ) t (ξ σ ) + 1 2 (a -b)(G i ) v (ξ σ ), I 2 = 1 2 (a -b)(G i ) v (ξ σ ), I 3 = (λ i -γ i )(G i ) p (ξ σ ).
Our goal is to prove that

I 1 + I 2 + I 3 ≥ ϕ(η(s α )η(s α )) 2 (t α -s α + β). (4.18)
To obtain our result, we need the following lemma whose proof is postponed.

Lemma 4.4. For all σ ∈ [0, 1], we have

η(s α ) ≤ G i (ξ σ ) ≤ η(s α ) + ϕ(η(s α )). ( 4

.19)

We prove now (4.18). We denote by t σ = t α + σ(s α -t α ). Using (4.19) and (G3), we have

I 1 ≥ 1 2 (s α -t α )(G i ) v (ξ σ )ψ(η(t σ )) + 1 2 (a -b)(G i ) v (ξ σ ) ≥ 1 2 (G i ) v (ξ σ ) (ψ(η(t σ ))(s α -t α ) + (a -b)) ≥ 1 2 (G i ) v (ξ σ )(t α -s α )(η(s α ) -ψ(η(t σ ))) (4.20)
where we use that a -b ≥ η(s α )(t α -s α ). Using (4.19), we have

ψ(η(t σ )) < ψ(η(s α )) ≤ ψ(G i (ξ σ )) ≤ ψ(η(s α ) + ϕ(η(s α ))) = η(s α )
where we use that η is decreasing, ψ is increasing and ψ(x + ϕ(x)) = x. Injecting this inequality in (4.20), we get that

I 1 ≥ 0.
Concerning I 2 , we use (G1) to obtain

I 2 ≥ 1 2 η(s α )(t α -s α + β)ϕ(η(s α )).
We will show now that I 3 ≥ 0. Using (4.16), we have that

p + γ i -b ≥ p + i (t α + σ(s α -t α ), 0). This implies that p + γ i + σ(λ i -γ i ) -b -σ(a -b) ≥ p + i (t α + σ(s α -t α ), 0)
and therefore (G i ) p (ξ σ ) ≥ 0. Hence, I 3 ≥ 0. Finally, injecting (4.18) in (4.17), we obtain

δ (t * -t α ) 2 + (t α -s α + β) η ′ (s α ) + η(s α )ϕ(η(s α )) 2 ≤ 0
which gives a contradiction and the proof is now complete.

Case 2:

x α ̸ = 0. In this case, there exists i ∈ {1, ..., N } such that x α = x α,i e i ∈ J * i . We then define M α,ε without the function Ψ(y) -φ(x). The supremum M α,ε is reached at some point (s, t, x, y) with x, y ∈ J * i for ε small. Writing the viscosity inequalities, we get equation (not necessarily the same i) (4.17) with ξ σ = (t α + σ(s α -t α ), x α , b + σ(a -b), p + σ2αx α,i ) and σ ∈ [0, 1]. We get that

d dσ G i (ξ σ ) = (s α -t α )(G i ) t (ξ σ ) + (a -b)(G i ) v (ξ σ ) + (2αx α,i )(G i ) p (ξ σ ) = I 1 + I 2 with    I 1 = (s α -t α )(G i ) t (ξ σ ) + 1 2 (a -b)(G i ) v (ξ σ ), I 2 = 1 2 (a -b)(G i ) v (ξ σ ) + (2αx α,i )(G i ) p (ξ σ ).
We proceed as in case 1.4. The difference comes from the term (2αx α,i )(G i ) p (ξ σ ). Using Lemma 4.4, we can use (G1) and we get

I 2 ≥ η(s α )ϕ(η(s α ))(t α -s α + β) 2 + (G i ) v (ξ σ ) δ 2t * -2αx α,i .
Finally, using that αx α,i → 0 as α → 0, we have δ 2t * -2αx α,i ≥ 0 for α small enough.

We turn now to the proof of Lemma 4.4.

Proof. Firstly, we will show that for σ ∈ [0, 1],

G i (ξ σ ) ≤ η(s α ) + ϕ(η(s α )).
We have

G i (ξ 1 ) ≤ η(s α ) -η ′ (s α )(t α -s α + β) ≤ η(s α ) + ϕ(η(s α ))η(s α ) 2 (t α -s α + β) ≤ η(s α ) + ϕ(η(s α ))η(s α ) 2 1 - δ t * < η(s α ) + ϕ(η(s α )) (4.21)
where we use that 1

≥ a -b ≥ η(s α )(t α -s α ) + δ t *
. By contradiction, assume that there exists

σ ∈ [0, 1) such that G i (ξ σ ) > η(s α ) + ϕ(η(s α )). We define σ = sup{σ ∈ [0, 1] such that G i (ξ σ ) > η(s α ) + ϕ(η(s α ))} .
By (4.21), we have that σ < 1. By definition of σ, we also have

G i (ξ σ ) = η(s α ) + ϕ(η(s α )).
Using the continuity of G i , we deduce that there exists σ ∈ (σ, 1] such that for σ ∈ [σ, σ], we have

G i (ξ σ ) ∈ [η(s α ), η(s α ) + ϕ(η(s α ))] and G i ( σ) < G i (σ).
However, arguing as above, we can prove that d dσ G i (ξ σ ) ≥ 0 for σ ∈ [σ, σ]. This gives a contradiction and therefore G i (ξ σ ) ≤ η(s α ) + ϕ(η(s α )). Let us now prove that

G i (ξ σ ) ≥ η(s α ).
The proof is similar to the above case. We have G i (ξ 0 ) > η(s α ). We define In fact, for the Hamiltonians H i (p) = |p| mi , taking A = 0, we can take c = 0 and in this case η is defined on (0, +∞) and η(+∞) = 0. It may be useful to exploit (4.22) in any future study of the large time behavior of viscosity solutions of (1.1). In fact, having (4.22) in the case of continuous Hamilton-Jacobi equation was a key idea to obtain the large time behavior of the viscosity solutions, see Corollary 6 in [START_REF] Barles | On the regularizing effect for unbounded solutions of firstorder hamilton-jacobi equations[END_REF] and Theorem 10.6 in [START_REF]An introduction to the theory of viscosity solutions for first-order hamilton-jacobi equations and applications[END_REF]. However, we believe that we can not easily extend the aforementioned results. In fact, a crucial idea in [START_REF] Barles | On the large time behavior of solutions of hamiltonjacobi equations[END_REF] was the use of the periodicity assumption to construct the ergodic constant. This property is missing even in the case of Hamiltonian depending only on the gradient variable.

σ = inf{σ > 0 such that G i (ξ σ ) < η(s α ) } . We have σ > 0, G i (ξ σ ) = η(s α ) and G i (ξ σ ) ∈ [η(s α ), η(s α )+ϕ(η(s α ))] for σ ∈ [0, σ]. The inequality G i (ξ 0 ) > G i (ξ σ ) gives a contradiction.

Local Lipschitz bound

The proof of Proposition 2.2 is a consequence of the following lemma.

Lemma 5.1. Let z = z i e i ∈ J * i and s ∈ (0, t * ). Let r > 0 be such that z i -r > 0 and the ball B r (z) ⊂ (0, t * ). There exists a constant K i depending on z i , t * and r such that for all x, y ∈ B r (z) and for all t ∈ B r (s), we have |u(t, x) -u(t, y)| ≤ K i d(x, y).

(5.1)

Proof of Lemma 5.1. From Lemma 4.3, we have for t ∈ (0, t * ), x = x i e i , y = y i e i and x i > y i , u(t, x) -u(t, y) ≥ -η i (t)Kd(x, y) ≥ -η i (s -r)Kd(x, y)

where K = max (τ,l)∈Br(s)×Br(z) e u (τ,l) .

It remains to prove the upper bound. Let P i > max (t,x)∈Br(s)×Br(z)

p + i (t, x) such that for all (t, x) ∈ B r (s) × B r (z), H i (t, x, P i ) > η(s -r)K.

Let (t 0 , y) ∈ B r (s) × B r (z). Let ν > 0 and define M = sup t∈Br(s),x=xiei,x∈Br(z),yi<xi<zi+r u(t, x) -u(t 0 , y) -P i (x i -y i ) -(t -t 0 ) 2 2ν .

We claim that M ≤ 0. Assume by contradiction that M > 0. For δ > 0 small enough, we define M δ = sup t∈Br(s),x=xiei,x∈Br(z),yi<xi<zi+r u(t, x) -u(t 0 , y)

-P i (x i -y i ) - (t -t 0 ) 2 2ν - δ z i + r -x i .
Classically, M δ > 0 for δ small enough and is reached at some point (t, x). If x i = y i , we get δ z i + r < 0. If x i = z i + r, we get a contradiction using that δ z i + r -x i → +∞. We deduce that y i < x i < z i + r. On the other hand, using the continuity of u, there exists a constant C r,s,z such that 0 < M δ ≤ C r,s,z -(t -t 0 ) 2 2ν .

We deduce that t ∈ B r (s) for ν small enough. Writing the viscosity inequality, we get

H i t, x, P i + δ (z i + r -x i ) 2 ≤ t 0 -t ν .

Remark 4 . 5 .

 45 We remark that the function η depends on the time T . The main reasons of this dependence are two: firstly, we need that η(t) > c > max t∈[0,T ] A(t). Secondly, we need to satisfy assumption (H2) in which the bounds of the solution u on [0, T ] are involved. For example, if we consider the model Hamiltonians (3.10) with A i = 1,f i = 0 and A(t) = A > 0, we obtain H i (p) = |p| mi . Although that the Hamiltonians are not space-time depending, we still needc ≥ A and η(s) = 2 C(m -1)s with C < (C u ) m.The constant C is needed to satisfy assumption (H2). However, if u is bounded in R + × {0}, we can construct η without T -dependence (since C u does not depend on T ). Let us mention that if u(•, 0) is bounded for t > 0, then we can obtain our main result (Theorem 4.2) for t > 0 and η(+∞) = 0.(4.22)

Lebanese American University, Department of computer science and mathematics, Byblos campus, P.O. Box 36, Byblos, Lebanon

Normandie Univ, INSA de Rouen, LMI (EA

-FR CNRS 3335), 76000 Rouen, France, 685 Avenue de l'Université, 76801 St Etienne du Rouvray cedex. France

Using Theorem 2.1, we have t 0 -t ν ≤ η(t)K. Therefore, we obtain

We get a contradiction using the definition of P i . Therefore, M ≤ 0, and taking t = t 0 , we get the desired result. Finally, we can take K i = max (P i , η i (s -r)K).

Proof of Proposition 2.2. The proof of this proposition can be easily derived from the previous lemma. In fact, let us consider the ball B r (0) and take x, y ∈ B r (0) such that x ∈ J i and y ∈ J j with i ̸ = j. Using (5.1), we have