Le Rolle
email: mireia.calvo-gonzalez@univ-rennes1.fr

Heart rate complexity analysis in Brugada syndrome during physical stress testing

Keywords: Brugada syndrome, heart rate dynamics, autonomic nervous system, risk stratification, stress test

Symptoms such as ventricular arrhythmias in Brugada syndrome (BS) typically occur at rest, especially during sleep, suggesting that the autonomic nervous system (ANS) function may be relevant in the arrhythmogenesis of the disease. The aim of this work was to assess the ANS response captured by a non-linear heart rate variability (HRV) measure in 69 patients diagnosed with BS, who underwent a standardized physical stress test. Heart rate complexity (HRC) was evaluated by the power-law scaling analysis (β slope) during rest, exercise, recovery and rest post-recovery, in order to discriminate between symptomatic and asymptomatic BS patients. Symptomatic patients showed a significant reduction in HRC in comparison to asymptomatic subjects, after exertion (p=0.015); during the whole recovery period (p=0.023), and, in particular within the passive recovery phase (p=0.025), as well as during rest post-recovery (p=0.022). Based on these results, symptoms could be associated with a lower ANS complexity during the stress test stages where parasympathetic activity is predominant. Therefore, the proposed HRV indicators could be of help in the risk stratification of asymptomatic patients.

Introduction

Brugada syndrome (BS) is an inherited disorder characterized by a typical electrocardiographic pattern presenting a distinct ST-segment elevation in right precordial leads, associated with a high risk for unexpected sudden cardiac death (SCD) due to malignant ventricular arrhythmias in absence of any apparent structural cardiopathy. Since its description in 1992 as a new cardiac disorder [START_REF] Brugada | Right bundle branch block, persistent {ST} segment elevation and sudden cardiac death: A distinct clinical and electrocardiographic syndrome: A multicenter report[END_REF]), BS has raised a great interest because of its high incidence, especially in far eastern countries, and its relation to sudden death in young adults and, less frequently, in infants and children. It has been estimated that BS are responsible for 4-12% of the total amount of SCD and for 20% of SCD in patients without structural cardiopathy [START_REF] Brugada | Right Bundle-Branch Block and ST-Segment Elevation in Leads V1 Through V3 : A Marker for Sudden Death in Patients Without Demonstrable Structural Heart Disease[END_REF][START_REF] Antzelevitch | Brugada Syndrome: Report of the Second Consensus Conference: Endorsed by the Heart Rhythm Society and the European Heart Rhythm Association[END_REF]. Although several indicators have been assessed for the prediction of cardiac events in BS, studies based on the largest clinical series only proved prior symptoms and spontaneous type-1 Brugada ECG pattern as consistent predictors of ventricular fibrillation (VF) episodes [START_REF] Probst | Long-term prognosis of patients diagnosed with Brugada syndrome results from the FINGER Brugada Syndrome Registry[END_REF]. Thus, risk stratification in BS still remains controversial, especially for asymptomatic patients. However, by relating previous symptomatology with the occurrence of cardiac events in BS, the information provided from new indicators in symptomatic and asymptomatic patients can be compared and related to the prediction of VF events.

Since ventricular arrhythmias in BS typically occur at rest and mainly at night, it is commonly assumed that vagal activity may play a determinant role in the arrhythmogenesis of this disorder [START_REF] Matsuo | The circadian pattern of the development of ventricular fibrillation in patients with Brugada syndrome[END_REF][START_REF] Kies | Abnormal Myocardial Presynaptic Norepinephrine Recycling in Patients With Brugada Syndrome[END_REF]). Indeed, previous works, focused on the analysis of the autonomic nervous system (ANS) function using positron emission tomography, have identified sympathetic dysfunctions on BS patients [START_REF] Kies | Abnormal Myocardial Presynaptic Norepinephrine Recycling in Patients With Brugada Syndrome[END_REF][START_REF] Wichter | Cardiac Autonomic Dysfunction in Brugada Syndrome[END_REF][START_REF] Bigi | Significance of cardiac autonomic neuropathy in risk stratification of Brugada syndrome[END_REF][START_REF] Paul | Autonomic Dysfunction in Patients with Brugada Syndrome: Further Biochemical Evidence of Altered Signaling Pathways[END_REF]. Also, many previous works have studied the autonomic function of BS patients through HRV analysis [START_REF] Krittayaphong | Heart rate variability in patients with Brugada syndrome in Thailand[END_REF][START_REF] Hermida | Decreased nocturnal standard deviation of averaged NN intervals[END_REF][START_REF] Kostopoulou | Disorders of the Autonomic Nervous System in Patients With Brugada Syndrome: A Pilot Study[END_REF][START_REF] Nakazawa | Autonomic imbalance as a property of symptomatic Brugada syndrome[END_REF][START_REF] Makimoto | Augmented ST-Segment Elevation During Recovery From Exercise Predicts Cardiac Events in Patients With Brugada Syndrome[END_REF][START_REF] Tokuyama | Deterioration of the circadian variation of heart rate variability in Brugada syndrome may contribute to the pathogenesis of ventricular fibrillation[END_REF][START_REF] Pierre | Abnormal nocturnal heart rate variability and QT dynamics in patients with Brugada syndrome[END_REF]). However, these works are based on linear, time and frequency-domain methods, being unable to capture the complex, non-linear effects on heart rate (HR) regulation. We hypothesize that the evaluation of HR dynamics on BS patients through non-linear fractal methods may provide additional information for risk stratification in this population.

Moreover, the autonomic response can be better assessed by stimulating the ANS through controlled maneuvers such as physical exercise. Exertion provokes a sympathetic activity increase and a parasympathetic inhibition, leading to higher HR. Conversely, post-exercise cardiodeceleration is adjusted by a continued increase in parasympathetic activity [START_REF] Imai | Vagally mediated heart rate recovery after exercise is accelerated in athletes but blunted in patients with chronic heart failure[END_REF], as well as a progressive sympathetic withdrawal [START_REF] Savin | Autonomic contribution to heart rate recovery from exercise in humans[END_REF]. Since traditional HR analysis at rest commonly fails to distinguish between symptomatic and asymptomatic BS patients, we explored HR changes during physical stress testing. Indeed, [START_REF] Makimoto | Augmented ST-Segment Elevation During Recovery From Exercise Predicts Cardiac Events in Patients With Brugada Syndrome[END_REF] published the first work reporting the interest of exercise to predict ventricular arrhythmias in BS patients, focusing on recovery after exertion.

Approaches using the power-law scaling exponent (β) obtained from the regression line of Power Spectral Density (PSD) vs. frequency plot on a bi-logarithmic scale have shown decreased complexity in patients with cardiovascular disorders and sudden cardiac death events [START_REF] Huikuri | Fractal correlation properties of RR interval dynamics and mortality in patients with depressed left ventricular function after an acute myocardial infarction[END_REF][START_REF] Huikuri | Power-law relationship of heart rate variability as a predictor of mortality in the elderly[END_REF][START_REF] Lombardi | Sudden cardiac death: role of heart rate variability to identify patients at risk[END_REF][START_REF] Gomis | Non-linear dynamic analysis of the cardiac rhythm during transient myocardial ischemia[END_REF][START_REF] Maestri | Nonlinear indices of heart rate variability in chronic heart failure patients: redundancy and comparative clinical value[END_REF][START_REF] Tapanainen | Fractal analysis of heart rate variability and mortality after an acute myocardial infarction[END_REF]. However, we are not aware of any study that has analyzed the HR dynamics' complexity (HRC) of symptomatic and asymptomatic BS patients under conditions of physical activity and recovery.

Thus, the aim of this manuscript was to compare, according to symptomatic status, the autonomic cardiac control captured by this non-linear complexity indicator from RR series, in 69 BS patients, for each phase of a standardized physical stress test. The role of complexity assessment would be helpful in risk stratification of asymptomatic BS patients for whom the decision to implant an Implantable Cardioverter Defibrillator (ICD) based on the benefit/risk for the patient is particularly complex.

Methods

Study population

Sixty-nine Brugada syndrome patients underwent a physical exercise stress test while continuous ECG recordings were acquired. They were enrolled in 6 French centres: Rennes, Saint Pierre de la Réunion, Nantes, Bordeaux, Brest and La Rochelle. After approval by the ethics committee of each centre, all participants provided informed consent to take part in the study, which was performed in accordance with the clinical research and ethics recommendations. Patients' age ranged from 19 to 73 years old (45.06 ± 12.79 years old) and 78.3% were males. Table 1 summarizes participants' baseline characteristics, where no significant differences in age and gender between patients included in the symptomatic and asymptomatic groups can be noted. Sixteen patients had the following documented symptoms: syncope (56.25%), cardiac arrest (25%), dizziness (18.75%) and, less frequently, palpitations (12.50%) and nocturnal convulsions (6.25%).

ICD implantation had been performed in 15 of 53 (28.3%) asymptomatic patients, based on a positive EPS (Electrophysiological Study) test, whereas all symptomatic patients were ICD carriers. Since no relevant cardiac episodes occurred during the test, these devices caused no significant effects on the recordings. Among 54 patients (13 were symptomatic) in whom genetic screening was performed, an SCN5A mutation was found in 19 (35.2%), from whom 5 were symptomatic.

Signal acquisition and test

Patients participated in a clinical protocol recommended by the American Heart Association [START_REF] Gibbons | ACC/AHA 2002 guideline update for exercise testing: summary articleA report of the American college of cardiology/American heart association task force on practice guidelines (committee to update the 1997 exercise testing guidelines) 12345[END_REF] described as triangular stress test since the load was increased until it reached the patient's maximum potential. It was performed on a cyclo ergometer (Ergoline 900 Egamed, Piestany, Slovakia) and divided in the following phases:

• Resting phase: supine rest for 8 minutes and rest sitting on the bike for 2 more minutes.

• Exercise phase:

-Warm-up phase: for men, initial load of 50 watts (W); for women, initial load of 30 W, both for 2 minutes. -Incremental exercise phase: for men, initial load of 80 W for 2 minutes and then incrementing 20 W every 2 minutes; for women, initial load of 50 W increasing load 20 W every 2 minutes.

• After exercise phase:

-Recovery phase: * Active recovery phase: for men, fixed load of 50 W; for women, fixed load of 30 W, both for 3 minutes. * Passive recovery phase: total cessation of effort for 3 minutes.

-Rest post-recovery phase: supine rest until reaching 50 minutes of total recording, never exceeding 30 minutes of resting after recovery.

The purpose of the test was to reach at least 80% of the theoretical maximum heart rate, defined by the formula M HR = 220 -age [START_REF] Fox | Physical activity and the prevention of coronary heart disease[END_REF]).

The standard 12-lead ECG recording was collected and analyzed by the central board (Centre Hospitalier Universitaire de Rennes). It was acquired with the Holter monitor (ELA medical, Sorin Group, Le Plessis Robinsson, France) and continuously recorded at rest (from 10 minutes before warm-up phase beginning), during incremental exercise, recovery (6 minutes) and rest post-recovery. Since each test differed in the incremental exercise and rest post-recovery periods, the duration for these phases was fixed to 7 and 10 minutes, respectively, corresponding to the smallest duration found in the clinical series. Sampling frequency was set at 1000 Hz. Heart rate was continuously displayed and blood pressure measured with the sphygmomanometric method at rest, every two minutes during exercise, maximum effort and recovery.

RR series were obtained using a noise-robust wavelet-based algorithm for QRS complex detection and subsequent R-wave peak location [START_REF] Dumont | Improving ECG Beats Delineation With an Evolutionary Optimization Process[END_REF]. Since the result is a non-uniformly sampled RR series, a cubic splines interpolation was applied to signals in order to obtain regularly sampled data at a rate of 4 Hz.

Power-law Scaling Analysis

The RR series in healthy subjects fluctuates in a complex way with a non-uniform pattern that has been associated with fractal behavior and, on a broadband spectrum, with a 1/f-like distribution (inverse power law). The β index measures the powerlaw behavior as the slope of the regression line calculated by the least-squares method on the linear portion of PSD against frequency plots on a bi-logarithmic scale. Several studies have employed the frequency band 0.0001-0.01 Hz to analyze long-term fluctuations when evaluating sudden cardiac death and myocardial infarction [START_REF] Huikuri | Power-law relationship of heart rate variability as a predictor of mortality in the elderly[END_REF][START_REF] Tapanainen | Fractal analysis of heart rate variability and mortality after an acute myocardial infarction[END_REF]. However, for short-term recordings, this technique has also been applied using the 0.003-0.1 Hz band [START_REF] Lombardi | Heart rate variability patterns before ventricular tachycardia onset in patients with an implantable cardioverter defibrillator[END_REF].

PSD was estimated using the Burg method [START_REF] Burg | Maximum Entropy Spectral Analysis[END_REF]) with an optimized order for each signal, based on the Broersen's combined information criterion [START_REF] Broersen | Finite sample criteria for autoregressive order selection[END_REF], to obtain the autoregressive model parameters. Since we analyzed short segments of the RR series, the scaling exponent β was extracted from the frequency band 0.003-0.1 Hz (-2.5 to -1 on a logarithmic scale).

Statistical Analysis

Comparisons between symptomatic and asymptomatic patients during different phases of the test were evaluated by Mann-Whitney U non-parametric tests. The analysis was made using the commercially available statistical program SPSS 21.0 (IBM, Chicago, USA) and setting the level of significance at p < 0.05.

Binary logistic regression models were applied to variables that significantly differed between groups in order to predict the probability of experiencing symptoms, as a function of HRC measures.

Results

Table 2 summarizes the p-values obtained when comparing the β slope analyzed across the different phases of the test between symptomatic and asymptomatic BS patients, as well as the R-squared (R 2 ) as a fitting quality measure between the PSD and its fitted model associated with each phase. In general, linear models reached an acceptable goodness of fit, where R 2 > 0.86. The lowest fitting qualities were found for exercise and after recovery phases, which were then analyzed in smaller segments with higher R 2 , always greater than 0.92. Symptomatic patients showed a significant reduction in HRC, as estimated by the β slope, in comparison to asymptomatic subjects after exertion (p=0.015); during recovery (p=0.023), and more specifically during the passive recovery phase (p=0.025), and at rest post-recovery (p=0.022). Figure 3 reflects how symptomatic patients showed significantly decreased values in both recovery and post-recovery phases, with respect to asymptomatic subjects. Table 3 collects the odds ratio (OR), 95% confidence intervals (CI) and p-values associated with each univariate model applied to variables that significantly differed between groups. During the rest post-recovery phase, β turned to be an independent predictor in univariate analysis with a 9.7% (p=0.030) decreased risk per one point increase. After exercise, β was also a significant predictor of symptoms, decreasing risk by a 18.9% per each point increase (p=0.049).

Finally, table 4 summarizes the mean ± standard deviation values for symptomatic and asymptomatic patients and summarizes p-values obtained in each significant case, showing higher values of HRC for asymptomatic BS patients.

Discussion

In this study, non-linear HRV indicators were compared in 69 BS subjects who took part in a standardized physical stress test, in order to discriminate between symptomatic and asymptomatic patients. Although previous works have already analyzed the autonomic response of BS patients, this study presents, to our knowledge, the first evaluation of HR dynamics captured by short-term power-law scaling analysis, under conditions of physical activity and recovery. This HRV marker has been employed to assess HRC in other cardiac disorders, such as chronic heart failure or transient myocardial ischemia, demonstrating decreased values of complexity in worse cardiac health conditions [START_REF] Huikuri | Fractal correlation properties of RR interval dynamics and mortality in patients with depressed left ventricular function after an acute myocardial infarction[END_REF][START_REF] Huikuri | Power-law relationship of heart rate variability as a predictor of mortality in the elderly[END_REF][START_REF] Lombardi | Sudden cardiac death: role of heart rate variability to identify patients at risk[END_REF][START_REF] Gomis | Non-linear dynamic analysis of the cardiac rhythm during transient myocardial ischemia[END_REF][START_REF] Maestri | Nonlinear indices of heart rate variability in chronic heart failure patients: redundancy and comparative clinical value[END_REF][START_REF] Tapanainen | Fractal analysis of heart rate variability and mortality after an acute myocardial infarction[END_REF].

The results obtained show a significant reduction in β after exertion, during recovery and rest post-recovery in symptomatic patients, with respect to asymptomatic patients. Therefore, subjects with higher risk showed lower complexity in HR dynamics, as already published for other cardiac pathologies.

Those statistically significant variables were then included in univariate logistic regression models in order to estimate their discrimination power. After recovery, β turned to be an independent predictor of symptomatology, supporting the idea that decreased complexity could be related to a greater risk in Brugada syndrome. After exercise, β was also a significant predictor, relating a decreased power-law scaling exponent to a worse prognosis. It should be noted that although other indicators described in the literature assessing HRV non-linear properties were analyzed, differences between symptomatic and asymptomatic BS patients were not as significant as those captured by the β slope.

Despite DFA, the Hurst exponent, Sample Entropy and Higuchi's Fractal Dimension capture similar HRV properties [START_REF] Acharya | Heart rate variability: a review[END_REF][START_REF] Baumert | Scaling characteristics of heart rate time series before the onset of ventricular tachycardia[END_REF], in our clinical series, only the β slope demonstrated a potential impact on the prediction of risk in BS patients.

The autonomic nervous system is governed by complex interactions resulting from feedback loops of non-linear systems, that enable the organism to adapt to abnormal physiological conditions, such as metabolic variations, stress or disease. According to previous works, a reduction in the complexity of physiological control systems occurs with age and disease [START_REF] Goldberger | What is physiologic complexity and how does it change with aging and disease?[END_REF][START_REF] Rabinovich | The role of chaos in neural systems[END_REF][START_REF] Magrans | Multifractal and nonlinear assessment of autonomous nervous system response during transient myocardial ischaemia[END_REF]. Fractal power-law index β decreased significantly in symptomatic patients after exertion and recovery, when parasympathetic activity is predominant. Indeed, [START_REF] Makimoto | Augmented ST-Segment Elevation During Recovery From Exercise Predicts Cardiac Events in Patients With Brugada Syndrome[END_REF] analyzed the parasympathetic reactivation in BS patients by computing heart rate recovery (HRR) as the heart rate decay in the first minute after peak exercise cessation, concluding that a higher parasympathetic activity was related to the occurrence of cardiac events in BS. Thus, recovery stages seem to unmask significant autonomic differences related to symptomatology, also captured by the non-linear HRV indicators here presented.

The results found in our work indicate a loss of cardiac rhythm complexity after exertion and recovery in symptomatic BS patients. This may be explained by the fact that a decreased organism's adaptability leads to a reduced complexity of the underlying biological control systems [START_REF] Smith | Molecular endocrinology and physiology of the aging central nervous system[END_REF]. Moreover, this altered autonomic cardiac system dynamics when vagal tone is predominant confirms previous works that highlight the relevant role of autonomic imbalance in the genesis of life-threatening cardiac arrhythmias in BS [START_REF] Matsuo | The circadian pattern of the development of ventricular fibrillation in patients with Brugada syndrome[END_REF][START_REF] Kies | Abnormal Myocardial Presynaptic Norepinephrine Recycling in Patients With Brugada Syndrome[END_REF]). Indeed, polymorphic ventricular arrhythmias (PVA) and sudden cardiac death (SCD) mainly occur at rest and especially during sleep, when parasympathetic activity is predominant. Data on cardiac ANS analyzed by positron emission tomography has confirmed that BS patients display some autonomic dysfunction, with increased presynaptic norepinephrine (NE) recycling and thus a reduction in the concentration of NE at the synaptic cleft. This imbalance is thought to promote arrhythmogenicity by decreasing intracellular levels of 3-5-cyclic adenosine monophosphate (cAMP) and thus disrupting the subsequent signaling pathway with potential implication for arrhythmogenesis [START_REF] Wichter | Cardiac Autonomic Dysfunction in Brugada Syndrome[END_REF][START_REF] Kies | Abnormal Myocardial Presynaptic Norepinephrine Recycling in Patients With Brugada Syndrome[END_REF]). This may lead to an enlarged transmural electrical dispersion that results in an increased probability of ventricular arrhythmias. In addition, the autonomic imbalance might be aggravated with the presence of physiological downregulation of adrenergic activity, justifying the occurrence of PVA and SCD at rest or during sleep [START_REF] Priori | Clinical and genetic heterogeneity of right bundle branch block and ST-Segment elevation syndrome A prospective evaluation of 52 families[END_REF][START_REF] Yan | Cellular basis for the Brugada syndrome and other mechanisms[END_REF]. Other studies based on the full stomach test demonstrated that vagal stimulation induced by a large meal not only increased distinctive Brugada-like ECG changes but also lead to prognostic tendencies, since patients with a positive full stomach test were more frequently related to life-threatening events than patients with a negative test result [START_REF] Miyazaki | Autonomic and antiarrhythmic drug modulation of ST segment elevation in patients with Brugada syndrome[END_REF][START_REF] Chalvidan | Near fatal electrical storm in a patient equipped with an implantable cardioverter defibrillator for Brugada syndrome[END_REF].

Thus, findings on the present study further support this hypothesis and provide additional evidence for the role of autonomic imbalance in the pathophysiology of the disease. Although we did not find significant differences between groups in metrics characterizing the parasympathetic component that directly relate the loss of complexity with an increase in vagal tone, this study highlights the importance of physical stress testing in order to unmask relevant autonomic changes that could not be captured at rest, before the test.

Therefore, measurement of HRV not only by conventional but also by new nonlinear methods during physical stress testing could be a helpful tool for documenting neural regulation alteration related to symptoms, which can have a potential impact on the progress of risk stratification in BS patients.

However, a relevant limitation of the study is that the clinical value of the analyzed indicators can only be proved if a close relationship between HRV variations and ventricular arrhythmias is identified. Here, symptomatic patients had documented symptoms but since no cardiac events of ventricular origin were induced during the study; we cannot directly relate HRC alterations in symptomatic patients to a greater risk for suffering symptoms.

Conclusion

We compared non-linear HRV indicators in 69 BS patients who took part in a standardized physical stress test in order to distinguish between symptomatic and asymptomatic groups.

Statistically significant differences were found during and after recovery, showing lower HR dynamics' complexity in symptomatic patients under conditions of recovery and rest after recovery.

Although this study is based on a relatively small population of 69 subjects and, thus, conclusions on symptomatic status predictions cannot be inferred, the identified differences between symptomatic and asymptomatic BS patients indicate relevant tendencies of clinical importance. Therefore, significant HRV results could be of potential help in risk stratification of asymptomatic patients for whom the decision to implant a cardioverter defibrillator is still contentious.

Figure 1 .

 1 Figure 1. Typical RR series obtained from a BS patient during stress test, showing each phase: (a) resting, (b) warm-up, (c) incremental exercise, (d) active recovery, (e) passive recovery and (f) rest post-recovery phases.

Figure 2

 2 Figure 2 displays an example of fit for a symptomatic and an asymptomatic patient after exertion, so as to represent the fitting quality, as well as the β slope decay experienced by the symptomatic case with respect to the asymptomatic patient.

Figure 2 .

 2 Figure 2. Linear fit for β extraction after exercise in a symptomatic (left) and in an asymptomatic (right) BS patient.

Figure 3 .

 3 Figure 3. Group average results for β during resting, exercise, recovery and postrecovery phases; for symptomatic and asymptomatic BS patients (*p <0.05).

Table 1 .

 1 Patients' baseline characteristics.

	All patients Symptomatic Asymptomatic p-value
	(n=69)	(n=16)	(n=53)

Table 2 .

 2 P-values and R 2 obtained when comparing symptomatic and asymptomatic patients along the test.

		β	R 2
	Resting	0.320	0.91
	Exercise	0.418	0.86
	Warm-up	0.135	0.92
	Incremental exercise	0.065	0.94
	After exercise	0.015* 0.94
	Recovery	0.023* 0.95
	Active recovery	0.118	0.96
	Passive recovery	0.025* 0.96
	Rest post-recovery	0.022* 0.89

*p < 0.05, Mann-Whitney U test.

Table 3 .

 3 Odds ratios, 95% confidence intervals and p-values associated with significant HRC indicators in univariate analysis.

		Odds Ratio	95% CI	p-value
	β, recovery	0.180	0.028 -1.146	0.069
	β, passive recovery	0.230	0.029 -1.812	0.163
	β, rest post-recovery	0.097	0.012 -0.803 0.030*
	β, after exercise	0.189	0.036 -0.992 0.049*

*p < 0.05, Mann-Whitney U test.

Table 4 .

 4 Mean ± standard deviation, for symptomatic and asymptomatic patients, and associated p-values of significant variables.

		Symptomatic (n=16) Asymptomatic (n=53) p-value
	β, recovery	-1.866 ± 0.258	-1.685 ± 0.356	0.023
	β, passive recovery	-1.773 ± 0.292	-1.655 ± 0.288	0.025
	β, rest post-recovery	-1.853 ± 0.207	-1.651 ± 0.334	0.022
	β, after exercise	-1.706 ± 0.335	-1.468 ± 0.423	0.015

of arrhythmogenesis associated with ST-segment elevation'. Circulation 100(15):1660-1666.

Acknowledgments

This work was supported by the French Ministry of Health (Programme Hospitalier de Recherche Clinique -PHRC Regional). MC thanks la Caixa Foundation (Spain) and DR acknowledges Lefoulon-Delalande Foundation (France) for financial support.


Juliette Blanchet
email: juliette.blanchet@univ-grenoble-alpes.fr

Claire Aly

Theo Vischel

Gérémy Panthou

Youssouph Sané

Mariane Diop

Théo Vischel

Diop Kane

Trend in the Co-Occurrence of Extreme Daily Rainfall in West Africa Since 1950

come   L'archive ouverte pluridisciplinaire

Introduction

The Sahelian region is one of the world's region with the strongest signal of climate change over the last sixty years [START_REF] Hulme | African climate change: 1900-2100[END_REF][START_REF] Dai | The recent Sahel drought is real[END_REF][START_REF] Gallego | An instrumental index of the West African Monsoon back to the nineteenth century: An Instrumental Index of the West African Monsoon[END_REF]. This African band makes the transition between the dry Saharan Desert towards the North and the wet Sudanian region towards the South. It experiences a monsoon period from June to September with very variable and intermittent rainfalls in the form of storms, called meso-scale convective systems [START_REF] Mathon | Mesoscale Convective System Rainfall in the Sahel[END_REF]. The occurrence and intensity of these storms condition the hydrology in the region [START_REF] Lebel | Seasonal cycle and interannual variability of the Sahelian rainfall at hydrological scales[END_REF][START_REF] Vischel | Assessing the water balance in the Sahel: Impact of small scale rainfall variability on runoff. Part 2: Idealized modeling of runoff sensitivity[END_REF]. They can lead to situations of either droughts or floods, with large consequences for such a vulnerable population to climate hazard [START_REF] Tarhule | Damaging Rainfall and Flooding: The Other Sahel Hazards[END_REF][START_REF] Di Baldassarre | Flood fatalities in Africa: From diagnosis to mitigation[END_REF].

Several studies have shown that the recent climate in the Sahelian region can be divided into three periods [e.g. [START_REF] Lebel | Recent trends in the Central and Western Sahel rainfall regime (1990-2007)[END_REF][START_REF] Nicholson | The West African Sahel: A Review of Recent Studies on the Rainfall Regime and Its Interannual Variability[END_REF]. The first one, from 1950 to 1969, corresponds to a relatively humid period, which was followed by a long and intense drought period between the 1970s and 1980s. Since the 1990s the mean annual rainfall levels have slightly increased with respect to the dry period. This has often led to consider the last two decades as a recovery period [e.g. [START_REF] Nicholson | On the question of the âĂIJrecoveryâĂİ of the rains in the West African Sahel[END_REF][START_REF] Sanogo | Spatiotemporal characteristics of the recent rainfall recovery in West Africa[END_REF]. However annual rainfall remains on average much lower it was during the humid period and it features a strong interannual variability. Moreover the recovery is not consistant over the region with some significant contrasts between the Central Sahel (Niger, Burkina, Eastern Mali), which has progressively become wetter, and the Western Sahel (Western Mali, Senegal), which has remained much in deficit [START_REF] Lebel | Recent trends in the Central and Western Sahel rainfall regime (1990-2007)[END_REF]. In the Central Sahel, the strong decadal variability has been associated with significant changes in the amount of rainfall at mesoscale. The rainfall deficit of the dry period is mainly attributable to a significant decrease in the frequency of storms with no significant decrease in intensity [Le [START_REF] Barbé | Rainfall Variability in West Africa during the Years 1950-90[END_REF][START_REF] Bell | Integration of weather system variability to multidecadal regional climate change: The West African Sudan-Sahel zone, 1951-98[END_REF]. On the contrary, over the last two decades, the occurrence of storms did not increase significantly compared to the dry period, but the contribution of the most intense storms to the rainfall regime reached unprecedented levels [START_REF] Panthou | Characterising the space-time structure of rainfall in the Sahel with a view to estimating IDAF curves[END_REF]. The fact that Sahelian rainfall has turned into a more extreme regime is in phase with the intensification of the water cycle expected at the global scale. The rise in intense storms in the Sahel was recently related to the global warming and is thus likely to be a consistent feature of the Sahelian rainfall regime over the long term [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF].

The present study brings new insight into the evolution of extreme precipitation in the Sahel through the statistical analysis of the co-occurrence of extreme daily rainfall, i.e. of daily rainfall exceeding some large level. The co-occurrence is studied for the period 1950-2014 over two contrasted regions in the Western (Senegal) and the Central Sahel.

The results lead us to speculate on possible changes in intense storm properties in light of the decadal variability and regional disparities.

A suitable theoretical framework for assessing the frequency of rainfall exceeding some large level and their evolution is provided by univariate extreme value statistics [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. [START_REF] Panthou | From pointwise testing to a regional vision: An integrated statistical approach to detect nonstationarity in extreme daily rainfall. Application to the Sahelian region[END_REF] uses this framework to show temporal nonstationarity in the local occurrence of extremes in the Sahelian region, i.e. in their marginal distributions. However this gives only indication about the evolution of extreme rainfall at a given site. The spatial (i.e. surfacic) scale is, by definition, missing. To cope with the spatial variability of heavy rainfall, areal data are required, which is often lacking. Thus one usually resorts to the spatial interpolation of point data supplied by raingages [START_REF] De Michele | IDAF (Intensity-duration-areafrequency) curves of extreme storm rainfall: a scaling approach[END_REF][START_REF] Ceresetti | Multiscale Evaluation of Extreme Rainfall Event Predictions Using Severity Diagrams[END_REF][START_REF] Panthou | Characterising the space-time structure of rainfall in the Sahel with a view to estimating IDAF curves[END_REF], whose quality for representing the most extreme events may be questioned [START_REF] Delrieu | Geostatistical radarâĂŞraingauge merging: A novel method for the quantification of rain estimation accuracy[END_REF].

Another way of looking upon the spatial variability of heavy rainfall is to study the probability of different locations to co-experience some extreme event. A suitable framework for this is provided by the theory of max-stable processes for spatial extremes [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]. Max-stable processes have been used recently for studying the co-occurrence in threshold exceedances of wave height [START_REF] Raillard | Modeling extreme values of processes observed at irregular time steps: Application to significant wave height[END_REF], rainfall [START_REF] Thibaud | Threshold modeling of extreme spatial rainfall[END_REF][START_REF] Huser | Space-time modelling of extreme events[END_REF] or river discharges [START_REF] Asadi | Extremes on river networks[END_REF].

However these studies assume temporal stationarity in the probability of co-occurrent extremes. Nonstationarity in the co-occurence of extreme temperature has been studied in [START_REF] Reich | A Hierarchical Model for Serially-Dependent Extremes: A Study of Heat Waves in the Western US[END_REF] but considering temporal co-occurrence, not spatial one. In a spatial context, it has been studied in [START_REF] Westra | Detection of non-stationarity in precipitation extremes using a max-stable process model[END_REF]; [START_REF] Shang | El Niño-Southern Oscillation influence on winter maximum daily precipitation in California in a spatial model[END_REF]; [START_REF] Zhang | Max-stable based evaluation of impacts of climate indices on extreme precipitation processes across the Poyang Lake basin, China[END_REF] for extreme precipitation by using nonstationary marginal distributions in maxstable processes or in [START_REF] Nicolet | Decreasing spatial dependence in extreme snowfall in the French Alps since 1958 under climate change[END_REF]; [START_REF] Shaby | Bayesian spatial extreme value analysis to assess the changing risk of concurrent high temperatures across large portions of European cropland[END_REF] by studying the evolution of the max-stable spatial dependence of extreme snowfall and extreme temperature, respectively. However these studies consider concomitance in annual maxima, whose explanatory power is rather limited because it only informs on whether large maxima tend to co-occur the same year, without accounting for their timing within the year. Studying trend in the co-occurrence of threshold exceedances aims precisely at overcoming this drawback.

The data are presented in Section 2 while Section 3 describes how is modeled and inferred the probability of concomitant exceedances, which is based on an anisotropic max-stable model. Section 4 assesses the probability of concomitant large exceedances in the two regions and their evolution since the 50s in a 20-year moving window approach.

Section 5 discusses the results through a possible link between changes in the probability of concomitant exceedances and changes in the properties of extreme storms.

Data

Two climatologically contrasted Sahelian regions are used in this study: Senegal and the Central Sahel. The instrumented area in Senegal covers a surface of about 600 × 600 km 2 , as displayed in Figure 2. The network comprises 38 stations with daily measurements since 1950. The instrumented area in the Central Sahel is divided into two networks. The first one covers about 1200 × 600 km 2 and straddles Niger, Burkina Faso, Togo and Benin (Figure 2). 44 stations are dispatched over this area, with daily measurements since 1950. The second network is the AMMA-CATCH rainfall network in Niger [Lebel et al., 2009]. It covers about 110 × 120 km 2 in the region of Niamey. It comprises 30 stations with 5' measurements since 1990, which we aggregate at daily scale.

We restrict the analysis to the four monsoon months of June-July-August-September (JJAS) and to the period 1950-2014, which is the period presenting the most complete data. The mean percentage of missing data over this period is 4% in Senegal, with a maximum of 15%, and 3% in the Central Sahel, with a maximum of 12%. Considering all pairs of stations, the mean percentage of days with at least one concurrent missing value is 7% in Senegal and 6% in the Central Sahel. However more and more lack is found since the 90s. The Senegalese network comprises 38 stations with less than 15% of missing data in the 20-year window 1950-1969, but this decreases to 23 stations in 1995-2014. Similarly, the Central Sahelian network comprises 44 stations with less than 15% of missing data in 1950-1969 but 34 such stations in 1995-2014. In the moving window study of Section 4.2, we will only consider, for a given window, the stations comprising less than 15% of missing data over this window, meaning that actually less stations will be used in the most recent period to document co-occurrence. A separate analysis restricted to the stations having less than 15% of missing data on every 20-year window (23 stations in Senegal and 30 stations in the Central Sahel) shows that the results of the present study are not dependent on the evolution of the network though years (not shown). For the small-scall network of AMMA-CATCH, the 1990-2014 period is considered. Within this period, each station has less than 10% of missing data. Due to the shorter measurement period, these data will not be considered in the non-stationarity study of Section 4.2, but the higher spatial density of the network will provide valuable information of cooccurrence for small interdistances in Section 4.1.

Daily rainfall in JJAS can be considered as identically distributed. Moreover the little extension (less than 10km large) and short duration (maximum a few hours) of the strong convective cells induce a very strong spatio-temporal variability of rainfall [see Figure 4 of Vischel and Lebel, 2007, for example]. This implies in particular that daily amounts can be considered as almost independent from one day to another, as validated in [START_REF] Ali | Utilisation optimale de l'information pluviométrique des MCGA aux échelles hydrologiques au Sahel = Optimal use of GCM rainfall information at the hydrological scale in the Sahel[END_REF] and [START_REF] Gerbaux | The sensitivity of Sahelian runoff to climate change / Sensibilite au changement climatique du ruissellement au Sahel[END_REF] for reproducing the interannual and intraseasonal variability of Sahelian rainfields. Serial independence of extremes -which are the main object of this study-can be quantitatively assessed by computing the mean number of consecutive values that exceed high thresholds. Independent exceedances correspond to a mean value of 1 and the closer to 1, the less dependent the exceedances. Using the es-timation method of [START_REF] Ferro | Inference for Clusters of Extreme Values[END_REF], we find that the average mean cluster size is 1.2 in Senegal, 1.14 in Central Sahel and 1.05 for AMMA-CATCH network, which confirms that assuming independence of threshold exceedances is reasonable for the three networks.

As shown in Figure 2, a drawback of the studied networks might be that that none of them samples uniformly the space and the azimuths. The Sahelian stations are not uniformly distributed over the region, with a bigger density of stations close to the Atlantic coast. Stations of the Central Sahelian network are fairly uniformly dispatched but the region features a rectangular shape, with a N-S extension more limited than the E-W one.

However a separate analysis performed on a subset of stations distributed quite uniformly on squared domains in both regions showed that actually the results of this study are quite robust to the shape and spatial distribution of the networks (not shown).

Method

Probability of interest

Let X jt be the random variable of daily rainfall at station j and day t, for t = 1, . . . , T.

Daily rainfalls X j1 , . . . , X jT are assumed to be independent and identically distributed, i.e. to be a random sample of a variable X j representing daily rainfall a given day at station j, with distribution function F j .

Let R be a given region. We are interested in the probability of concomitant exceedance across R, i.e. in pr(X j > x j for all j ∈ R)

(1) for x j s large. Eq. ( 1) is the probability that a large level is exceeded simultaneously everywhere across R, due to a very convective activity making storms to affect all of R the same day. By large we mean here larger than the α-quantile of the marginals distributions, i.e. for x j > u j = F -1 j (α) with α large (e.g. α = 0.97 in Section 4). Note that knowing (1) for all x j large is equivalent to knowing pr(h j (X j ) > x j for all j ∈ R)

for all x j large, with h j any positive bijective function. It is convenient to consider h j so that all h j (X j ) have the same margins. In particular, with h j (.) = -1/log{F j (.)}, the trans-formed variables

Y j = h j (X j ) (3)
are unit-Fréchet distributed, i.e. pr(Y j ≤ y) = e -1/y . This framework turns out to be convenient to work with for spatial extremes, so in the rest of this paper we will be interested in the probability of simultaneous exceedance for the transformed variables Y j , i.e.

pr(Y j > y j for all j ∈ R)

for all y j s large. The probability ( 1) is then simply obtained as (4) taken in y j = -1/log F j (x j ).

By transforming the X j s into Y j s and studying (4) instead of ( 1), we assume that the marginal distributions F j s are known. Actually, since we are interested in ( 1) when all x j are larger than some large quantile u j , only the tail of the marginal distribution is needed.

The theory of univariate extremes [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]) tells us that if u j is large enough, we can assume that peaks-over-threshold follow a GPD (Generalized Pareto Distribution), i.e.

that for x > u j

F j (x) = pr(X j ≤ x) = 1 -(1 -α) 1 + ξ j x -u j σ j -1/ξ j , (5)
where σ j > 0 is the scale parameter and ξ j is the shape parameter. Thus a first step of our analysis is to estimate, for each station j, the GPD parameters σ j and ξ j based on the data exceeding the chosen threshold u j , then use these estimates to transform the data x j exceeding the threshold into y j = -1/log F j (x j ). The data below the threshold are set to v = -1/log α and tagged as "censored". They will be used to estimate the model (4) in a censored likelihood approach; see Section 3.2.4.

Modeling of co-occurrence

Asymptotic dependence and asymptotic independence

Here we describe how is modeled the probability of co-exceeding some large thresholds ( 1). This implies modeling the dependence between variables at large levels. Dependence models for extremes can be split into two families [START_REF] Davison | Geostatistics of Dependent and Asymptotically Independent Extremes[END_REF]: asymptotic dependence models and asymptotic independence models. Asymptotic independence models assume dependence of variables at moderate levels but independence as the levels approach the supremum of the support of the distribution. In particular the probability of one variable being larger than its η-quantile conditionally on another variable being larger than its η-quantile,

χ i j (η) = pr(X j > F -1 j (η)|X i > F -1 i (η)), (6)
converges to 0 as η → 1, whatever the dependence at sub-extremal levels (provided the two variables are not almost surely equal), as illustrated Figure 1. The Gaussian dependence model, for example, belongs to this class as soon as its correlation coefficient is not 1, no matter how close to 1 [START_REF] Coles | Dependence Measures for Extreme Value Analyses[END_REF]. Asymptotic dependence models, on the contrary, assume that some dependence remains at infinite levels, i.e. that χ i j (η) converges to a positive value χ + i j as η → 1, as illustrated Figure 1. In practice distinguishing between weakly asymptotically dependent variables ( χ + i j close to 0) and asymptotically independent variables reveals difficult [START_REF] Coles | Dependence Measures for Extreme Value Analyses[END_REF], and this is confirmed on our data. A separate analysis revealed that, when considering annual maxima, the asymptotically dependent Brown-Resnick process (see Section 3.2.2) and the asymptotically independent version of [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF] are actually pretty equivalent. Comparison of the values of CLIC [START_REF] Davison | Geostatistics of Dependent and Asymptotically Independent Extremes[END_REF] in the stationary and non-stationary cases gave preference to either model depending on the period, in equal proportion, with always very similar values of CLIC. Indeed the estimated probabilities of co-occurrence, e.g. of χ i j (η) in ( 6), are pretty similar whatever model when η lies in the observed range of values, making selection difficult. Extrapolation at much larger levels obviously differ but according to CLIC no clear preference is given to either family of models. In this article, we have chosen to consider the asymptotically dependent case, discarding asymptotically independent models for the sake of concision, but also because, from a physical point of view, we don't see why more and more extreme storms should be less and less extended, as assumed by asymptotically independent models. However we keep in mind that our choice may not be optimal at extrapolation, and so that estimations of e.g. χ i j (η) might be only valid in the observed range of values, in which case the choice of model seems to have little impact.

Brown-Resnick model

Let assume that the daily rainfall variables X j s, or equivalently the Y j s of (3), are asymptotically dependent. Following de Haan and [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], chapter 9, well-founded asymptotically dependent models for maxima belong to the family of max-stable processes. [START_REF] Thibaud | Threshold modeling of extreme spatial rainfall[END_REF] justify that this theory can be used for large exceedances.

Various parametric models of max-stable processes have been proposed in the literature.

The most popular ones are the Smith process, the Schlather process [START_REF] Schlather | Models for Stationary Max-Stable Random Fields[END_REF], the Asymptotically dependent model

θ ij = 1 θ ij = 1.2 θ ij = 1.4 θ ij = 1.6 θ ij = 1.8 θ ij = 2 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 η χ ij (η) Asymptotically independent model Figure 1. Theoretical conditional exceedance probabilities χ i j (η) = pr(X j > F -1 j (η)|X i > F -1 i (η)) for (left)
an asymptotically dependent model corresponding to a max-stable process with extremal coefficient θ i j = 1, 1.2, 1.4, 1.6, 1.8, 2; (right) an asyptotically independent model constructed as the inverted max-stable process of the left panel [START_REF] Wadsworth | Dependence modelling for spatial extremes[END_REF].

Brown-Resnick process [START_REF] Kabluchko | Stationary Max-Stable Fields Associated to Negative Definite Functions[END_REF], the Geometric Gaussian process [START_REF] Davison | Statistical Modeling of Spatial Extremes[END_REF] and the Extremal-t process [START_REF] Opitz | Extremal processes: Elliptical domain of attraction and a spectral representation[END_REF]. Few comparative studies exist in the literature but the overall consensus is to prefer the Brown-Resnick, Extremal-t and Geometric Gaussian processes towards the Smith and Schlather processes [START_REF] Davison | Statistical Modeling of Spatial Extremes[END_REF][START_REF] Gaume | Mapping extreme snowfalls in the French Alps using max-stable processes[END_REF], with actually very similar performances for the three first ones but a slightly better robustness of the Brown-Resnick process [START_REF] Nicolet | A multi-criteria leave-two-out cross-validation procedure for max-stable process selection[END_REF]. This lead us to consider the Brown-Resnick process for this study. Alternative options to maxstable modeling could be to consider generalized Pareto processes [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF].

Huser and Davison [2013] derive a closed-form expression for the multidimensional joint distribution of Brown-Resnick process. In particular, the bivariate distribution function is given by [START_REF] Kabluchko | Stationary Max-Stable Fields Associated to Negative Definite Functions[END_REF]:

pr(Y i ≤ y i , Y j ≤ y j ) = exp - 1 y i Φ a i j 2 + 1 a i j log y j y i - 1 y j Φ a i j 2 + 1 a i j log y i y j , (7)
for y i and y j > v, where Φ is the standard Gaussian cumulative distribution function,

a i j = {2γ(||s i -s j ||)}
γ(||s i -s j ||) = ||s i -s j || λ κ , (8)
with λ > 0 and 0 < κ ≤ 2 unknown parameters to estimate.

Combining (3) and ( 7) gives that, for any η ≥ α, the probability (6) of one variable being larger than its η-quantile conditionally on another variable being larger than its ηquantile is

χ i j (η) = = (1 -2η + η θ i j )(1 -η) -1 , (9)
as illustrated Figure 1.

The scalar θ i j in ( 9) is defined by

θ i j = 2Φ{γ(||s i -s j ||) 1/2 2 -1/2 } ∈ [1, 2]
and is termed the extremal coefficient. For the semivariogram γ defined by ( 8),

θ i j = 2Φ ||s i -s j || κ/2 λ -κ/2 2 -1/2 . ( 10
)
Whatever choice of γ, when θ i j = 1, X i and X j are completely dependent, and χ i j (η) in ( 9) is uniformly 1. When θ i j = 2, X i and X j are independent and χ i j (η) = 1 -η. Besides,

whatever θ i j ∈ [1, 2], χ + i j := lim η→1 χ i j (η) = 2 -θ i j . (11)
Actually, whatever θ i j , χ i j (η) can be approximated by χ + i j with less than 10 -2 error as soon as η ≥ 0.99, and in particular for the JJAS season as soon as we consider exeedances of the yearly return level. However, let us recall that extrapolation of χ i j (η) beyond the range of observed values is very uncertain for our data because it relies on the assumption of asymptotic dependence, which we were unable to clearly validate. Thus, rather than a limit at infinite levels, χ + i j will be interpreted in the rest of this article as the conditional probability of concomitant large exceedances, the term "large" embracing roughly the range of observed values (i.e. up to about the 100-year return level).

Accounting for anisotropy

According to (10) combined with ( 9), the Brown-Resnick is an isotropic model of concomitance: any two pairs of stations at the same distance apart have the same conditional probability of concomitant exceedance, whatever η. This means that the set of sites s j such that χ i j (η) = p, are circles centered on station i. However assumption of isotropy is often violated for hydrometeorological variables, which tend to show increased dependence in the direction of the governing wind fields [START_REF] Ali | Invariance in the Spatial Structure of Sahelian Rain Fields at Climatological Scales[END_REF][START_REF] Blanchet | Spatial modeling of extreme snow depth[END_REF][START_REF] Gaume | Mapping extreme snowfalls in the French Alps using max-stable processes[END_REF][START_REF] Nicolet | Decreasing spatial dependence in extreme snowfall in the French Alps since 1958 under climate change[END_REF]. In order to account for such anisotropy, we replace in ( 8) and ( 10) the isotropic distance by the anisotropic distance defined as

||s i -s j || 2 = (s j -s i ) M M(s j -s i ) (12)
where s i and s j are 2D-coordinates vectors, s is the transpose of s and

M = cos ψ sin ψ -b sin ψ b cos ψ (13)
with b > 1 the elongation coefficient and ψ ∈ [-π 2 , π 2 ] the angle. Now the set of sites s j such that χ i j (η) = p are ellipses centered on station i, with the major axis oriented along ψ and an elongation (ratio of the major and minor axes) equal to b.

Inference

Model (7) involves four parameters: λ and κ driving the strength of dependence along the major axis at a fixed distance, b driving the strength of anisotropy and ψ giving the direction of major dependence. Assuming temporal independence of the y's, but spatial dependence, the log-likelihood of the model is written:

(λ, κ, b, ψ) = T t=1 log g(y 1t , . . . , y N t ) ( 14
)
where g is the multivariate density of Brown-Resnick model. [START_REF] Wadsworth | Efficient inference for spatial extreme value processes associated to log-Gaussian random functions[END_REF] gives a closed form expression for g, however its computation results in a combinatorial explosion [START_REF] Davison | Geostatistics of extremes[END_REF][START_REF] Castruccio | High-Order Composite Likelihood Inference for Max-Stable Distributions and Processes[END_REF]. It is possible to circumvent this issue by making estimation based on the pairwise log-likelihood [START_REF] Varin | AN OVERVIEW OF COMPOSITE LIKELI-HOOD METHODS[END_REF][START_REF] Padoan | Likelihood-Based Inference for Max-Stable Processes[END_REF] 1

(λ, κ, b, ψ) = T t=1 N -1 i=1 N j=i+1 log g i j (y it , y jt ), (15)
where g i j is the bivariate density of (Y i , Y j ), i.e. associated to (7), and N is the number of stations.

In applying pairwise likelihood we must account for the fact that exceedances may occur in both variables, in one variable or in neither, whereas the bivariate density associated to (7) is only valid when both variables exceed the threshold v. To do so we apply the censoring approach described by [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF], Section 8.3.1, and used in the context of threshold exceedances of spatial extremes in [START_REF] Thibaud | Threshold modeling of extreme spatial rainfall[END_REF]; [START_REF] Bacro | Estimation of spatial max-stable models using threshold exceedances[END_REF]; [START_REF] Huser | Space-time modelling of extreme events[END_REF]; [START_REF] Raillard | Modeling extreme values of processes observed at irregular time steps: Application to significant wave height[END_REF], for example. Writing G i j the bivariate distribution, valid only when both variables exceed v, the likelihood contribution of sites i and j in ( 15) is:

g i j (y it , y jt ) =                      ∂ i ∂ j G i j (y it , y jt ), y it > v, y jt > v; ∂ i G i j (y it , v), y it > v, y jt ≤ v, ; ∂ j G i j (v, y jt ), y it ≤ v, y jt > v; G i j (v, v), y it ≤ v, y jt ≤ v; (16)
∂ i G i j denotes the differentiation with respect to y i in (7). No analytical expression for the maximum of ( 15) under ( 16) is available but maximization may be performed numerically (e.g. quasi-Newton method).

Workflow

Two frameworks are considered in this study. The first one assumes stationarity over the observation period. The workflow consists of the following steps:

• For each station j, the empirical α-quantile is computed. The marginal GPD distribution F j above this threshold is estimated by maximum likelihood.

• The data x jt exceeding the threshold are transformed into unit-Fréchet variates

y jt = -1/log F j (x jt ).
For the data x jt below the threshold, the exact value of y jt is not necessary since to compute (16) one only needs to know that y jt is in this case below v.

• The anisotropic Brown-Resnick model for the region R of interest is estimated based on the {y jt , j ∈ R, t = 1, . . . T }, by maximizing the pairwise log-likelihood (15) under (7) with censoring ( 16).

This gives one anisotropic Brown-Resnick model for each of the three networks. Determining the threshold above which this workflow may be applied is a bias-variance tradeoff. Too low a threshold is likely to violate the assumptions of extreme value theory, leading to bias. Too high a threshold implies too few excesses to estimate the model, leading to high variance. In this article, the threshold is set to the empirical 97% quantile of daily rainfall at each station, i.e. α = 0.97. This is slightly higher than the thresholds of [START_REF] Thibaud | Threshold modeling of extreme spatial rainfall[END_REF] (90%-quantile) and Huser and [START_REF] Huser | Space-time modelling of extreme events[END_REF] (95%-quantile) for example but we use here longer time-series, so the number of exceedances is roughly similar. Setting the threshold to the 95%-quantile leads anyway to similar results (not shown).

The second framework considers 20-year moving windows from 1950-1969 to 1995-2014. The workflow consists of the same steps but applied to each window separately (i.e. considering only the days t within the window). This means that the threshold is recomputed for each window and thus different thresholds are used for different windows.

However the models considered here rely on asymptotic theory, so they should apply equally for any large enough threshold. We checked that using for each window the 97%quantile of the entire period does not affect the results, but this implied considering different amounts of data depending on the window, which is not recommended for model comparison and in particular for what matters uncertainty, so we prefer showing results of the window-dependent threshold case. Applying the above workflow gives us one Brown-Resnick model for each network and each window. Given the shortness of its observation period (25 years), the moving-window framework is not applied to the AMMA-CATCH network.

Results

Stationary case

Marginal distributions

Although marginal distributions are not the scope of this paper, they are interesting to study to better understand the climatology of rainfall in the region. Coordinates are UTM coordinates. gal are on average 20% larger than in the Central Sahel. This superiority applies actually for all K ≥ 1 year, with greater differences as the return period increases. The 50-year return level map in the Central Sahel shows a similar pattern to the 20-and 100-year return level maps produced in [START_REF] Panthou | Extreme rainfall in West Africa: A regional modeling[END_REF] over the same area but considering annual maxima. The fact that it differs from the 50-year return level map of Senegal both in terms of magnitude and shape is a novel result.

As shown in Table 1, the proportion of rainfall days (b) explains more the annual totals (a) than the mean of non-zero daily rainfalls (c), as already noted in Le [START_REF] Barbé | Rainfall climatology of the HAPEX-Sahel region during the years 1950-1990[END_REF]; Le [START_REF] Barbé | Rainfall Variability in West Africa during the Years 1950-90[END_REF]. The return level map of the Central Sahel (d) is more correlated to the proportion of rainfall days (b) than to the mean of non-zero daily rainfalls (c), which is more patchy. This suggests that the N-S increase of extreme rainfall intensity is more driven by the occurrence of the rainfall systems than by their own intensity. In Senegal, correlation between the return level map (d) and the proportion of rainfall days (b) is even larger. However the NE-SW gradient of return level maps (d) differs from the N-S gradient characterizing the mean number of rainfall days (b), while the mean of non-zero daily rainfall (c) displays a slight east-west distortion near the coast. This suggests that the spatial organization of extreme rainfall in Senegal might be a combined effect of rainfall occurrence and influence of oceanic moisture entries. This additional local oceanic moisture might also explain why extremes are more intense in Senegal than in the Central Sahel.

Co-occurrence of extremes

Fig. 3 shows the empirical estimates of χ i j (0.99), the conditional probability of coexceeding the 99%-quantile, for all the pairs of stations (i, j) lying in either region. The value of p = 0.99 is chosen for illustration with the aim of being far enough in the tail of the distribution while having enough points to make robust empirical estimations. Remind that the 99% quantile corresponds almost to the 1-year return level in JJAS. Estimates are plotted as a function of the anisotropic distance ( 12) obtained from the estimated Brown-Resnick model of Section 3.2.2 (i.e. with the estimated matrix M in ( 12))

for the 1950-2014 period. Each plot contains N(N -1)/2 points corresponding to the number of pairs of stations. The scatter plot of empirical estimates versus anisotropic distance follows relatively well the curve of χ i j (0.99) predicted by the estimated Brown-Resnick model, which is given by combining ( 9), ( 10), ( 12) and ( 13). Let us recall that the model is fitted on the rainfall data directly, and not on the empirical values of χ i j (0.99) which are uncertain and are only shown here to help judging the quality of the fit. All in all, goodness-of-fit seems satisfactory. Furthermore, in the Central Sahel case, the empirical estimates of the AMMA-CATCH network are also relatively well aligned along the predicted curve, although the estimation is made independently of these data. However the estimated model tends to slightly underestimate the probability of conditional exceedances at very small distance. This may be due to the combination of several factors: first, the fact that AMMA-CATCH network covers a much smaller domain and dependence in the Niamey region may be slightly different than in the whole Central Sahel box. Related to this is the fact that few pairs are located at short distances in the Central Sahel network, and this may induce some bias for estimating the small scale dependence. Add to this the fact that AMMA-CATCH network covers only the most recent period and cooccurrence of extremes may have changed since 1950. This question will be investigated in Section 4.2.

In order to allow comparison of the estimated probability of concomitant exceedances in the two region, we plot in Fig. 4 the values of χ + i j as a function of the difference in coordinates (s js i ), in the two regions, adding the plot obtained for Niamey region from AMMA-CATCH network. Remind that χ + i j is an approximate of the conditional probability of exceeding the K-year return level, for any K ≥ 1 but K ≤ 100 to embrace roughly the observed range of levels. For shortness, it is referred below as the conditional probability of concomitant large exceedances. As mentioned in Section 3.2.3, contours of χ + i j = p delineate circles or ellipses, depending on whether dependence is respectively isotropic or anisotropic. Fig. 4 shows that a slight anisotropy is found in Senegal. The conditional probability of concomitant large exceedances is slightly larger in the E-50-N direction (i.e. bearing 50 • anti-clockwise from East, with ψ of (13) estimated to 50 • ).

Anisotropy is more marked in the Central Sahel region with enhanced concomitance in the E-W direction ( ψ ≈ 0). At a given distance, Senegal shows a slightly larger probability of concomitance than the Central Sahel. For example there is on average 12% chance that a location in Senegal receives some extreme rainfall amount a day when extreme amount is received 100km away (i.e. the average χ + i j is 0.12 at 100km distance). In the Central Sahel, there is on average 9% chance. In Niamey region, we find a marked anisotropy, as in the Central Sahel case, however the direction of maximum probability of concomitance is shifted to E-15-N. The dependence is also larger than in the entire Central Sahel region with e.g. on average 15% chance of conditional large exceedances at 100km distance. As already mentioned, these differences may be due to the different domains, periods of observations and spatial scales the two networks allow to document. show the direction of maximum conditional probability of large exceedances, given by ψ of (13).

Moving window case

Marginal distributions

Fig. 5 shows the evolution of the regional averages of ( found in Senegal and the Central Sahel for the mean values (a to c). In Senegal, three periods are found for the mean annual rainfall (a), with a decrease from the 50s to the 70s, followed by a relative stationarity till the 80s and an increase since then. This is in line with the results of [START_REF] Lebel | Recent trends in the Central and Western Sahel rainfall regime (1990-2007)[END_REF] and [START_REF] Nicholson | The West African Sahel: A Review of Recent Studies on the Rainfall Regime and Its Interannual Variability[END_REF]. The increase since the 80s is due to an increase in both the proportion of wet days (b) and the mean rainfall during wet days (c). Despite the increase in the last decades, the levels of the mean values (a to c) at the end of the period have barely reached the stationary levels. In the Central Sahel, two periods are mainly found for the mean values (a to c). A decrease is found from the 50s to the 70s, as in Senegal, but followed by an increase since the 70s, which is particularly strong for the mean rainfall during wet days (c).

For the 50-year return level (d), similar evolutions are found for the two regions, with two main regimes. First, a decrease between the 50s and the 70s, with return levels 15% smaller in 1968-1977 than they used to be in 1950-1969. This period of decrease is followed by a period of increase in both regions since the 70s-80s, particularly in Senegal. The 50-year return levels in Senegal have exceeded the stationary levels since the 80s and almost reached the levels of the beginning of the period. In the Central Sahel, the stationary levels have barely been recovered so far. These results on extreme rainfalls evolution in the Central Sahel are coherent with [START_REF] Panthou | From pointwise testing to a regional vision: An integrated statistical approach to detect nonstationarity in extreme daily rainfall. Application to the Sahelian region[END_REF] who found a significant decrease of around 10% in the 2-year return level of the periods 1950-1969 and 1970-1990. Over the recent decade, the evolution of the K-year return level in the Central Sahel is also in line with the upward trend in daily rainfall maxima shown by [START_REF] Panthou | Characterising the space-time structure of rainfall in the Sahel with a view to estimating IDAF curves[END_REF].

Co-occurrence of extremes

Fig. 6 shows the ellipses of the contours { χ + i j = 0.1} predicted by the Brown-Resnick model estimated on 20-year moving windows for the two regions. Thus each ellipse delineates the set of points having 10% chance of exceeding some extreme level a day when some extreme level is exceeded in the center of the ellipses. To better assess changes in these ellipses, we also plot in Fig. 6 the length of the minor and major axis of the ellipses. In order to assess uncertainty in these values, a bootstrap method is applied.

For each window and each region, i) we sample the days with replacement; ii) we estimate the marginal GPD distributions on the bootstrap data exceeding the threshold u j used for the original data (thus the threshold itself does not change along the bootstrap procedure);

iii) we transform the data exceeding the threshold into unit-Fréchet variates; iv) we fit the anisotropic Brown-Resnick model on these new data. We apply this 200 times, giving 200 bootstrap estimations for each window. The confidence bands in Fig. 6 show the 0.025and 0.975-quantile of these estimations.

In Senegal, concomitance for the periods 1950-1969 to 1975-1994 is characterized by a preferred direction in the E-50-N direction, which is also the direction found in the stationary case (see Section 4.1.2). However concomitance of the largest rainfalls along this direction was more likely during this period than over all 1950-2014, as shown by the major axis compared to the stationary case in Fig. 6. This anisotropy can be considered as relatively significant since there is little overlap of the 95% confidence bands of the major and minor axis lengths.

For the periods 1976-1995 to 1981-2000, ellipses are less and less elongated, which means that within this period, the preferred direction of concomitance vanished. Finally, from 1982Finally, from -2001Finally, from to 1995Finally, from -2014, concomitance is characterized by a stability in the shape of the ellipses which are roughly isotropic (i.e. the lengths of the major and minor axes are equal). Thus since the 80s concomitance of extreme rainfall is as likely in any direction but weaker than it used to be, as shown by the red ellipses of Fig. 6.

In the Central Sahel, the periods from 1950-1969 to 1981-2000 are characterized by a relative stability in the shape of the ellipses of co-occurrence, with a preferred direction in roughly the E-W direction, as in the stationary case of Section 4.1.2. Judging by the confidence bands in Fig. 6, this preferred direction is highly significant. As for the Senegal, a change in the shape of the ellipses is found in the 80s but surprisingly with opposite consequences. The period from 1982-2001 to 1990-2014 is characterized by a strong enlargement of the ellipses of Fig. 6 in all directions, meaning that, contrary to Senegal, co-occurrence tends to be significantly stronger in all directions these last years.

Conclusions and Discussion

In this paper we study the co-occurrence of extreme daily rainfall and its evolution since 1950 in two Sahelian regions: the Western (Senegal) and the Central Sahel. Our approach is based on a max-stable modeling of rainfall threshold exceedances. This allows us to quantify the probability of experiencing concomitant extremes in these regions and to document, in a moving window approach, whether and to what extend this probability has evolved since 1950.

Some features of the evolution of extreme rainfall regime appear to be consistent for both studied Sahelian regions. They concern the evolution of the marginal distributions of extreme daily rainfall. A change in the 80s is found with a period of decreasing rainfall intensity between the 50s and the 70s, followed by a period of increasing intensity since the 80s. This confirms previous results found for the Central Sahel by [START_REF] Panthou | From pointwise testing to a regional vision: An integrated statistical approach to detect nonstationarity in extreme daily rainfall. Application to the Sahelian region[END_REF] and [START_REF] Panthou | Characterising the space-time structure of rainfall in the Sahel with a view to estimating IDAF curves[END_REF]. It is also in line with the increasing trend in extreme rainfall detected since 1982 by [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF] over the entire Sahel. However some results show contrasted evolutions between Senegal and the Central Sahel. In Sene-gal, a change is found in the 80s in the size and direction of the ellipses delineating the contours of co-occurrence at extreme levels, with preferred co-occurrence in the direction E-50-N before the 80s, and isotropic co-occurrence with lesser extent afterwards. In the Central Sahel, a change is also found in the 80s but with opposite consequences. Ellipses of co-occurrence show a greater extent after the 80s, while anisotropy remains in the E-W direction over the whole period.

These contrasts raise several questions about the evolution of extreme precipitation systems. In the Sahel, most of the rainfall is produced by Mesoscale Convective Systems (MCSs). MCSs are often organized in squall lines propagating from East to West. Under particular conditions, organized convective systems can be stationary (long-lasting systems) and may produce large rainfall amounts over a particular location [START_REF] Lafore | A multi-scale analysis of the extreme rain event of ouagadougou in 2009[END_REF][START_REF] Vischel | Precipitation extremes in the West African Sahel: recent evolution and physical mechanisms[END_REF]. Local convection can also occur, but this generally produces less rainfall. We refer to [START_REF] Laurent | How Important is the Contribution of the Mesoscale Convective Complexes to the Sahelian Rainfall?[END_REF] and [START_REF] Mathon | Mesoscale Convective System Rainfall in the Sahel[END_REF] for a review of MCSs climatology in the Sahel.

The observed evolution of the ellipses of co-occurrence in extremes can be interpreted in light of a change in typology of extreme precipitation systems, assuming some link between the two of them. First, one can postulate that anisotropic ellipses correspond to propagative systems, while isotropic ellipses correspond to stationary systems, and that the more elongated the ellipses, the more propagative the system. Second, larger probabilities of concomitant extremes (larger minor axis and/or bigger ellipses) can be expected to come from bigger extreme rainfall systems.

If these assumptions hold, then the increase in the size of the ellipses since the 80s in the Central Sahel would be a consequence of more propagative and slightly bigger MCSs. This is a plausible hypothesis since the more organized systems in this region [Mesoscale Convective Complex, as defined by [START_REF] Mathon | Mesoscale Convective System Rainfall in the Sahel[END_REF] have such spatial characteristics and are very rainy efficient [START_REF] Laurent | How Important is the Contribution of the Mesoscale Convective Complexes to the Sahelian Rainfall?[END_REF][START_REF] Mathon | Mesoscale Convective System Rainfall in the Sahel[END_REF].

The way we interpret the changes in ellipse features is supported by the recent results of [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF]. They show from Infra-Red (IR) satellite data (period 1982-2016) that the trend in extreme rainfall is associated with a significant increasing trend in the vertical development of the largest Sahelian MCSs, that are known to propagate faster, and to a smaller extent with a trend in their horizontal extension. This evolution of MCS features is retrieved in Fig. 6 where the trend in the major axis of ellipses in the Central Sahel is more pronounced than the trend in their minor axis.

Following the same reasoning in Senegal, the transition between anisotropic to isotropic ellipses could be the manifestation of (i) a transition between a large proportion of propagative systems to a large proportion of stationary systems, (ii) a more erratic propagative direction of systems, and/or (iii) more localized but intense convection.

The reasons for such a contrasted evolution of extreme rainy system features between Senegal and the Central Sahel still have to be determined. The main drivers of the rainfall intensification in the Sahel have been recently identified by [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF].

They incriminate the warming of the Saharan (in link to anthropogenic GHG emissions) since it intensifies convection within MCSs through increased wind shear and intrusion of the hot and dry Saharan air layer at the mid-troposphere. These factors tend to decrease the occurrence of small and localized systems (due to higher convective energy barriers initiating convection) but favor higher development and stronger storms once convection is initiated. While these processes support our results in the Central Sahel, they do not directly explain why storm properties evolved differently in the Western Sahel. [START_REF] Lebel | Recent trends in the Central and Western Sahel rainfall regime (1990-2007)[END_REF] have pointed out the presence of an East-West contrast over the period 1990-2007 with a Western Sahel remaining dryer than the Central Sahel for which a relative increase of mean annual rainfall was observed. This zonal rainfall dipole is also a consistent pattern in GCM rainfall projections [START_REF] Monerie | Expected future changes in the African monsoon between 2030 and 2070 using some CMIP3 and CMIP5 models under a medium-low RCP scenario[END_REF][START_REF] Gaetani | West African monsoon dynamics and precipitation: the competition between global SST warming and CO2 increase in CMIP5 idealized simulations[END_REF] and is expected to accentuate over the 21st century [START_REF] James | Process-based assessment of an ensemble of climate projections for West Africa[END_REF]. The warming of the Sahara is also thought to be responsible for this contrast especially through the enhancement of the Saharan Heat Low activity which simultaneously favors air subsidence in the West -preventing deep moist convection -and accentuates monsoon fluxes in the Central Sahel [START_REF] Monerie | Expected future changes in the African monsoon between 2030 and 2070 using some CMIP3 and CMIP5 models under a medium-low RCP scenario[END_REF][START_REF] James | Process-based assessment of an ensemble of climate projections for West Africa[END_REF][START_REF] Lavaysse | Recent climatological trend of the Saharan heat low and its impact on the West African climate[END_REF].

However the way these changes in monsoon circulation might affect extreme rainy systems properties differently in the West and the Central Sahel remains an open question.

In addition to these synoptic features, more local processes due to land-ocean interface can also not be excluded to have an influence on storm development and propagation.

While the main objective of the present paper was to apply a modern statistical frameworks to characterize the extreme rainfall co-occurrence, the way we interpret the results in light of changes in storm features and atmospheric environment calls for further

Figure 2 .

 2 Figure 2. (a): Annual rainfall, (b): proportion of rainfall days, (c): mean of non-zero daily rainfall and (d): 50-year return level of daily rainfall accumulation (mm) in Senegal and Central Sahel in JJAS season, for the period 1950-2014. The black crosses in the Central Sahelian map show AMMA-CATCH network.

Figure 3 .

 3 Figure 3. Empirical estimates of χ i j (0.99), the conditional probability of co-exceeding the 99%-quantile, for all pairs (i, j) lying in (left) Senegal and (right) the Central Sahel, as a function of the anisotropic distance (12). The red curve shows the values predicted by the estimated Brown-Resnick process for the 1950-2014 period. The horizontal lines at abscissa 1 and 0.01 correspond respectively to cases of complete dependence and independence. In the Central Sahel case, the gray points show the empirical estimates obtained for the AMMA-CATCH network covering the period 1990-2014.

Figure 4 .

 4 Figure 4. Values of χ + i j , the conditional probability of concomitant large exceedances, under Brown-Resnick model, as a function of the difference in UTM coordinates (s is j ), for sites lying in (left) Senegal, (middle) the Central Sahel, and (right) the region of Niamey. In the latter case, estimation is based on the 1990-2014 period while in the two former cases it is based on the 1950-2014 period. The gray dotted lines

  Fig. 5 shows the evolution of the regional averages of (a): mean annual rainfall, (b): mean proportion of rainfall days, (c): mean of non-zero daily rainfall and (d): 50-year return levels, estimated on 20-year moving windows from 1950-1969 to 1995-2014. The 95% confidence intervals are obtained from 200 bootstrap samples, resampling the days with replacement for each window and each station. Surprisingly contrasting results are

Table 1 .

 1 Pattern correlation table of (a): Annual rainfall, (b): proportion of rainfall days, (c): mean of nonzero daily rainfall and (d): 50-year return level of daily rainfall accumulation (mm) in Senegal and Central Sahel in JJAS season, for the period 1950-2014.

	353
	354
	355
	356

Acknowledgments

The research leading to these results has received funding from the NERC/DFID Future Climate For Africa programme under the AMMA-2050 project, grant numbers NE/M020428/1. This work was also supported by the French national programme EC2CO-LEFE "Recent evolution of hydro-climatic hazards in the Sahel: detection and attribution". The AMMA-CATCH rainfall data are freely available at http://bd.amma-catch.org/ except for the last three years which are subject to restricted access. The other daily data are subject to access restrictions according to the data policy defined by the National Weather Services in Senegal, Burkina Faso, Niger, Benin and Togo and are not freely available.

Central Sahel

To that end, the direct analysis of brightness temperature from IR satellite data as in [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF] or the analysis of IR-derived MCS tracking products [as those proposed by Fiolleau and Roca, 2013] might be considered. Further work is also required to better understand the atmospheric mechanisms responsible for the contrasted changes in storms properties. A first option to that end is to make an analytic approach similar to that proposed by [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF], including the analysis of regional contrasts in storms properties within the Sahel. A second possibility is to extend the case studies of extreme events such as [START_REF] Lafore | A multi-scale analysis of the extreme rain event of ouagadougou in 2009[END_REF] and [START_REF] Vischel | Precipitation extremes in the West African Sahel: recent evolution and physical mechanisms[END_REF] in Ouagadougou to other locations in the Sahel [see e.g. [START_REF] Engel | Extreme precipitation in the west african cities of dakar and ouagadougou -atmospheric dynamics and implications for flood risk assessments[END_REF]. Atmospheric simulation experiments, in particular high-resolution convection-permitting simulations, could also be very helpful to find some physical explanations of the observed trends in the co-occurrence of extremes.