Malika Chassan
email: malika.chassan@math.univ-toulouse.fr

Jean-Marc Azaïs
email: jean-marc.azais@math.univ-toulouse.fr

Guillaume Buscarlet

Suard

A proportional hazard model for storm occurrence risk

The aim of this paper is to give a precise estimation of the extreme magnetic storms frequency per time unit (year) throughout a solar cycle. An innovative approach based on a proportional hazard model is developed. Based on the Cox model, this method includes non-stationarity and covariate influence. The model assumes that the number of storms during a cycle is a non-homogeneous Poisson process. The intensity of this process can be expressed as the product of a baseline risk and a risk factor. In the Cox model, the baseline risk is a nuisance parameter. In our model, it is a parameter of interest that will be estimated. The risk factor depends on a covariate, the solar activity index. As in Extreme Value Theory (EVT) and especially in Peaks Over Threshold (POT) modeling, all the high level events are used to make estimations and the results are extrapolated to the extreme level events. This study highlights a strong correlation between the occurrence intensity of magnetic storms and their position on the solar cycle. The model can be used to forecast occurrence intensity for the current 24th solar cycle.

INTRODUCTION

Sun activity has a direct influence on Earth ionosphere and magnetic field. Earth is protected by its magnetosphere and most of the time, solar radiations have a weak impact. However, sometimes Earth is attained by a more important solar flare or a large coronal mass ejection. These phenomena cause severe disturbances of ionosphere properties and are analysed as ionospheric or geomagnetic storms.

Strongest storms can affect the navigation applications as the SBAS (Satellite-based augmentation systems) EGNOS (European Geostationary Navigation Overlay Service), the European counterpart of the WAAS (Wide Area Augmentation System), mainly in terms of integrity and continuity. For example, the period from October 19th to November 7th 2003 was particularly perturbed and the WAAS was disabled during about 30 hours [START_REF]Severe space weather events-understanding societal and economic impacts[END_REF]. This event, called the Halloween event, is composed of several major magnetic storms. Two of them (on 29-30 October) are among the strongest storms since 1932 and they are the first and the second largest storms of the solar cycle 23. During these two storms, the vertical position error limit for the Vertical NAVigation (VNAV) function exceeded the upper limit defined by the Federal Aviation Administration (50 meters) during more than 25 hours, making the WAAS unusable for aircraft precision approaches [START_REF] Weaver | Halloween space weather storms of[END_REF].

The need for investigation of extreme storm risk is demonstrated with these highly critical examples. The aim of this paper is to describe an innovative approach for the estimation of the probability of occurrence of severe ionospheric storms (per time unit (year)).

Many geomagnetic indices can be used to describe the ionosphere magnetic activity. Rifa [START_REF] Rifa | Études des relations entre indices solaires et géomagnetiques[END_REF], shows a high correlation rate between all these indices. Thus, only one index is used to make the majority of analyses, the three-hour ap index (for "planetary amplitude"). This index is obtained using measurements of magnetic field variations from 13 geomagnetic observatories over the world. Thus, the main advantage of this index is its global character. An extended description of this index and a review of the other geomagnetic indices are presented in Section 3.1.

The data used are retrieved from the National Geophysical Data Center of the NOAA (National Oceanic and Atmospheric Administration) [1]. They consist of 80 years of registration, including 7 complete solar cycles, from the 17th (on the general list), which starts on September 1933, to the 23th, which ends on December 2008. The ap index for the first three cycles of the data set is plotted in Figure 1 on page 3.

Intensive storms being scarce, classical statistical methods for probability estimation, such as empirical frequency estimation, are not accurate enough. In many domains, the probability of scarce extreme events can be estimated by extrapolation, using the Extreme Value Theory. Nevertheless, a direct application of EVT is not achievable here. Indeed, a first issue to investigate is the non-stationarity of the studied phenomenon. Since 1843, and the publication of the short article by Schwabe [START_REF] Schwabe | Solar observations during 1843[END_REF], it is known that Sun activity (characterized by the number of solar sunspots) is characterized by 11-year cycles. Corresponding cycles are observable in the ap index data set (Figure 1). Hüsler [START_REF] Hüsler | Extreme values of non-stationary random sequences[END_REF] extends the results of EVT for stationary processes to some non-stationary ones but no general theory exists for non-stationary processes ( [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF], Chapter 6).

However, various non-stationary extreme models are proposed, as the pragmatic approach of Coles (again Chapter 6 in [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]) or the work of Jonathan and Ewans [START_REF] Jonathan | Modeling the seasonality of extreme waves in the gulf of mexico[END_REF]. In this paper, authors model the seasonality of extreme waves in the Gulf of Mexico. Occurrence rate and intensity of storm peak events vary with season (the Generalized Pareto parameters are expressed as a function of seasonal degree using a Fourier form). Our model could have adopted these approaches but the following arguments show that it was not relevant.

The main problem to the classical application of EVT comes from the type of data used for this study. The ap index lies in the set {0, 2, 3, 4, [START_REF] Aalen | Survival and event history analysis. A process point of view[END_REF][START_REF] Andersen | Cox's regression model for counting processes: a large sample study[END_REF][START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF][START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Hall | Large sample optimality of least squares cross-validation in density estimation[END_REF][START_REF] Love | Long-term biases in geomagnetic k and aa indices[END_REF][START_REF] Schwabe | Solar observations during 1843[END_REF]22,27,32,39,48,56,67,80,94,111,132,154,179,207,236, 300, 400}. Thus, it is a discrete variable, taking a finite number of values. On one hand, the discrete case is explored in [START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF] where Ander-son provides an extension of classical convergence results to some discrete variables. However, in practice, working with a discrete support distribution can be troublesome in some applications. As example, Cooley, Nychka and Naveau [START_REF] Cooley | Bayesian spatial modeling of extreme precipitation return levels[END_REF] analyse hydrometric data but they found that the low precision of the measurements makes the data almost discrete. Hence, a sawtooth behavior of the parameter estimators was found with increasing threshold, which hampered to clearly select a value for the threshold. On the other hand, the fact that the ap index takes a finite number of values seems to be incompatible with the application of EVT. Indeed, an assumption used by Hüsler [START_REF] Hüsler | Extreme values of non-stationary random sequences[END_REF] in the study of exceedances behavior is violated when considering a variable with a finite support. Similarly, the result obtained by Anderson [START_REF] Anderson | Extreme value theory for a class of discrete distributions with applications to some stochastic processes[END_REF] assumes that the distribution support is unbounded. Hence, to the best of our knowledge, the application of EVT to the ap index data set is not achievable.

Thus, an innovative approach is developed to estimate the occurrence probability of severe magnetic storms. This model, based on proportional hazard modeling, is also based on EVT. As in the classical Peaks Over Threshold method, the estimation results are obtained using high and less rare events and extrapolated to extreme events.

The main contribution of this paper is the description of a new proportional hazard model, which is detailed in Section 2, along with the Cox method. Definitions of the handled objects and data pretreatments can be found in Section 3. The descriptions of parameter estimators are gathered in Section 4. Section 5 is dedicated to data applications and method extensions.

MODEL DESCRIPTION

The approach developed in this paper is based on the Cox proportional hazard model. Firstly, the general notion of proportional hazard model is presented and details on the particular Cox model are given. Secondly, the counting process formulation of these models is introduced. Finally, our model is described with all the modifications from the initial Cox model.

Cox model and Proportional hazard

First introduced in epidemiology [START_REF] Cox | Regression models and life-tables[END_REF], the Cox model is a survival model. Survival analysis consists in considering an individual (or a sample of individuals) from a starting time (t 0 = 0 without loss of generality) and observing it to detect the occurrence of an event of interest that might happen once.

The occurrence time can be represented by a random variable, T , which also corresponds to survival time (since t 0 = 0). The variable T is used to define the instantaneous risk (or hazard rate) for an individual. The hazard rate at time t, λ(t), is the probability that an event occurs in the infinitesimal interval [t, t + dt] for an individual who has not experienced the event until t:

λ(t) = lim dt→0 1 dt P(t ≤ T < t + dt|T ≥ t), (1)
where P(A) is the probability of the event A.

The Cox model expresses the instantaneous risk with respect to time t and the covariates (X 1 , ..., X p ) as:

λ(t, X 1 , ..., X p ) = λ 0 (t) exp( p i=1 β i X i ), (2)
where λ 0 (t) is the baseline risk.

A consequence of the model ( 2) is that the instantaneous risk is proportional to a baseline risk and that the coefficient is the exponential of a linear combination of the covariates. Thus, considering two individuals with covariates (X

(1) 1 , ..., X (1) p ) and (X (2) 1 , ..., X (2)
p ) the risk factor for the first patient, compared to the second, only depends on covariate effects:

λ(t, X (1) 1 , ..., X (1) p ) λ(t, X (2) 1 , ..., X (2) p ) = exp p i=1 β i (X (1) i -X (2) i ) .
In survival analysis, data used are often censored. In this paper, data are uncensored since only complete solar cycles are considered. Hence, the censoring question will not be developed further. For more details about survival analysis and proportional hazard models, see Aalen, Borgan and Gjessing [START_REF] Aalen | Survival and event history analysis. A process point of view[END_REF].

Counting process

Now, consider a recurrent event that can happen several times for the same individual and the counting process N (t) that counts the number of events up to t. In analogy with Formula (1), we define a proportional intensity model by:

P(N (t + dt) -N (t) = 1) = λ 0 (t) exp( p i=1 β i X i )dt.
Andersen and Gill, in [START_REF] Andersen | Cox's regression model for counting processes: a large sample study[END_REF], define a counting process formulation of the Cox model, including the possibility for the event to be recurrent.

Our model

In this paper, a proportional hazard model is applied. The event of interest is the occurrence of a storm and a solar cycle corresponds to an individual (indexed by the subscript j). Only one covariate is considered, the solar activity index of the cycle j, X j , and the parameter β models its influence. The model constructed has undergone meaningful modifications from the Cox model:

-observations are not censored; -the number of events (storm occurrences) up to time t, N (t), is considered. Hence the counting process is supposed to be an inhomogeneous Poisson process with intensity depending on time and covariates;

-the variable D j (length of cycle j, in years) is included as factor. Thus, the measurement unit is the number of events per time unit and not per cycle;

λ 0 (t) is not considered as a nuisance parameter but as a parameter to estimate;

-as in POT modeling in EVT, estimation is made using all the high level storms and an extrapolation to the storms of extreme level is applied using the parameter P 400 , the probability that a high level storm grows into a storm of level 4001 . Utilization of this parameter assumes that the level reached by a high storm does not depend on the instant of appearance. A chi-square independence test showed that this assumption is acceptable (Appendix A).

Thus, in our model, the number of observed storms (of high level) during cycle j up to time t, called N j (t), is supposed to be a non-homogeneous Poisson process of intensity λ j (t) such as: The baseline intensity λ 0 (t) takes into account the fact that storms occur more likely during the second half of the cycle. It is the main parameter.

λ j (t) = λ 0 (t)D j exp(βX j ), i.e. N j ([a, b]) ∼ P
Throughout this paper λ 0 (t) represents the baseline intensity of high level storms, λ j (t) represents the instantaneous intensity of high level storms, which takes into account the risk factor D j exp(βX j ) (for cycle j). The extrapolated instantaneous intensity for extreme level storms (level 400) is λ 400,j (t), i.e. λ 400,j (t) = λ j (t) × P 400 . Unless the indices j are specified, the parameters refer to a theoretical cycle with a mean solar activity (146.7 in the case of the ap data set).

DATA AND PRETREATMENTS

Prior to model parameter estimation, in this Section, the used data and their pretreatments are illustrated. The choice of the ap index is discussed and definitions of a storm as well as high and extreme levels are given.

Geomagnetic indices

The most commonly used geomagnetic indices are listed below, each with its advantages and disadvantages:

-the K index, introduced by Bartels in 1938, is given with a sampling period of three hours. Each measurement corresponds to the maximal variation of the horizontal component of Earth magnetic field (sum of maximum positive and negative deviations, in nano-Tesla) during the past three-hour interval, compared to a reference quiet day. Each value of the K index is a digit between 0 and 9. As its quasi logarithmic scale limits the use of the K index, especially for calculation of averages, additional indices were created.

-The a index (for "amplitude") is a linear transformation of K. The A index is the average of the 8 a-values per day. The K, a and A indices are computed by various observatories around the world and measurements have to be merged to describe the global behavior of ionosphere.

-Local variations are smoothed by considering weighted averages over the world, represented by the ap and Ap indices, where p stands for "planetary" and the capital letter A refers to the daily average. These two indices are computed using the K index values from 13 observatories distributed over Earth. Weights are fixed according to observatories latitudes. A planetary version of K also exists (Kp) and Table 1 gives the relation between ap and Kp. The Kp and ap indices have been available since 1932.

-The aa index (for "antipodal amplitude") is computed using the K index from only two nearly antipodal geomagnetic stations located in England and Australia. Thus, this index can not be considered as global as the ap index. However, it has the advantage of being available since 1868.

-The Dst (Disturbance storm time) index is restricted to the equatorial magnetic perturbation (see Figure 2). 57 years of data are available (against 80 for the ap index).

-The raw geomagnetic data are available for many geomagnetic observatories. Oldest data date back to 1883 for hourly values and to 1969 for 1 minute values. They consist of measurements of magnetic field variations by means of magnetometers. One disadvantage of these data is the presence of recording gaps (with gap lengths varying from one month to several years depending on the observatory). However, the principal drawback of these data is the level of required pre-treatments. See [1] for more details on geomagnetic indices.

The ap data set has the advantage of being larger than the Dst one and it does not exhibit recording gaps, contrary to raw geomagnetic data. Moreover, the chosen ap index is more global than the aa or Dst indices, as can be seen in Figure 2.

Declustering and storm definition

Ionospheric perturbations are classified in a standardized way using the ap index, according to In this work, a storm is defined by the event with the highest intensity even though a storm normally presents different periods of high intensity separated by less active ones. In order to detect the event of highest intensity, a declustering process (referred to as "runs declustering" in [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF], Chapter 5.3) is applied. Two parameters have to be set:

-a low level, the threshold above which a storm is considered to begin (typically 111, 132 or 154); -the run length r which is the minimal number of observations below the low level between two events for them to be considered as independent.

Thus, two exceedances of the low level separated by less than r measurements are considered to belong to the same cluster (same storm). For each cluster, the storm level is defined as the maximal level reached in the cluster. The first occurrence of this maximum is also recorded, representing the storm date. For a cluster, the storm length is the number of observations between the first up-crossing and the last down-crossing of the low level.

Durations of magnetic storms are highly variable, from 3 or 6 hours for a very high or extreme storm (level 300 or 400) to 90 hours for a low level storm. Due to the declustering in this work, only one-time events are considered (only the first maximum occurrence time is recorded). This is not incoherent since our analyses focus on strong storms, which are brief compared to longer lasting but weaker storms. Nevertheless, the duration of strong storms should be taken into account for the probability of occurrence definition (see Section 3.5).

Data summary

The ap data set contains 7 complete solar cycles and after declustering only 23 magnetic storms of level 400 are detected. The number of events is too low to estimate their frequency as a function of covariate. For storms of level 300, 44 events are detected and this is still insufficient for frequency estimation.

Consequently, as in Extreme Value Theory, less rare events are used to estimate influence of each covariate and an extrapolation of the results to extreme level is performed. In this paper, less rare events are the high level storms (but not necessary extreme). Thus, all the storms of level greater or equal to the low level parameter defined in the declustering process are used. For example, if the low level is 111, a "high level storm" is any storm of level 111, 132, 154, 179, 207, 236, 300 or 400. The "extreme level" is only 400. Mean probability of occurrence for each high level is given in Table 2 Besides the three-hour ap index, a covariate representing the solar activity of a cycle is available and is denoted as X (the unity is the number of sunspots). The solar activity of a cycle is the maximum of monthly Smoothed Sunspot Number (monthly SSN). For an easier interpretation of results, this covariate is centered. See [2] for more details on sunspot number. For the seven cycles of the data set, the solar activity is X = (119.2, 151.8, 201.3, 110.6, 164.5, 158.5, 120.8) with a mean of 146.7.

For each cycle, beginning and end dates are also available, as well as peak date. Peak date corresponds to the date of maximal solar activity within a cycle (characterized by the highest sunspot number). Lengths of cycles are gathered in the vector D, with D j being the length of the cycle j. All the covariates for the ap data set are summarized in Table 3. In addition, beginning date of the current solar cycle (24th) can be found in the last line of Table 3

Time Warping

Durations of the 7 complete solar cycles range from 9.7 to 12.6 years. Thus, in order to analyse all the 7 cycles together, a time warping is applied to each cycle: the position of a storm within a cycle is represented by a number between -0.5 and 0.5 where -0.5 is the beginning of the cycle, 0.5 its end and 0 its middle (peak). In Figure 3, the dash-dotted line represents the warped time for the first complete solar cycle. Notice that the addition of the cycle's length as factor in λ j (t) = λ 0 (t)D j exp(βX j ), assures that the measurement unit is the number of storms per year and not per cycle. Thus, the time standardization may not be biasing the intensity estimation. In the rest of this paper, unless otherwise indicated, warped time is considered.

Precision on probability of occurrence

Due to the application of a declustering process, a storm is now defined by three values: the maximal level, the first time when this maximum is attained and the length of the cluster. Using this modeling, we can estimate the probability: P 1 (t) = P(a storm of level 400 begins at time t), and we want to know the probability: P 2 (t) = P(a storm of level 400 is ongoing at time t).

In the whole data set, the level 400 is reached 29 times, but only 23 storms of level 400 are counted after declustering. Among these 23 storms, 17 reach the level 400 only one time and 6 remain at this level two consecutive measurement times. So, assuming that duration of an extreme storm is never exceeding 6 hours, and denoting p = P(storm stays at the level 400 two times), P 2 (t) can be expressed as:

P 2 (t) ≃ P 1 (t) + P 1 (t -1) × p.
Using a local stationarity argument, we can declare that P 1 (t -1) ≃ P 1 (t). Thus, an approximation of P 2 (t) is:

P 2 (t) ≃ P 1 (t) × (1 + p),
and using the ap data set, this is estimated by:

P 2 (t) ≃ P 1 (t) × (1 + 6/23).

ESTIMATION

P 400 and β

Since P 400 is independent of the position within the cycle (Appendix A), the empirical frequency is used

P 400 = #{storms of level 400} #{storms of level ≥ low level} .
A 95% confidence interval is computed via Gaussian approximation. With m = #{storms of level ≥ low level}, the number of storm with a level above the low level, we have:

P 400 ∈ P 400 ± 1.96 P 400 (1 -P 400 )/m .
For β, we use the fact that

N j = N j ([-0.5, 0.5]) ∼ P 1/2 -1/2 λ 0 (s)ds D j exp(βX j ) .
As in the Cox model, the sufficiency of the statistic N j is verified and β is estimated by its maximum likelihood estimator in a Poisson generalized linear model. A confidence interval is also computed. More details can be found in Appendix B.

Baseline intensity λ 0 (t)

For the estimation of the baseline intensity, a kernel estimator is used. The kernel estimator used here is a slightly modified version of the extension of the Hjort's estimator, proposed by Nielsen and Linton in [START_REF] Nielsen | Kernel estimation in a nonparametric marker dependent hazard model[END_REF]. In our model, the covariate does not depend on time. Thus, some terms of the Nielsen and Linton's estimator become constant. All constant terms are gathered in a normalization constant. Another difference is the fact that in our approach, the parameter of interest is λ 0 (t) instead of λ(t), but the normalization constant is set in order to have a consistent estimator.

Assuming that β is known, the estimator of λ 0 (t) is:

λ 0 (t) = C J j=1 1/2 -1/2 dN j (t -s)φ(s) = C J j=1 1/2 -1/2 N j (t -s)φ ′ (s)ds,
where J is the number of individuals (cycles), C a normalization constant and φ the kernel, verifying φ(±1/2) = 0 (for the integration by parts) and

1/2 -1/2 φ(s)ds = 1.
Bias and variance of this estimator are computed using step functions and by passage to limit. Let φ be a step function,

φ(s) = n i=1 a i 1 A i (s),
where 1 A i is the indicator function of the interval A i , the A i = [t i , t i+1 ] form a partition of [-1/2, 1/2] (we can assume t i < t i+1 without loss of generality) and a i are such that

1/2 -1/2 φ(s)ds = 1. Then, for each t ∈ [-1/2, 1/2], λ 0 (t) = J j=1 C 1/2 -1/2 dN j (t -s)φ(s) = C J j=1 a 1 N j ([t -t 2 , t -t 1 ]) + ... +a n N j ([t -t n+1 , t -t n ]) .
Thus, since N j ([a, b]) ∼ P Q j b a λ 0 (s)ds with Q j = D j exp(βX j ) and since EP(ξ) = VP(ξ) = ξ, we get:

E λ 0 (t) = C J j=1 Q j 1/2 -1/2 λ 0 (s)φ(t -s)ds, (3)
where E is the expected value and V the variance. Similarly, for the variance:

V λ 0 (t) = C 2 N j=1 Q j 1/2 -1/2 λ 0 (s)φ 2 (t -s)ds. (4)
Equations ( 3) and (4) are established for step functions. Extension to general functions is straightforward by a monotone convergence argument.

In case of a kernel concentrated around zero we obtain:

E λ 0 (t) ≃ C J j=1 Q j λ 0 (t).
Hence, the choice C = 1/ Q j is convenient and we get:

V λ 0 (t) ≃ 1 Q j λ 0 (t) 1/2 -1/2 φ 2 (s)ds.
In this paper, a Gaussian kernel is used, i.e.

φ(s) = 1 √ 2πh exp(- s 2 2h 2 ),
with h being the band-width parameter, determined later. Then, using the fact that

φ 2 (s) = 1 2 √ πh φ( √ 2s),
with φ( √ 2s) being the density function of a normal distribution N 0, (h/ √ 2) 2 and assuming sufficiently small h values, it follows that:

1/2 -1/2 φ 2 (s)ds ≃ +∞ -∞ φ 2 (s)ds = 1 2 √ πh .
To avoid edge effects, that make the estimated intensity tending artificially to zero at t = ±0.5, a periodization of data is applied before estimation process. The band-width parameter h is chosen by cross-validation with minimization of Integrated Square Error. See [START_REF] Bowman | An alternative method of cross-validation for the smoothing of density estimates[END_REF] or [START_REF] Hall | Large sample optimality of least squares cross-validation in density estimation[END_REF] for more details. Note that, as indicated in Section 3.5, the intensity estimated by kernel method does not correspond to the intensity we want to evaluate. Instead, it corresponds to the probability P 1 that a storm of level 400 begins at time t. Hence a correction is applied multiplying by λ 0 (t) by 29/23. Thus, the approximate confidence interval for λ 0 (t) is:

λ 0 (t) ∈   λ 0 (t) ± 1.96 1 Q j λ 0 (t) 2 √ πh   .
We recall that Q j = D j exp(βX j ). Thus, in practice, a plug-in estimator of Q j is computed, using β.

Simulations

Simulations are performed in order to assess the quality of estimation method. Sample size (J = number of cycles used) varies from 2 to 300. For each value of J, errors are computed over 700 estimation results. For a fixed number of cycles, a simulation is performed as follow:

-Firstly, X, D and the parameter β * are generated according to values of the ap data set. For example, for β * , we use a normal distribution N (0.006, 0.002 2 ).

-Secondly, an intensity function λ * 0 is created using Gaussian distributions and periodic functions.

-Finally, for each cycle, an inhomogeneous Poisson process of intensity λ * 0 (t)D j exp(β * X j ) is simulated and estimations are performed using our model.

Then, for each value of J, mean relative errors

E β-β * β * 2 and E || λ 0 -λ * 0 || 2 ||λ * 0 || 2 2
are computed and plotted in Figure 4. Relative square error for β is larger than for λ 0 (t) but it decreases more rapidly for small values of J. Before J = 6, the relative square error is on average greater than 30%. Indeed, for a very small number of cycles, the maximum likelihood estimator of β has a large variance and estimation is not very accurate. However, for 7 cycles (as in the ap data set) the mean relative square error is about 0.14, which is acceptable.

In Figure 5, an example of the baseline intensity estimated with our model and using a sample of 10 cycles is illustrated. In accordance with what is shown in Figure 4, one can see that kernel estimation of baseline intensity is rather accurate. Other kernels are tested (uniform, triangular and Epanechnikov) and results are not sensitive to kernel change, for both simulated intensities and intensities estimated from real data.

RESULTS

Baseline intensity λ 0 (t)

Figure 6 depicts the estimation of λ 0 (t) for a low level of 111 and a run length r = 7 including the confidence area (i.e. the intensity for all the storms of level greater or equal to 111, for a theoretical cycle with X = 146.7). The band-width parameter is selected by cross validation and is equal to 0.035. The baseline intensity is higher during the second half of cycle with a significant increase around t = 0, highlighting the difference between the two halves of a solar cycle. This behavior is coherent with empirical observations of frequency of geomagnetic storms (Figure 1).

P 400 and β

For P 400 , the results for different low levels are gathered in With a low level of 111, estimation of β gives: β = 0.0060 with the 95% C.I. [0.0039; 0.0083]. Although this value seems to be small, significance of β has been demonstrated by a likelihood ratio test. The test of β = 0 against β = β returns a p-value of 7.02 × 10 -7 . Thus, solar activity index X affects the number of storms occurring during a cycle. Graphically, the influence of solar activity index on the number of storms per cycle is observable in Figure 7.

λ(t), λ 400 (t) and relative risk

The extrapolation for intensity of extreme level storms is made by multiplying λ(t) by P 400 (with confidence interval). The final intensity obtained, λ 400 , is shown in Figure 8 and corresponds to the intensity of occurrence of extreme storms for a solar cycle with a mean solar activity of 146.7. Because of the properties of proportional hazard models, the occurrence intensity for a cycle of solar activity 146.7 corresponds to a risk factor of 1. Therefore, for a given solar cycle j 0 with a risk factor exp(βX j 0 ), the relative intensity with respect to a solar cycle of solar activity 146.7 is: λ 400,j 0 (t) = λ 0 (t) × P 400 × exp(β(X j 0 -146.7))

Using β and λ 0 , the risk factor of a cycle can be computed and the relative intensity evaluated. For example, compared to the average level of solar activity (146.7), a cycle with a high solar activity of 180 has a risk factor of exp( β(X j 0 -146.7)) = exp(0.0060 × (180 -146.7)) = 1.22 with the 95% C.I. [1. 14; 1.32]. This implies that the extreme storms occurrence probability is expected to be 1.22 times larger during this cycle.

Method sensitivity

Results presented in previous sections are given for a fixed low level (of 111). In the following, the sensitivity of the employed method is evaluated by studying its robustness to low level changes.To do so, results for three other low levels, 132, 154 and 236, are compared to that obtained with a low level of 111 (Figure 9). Figure 9 shows that no significant difference is observable between the three curves, showing that our model is rather stable. In particular, estimated intensities for low levels 132 and 154 are always in the 95% confidence interval of the intensity estimated with a low level of 111. Even if the low level is high (236) the curve shape remains coherent and most of the time within the confidence interval.

Model extensions

Estimations on half-cycles

Since a difference of intensity is clearly observable between the two halfcycles, a modified model is tested, where estimations are performed on each half. The number of storms for a cycle is supposed to be a non-homogeneous Poisson process, with a different intensity for each half and estimations are made separately. However, inconclusive results are returned. Indeed, because of the presence of a normalization constant different for each half, the occurrence intensity during the first half is higher than during the second one. Hence, this approach is not further explored.

Gradient utilization

In an alternative approach, the gradient is used to characterize the strength of a storm (instead of the ap index level). Gradients are computed over one time step (3 hours) and the storm gradient is defined as the maximal gradient attained during a storm. This approach is tested due to the occurrence of unexpected strong effects of low level storms caused by fast variations of the ap index. The same study is performed with this new definition of storm strength. The extreme gradient levels are those greater than 100 and the low one is 35. The estimation of β gives: Resulting values are similar to those obtained with the ap index. The estimated intensity for the storms of extreme gradient is plotted in Figure 10. The step between the two halves of the cycle is more apparent than for the previous estimation (Figure 8). Nevertheless, the gradient utilization exhibits a disadvantage. Since the ap index represents a maximum over a three-hour period, the two values of ap index used for the gradient computation can be separated by nearly 6 hours or only by few minutes. The real dates of these values are not known and the gradient is computed using a three-hour time step. However, the computed gradient gives an approximation of the variation speed of the ap index. Moreover, since gradient is used analogously to the ap index, the original model is still appropriate here.

Application to the aa index

In a further sensitivity study, the proportional hazard model is applied to the aa index.

aa index description

As the ap index, aa is derived from the K index but for aa, only 2 observatories are used instead of 13 as in case of ap (one observatory is located in England and the other one in Australia). Thus, these two indices are correlated ( [START_REF] Rifa | Études des relations entre indices solaires et géomagnetiques[END_REF]) but utilization of aa is motivated by the larger size of the data set. Indeed, this index has been available since 1868, and 12 complete solar cycles, from the 12th (which starts on December 1878) to the 23th, are covered (against 7 for the ap index). Especially due to its derivation process, the utilization of the aa index implies adjustments; the derivation process of aa is as follow:

-Firstly, for a given observatory, each magnetogram measure corresponds to a K value. For example, during an extreme storm with K = 9, the Canberra magnetometer measurement is greater than 450 nT (nanoTesla) and the Hartland one is greater than 500 nT. Then, for each value of K, the scaling factor r K is found using Table 2 in [START_REF] Love | Long-term biases in geomagnetic k and aa indices[END_REF] and does not depend on the observatory. The scaling factor for K = 9 is r K = 667 nT.

-Secondly, weighted average of the two rescaled values is computed. Weights depend on observatory latitude [START_REF] Love | Long-term biases in geomagnetic k and aa indices[END_REF] and since locations of observatories are changing during the considered period, the weights vary too, as indicated in Table 5 -Finally, the mathematical expression of aa is:

aa = 1/2[w Eng .r K Eng + w Aus .r K Aus ]. (5)
For an extreme storm with K = 9 measured after 1980 the corresponding aa index is aa = 1/2[1.059 × 667 + 1.084 × 667] ≃ 714.69 rounded to 715 (the maximal value of aa).

Definition of an extreme storm

The previous relation raises the question of rescaling low and extreme levels used for estimation. A particular attention is taken to the change of magnetic observatories over time, in both England and Australia. These changes are handled by creating two different levels for the periods 1868-1956 (cycles 12 to 18) and 1957-present (cycles 19 to 23) since for each period, weights are quite equivalent. In order to detect all the extreme storms, the lower level has been selected for each period.

An extreme storm is characterised by an ap index of 400 or a K index of 9 according to Table 1. When the corresponding extreme level for aa is computed using Equation ( 5), only 3 storms are detected as extreme since 1878 (the beginning of the first complete solar cycle) and only 2 since 1932 (the period when the ap index is also available). 23 extreme storms are counted with the ap index. These numbers are too different and comparison is not achievable.

Since the aa index is less global than the ap index, the following rule is adopted: considering the time period when the two indices are available, an extreme storm for the ap index should also be detected as extreme for the aa index, even if additional local events are detected. Under this condition, the aa extreme level of 500 seems to be adapted with 23 extreme storms detected for aa against 23 for the ap index (only the second period level is given, the first is automatically computed and is 472). Thus, an extreme storm for aa is defined as a storm of level greater or equal to 500 and the total number of extreme storms detected for the aa index is 33.

Results

The model is then applied to the aa data set with an extreme level of 500, a low level of 130 and a run length of 2. The low level selection is less important than the extreme one because of the model non-sensitivity to this parameter, that remains valid in this application. The low level of 130 is selected similarly to the extreme level. The solar activity index covariate X is now X =(74. The β value is in the 95% confidence interval of the estimation using ap. Since more cycles are available, the confidence interval has been now shortened. Influence of solar activity is a little bit more apparent (Figure 11). The final intensity estimation λ 500 is plotted in Figure 12. The empirical frequency of extreme storms is now 0.23 (since 1878), against 0.29 for the ap data set. This decrease can be explained by addition of solar cycles with a lower solar activity index (the mean solar activity since 1878 is 119.7, against 146.7 for the period 1932 to present).

Figure 12 should be read carefully since intensity curves correspond respectively to a theoretical cycle with a mean solar activity of 146.7 for the ap index and 119.7 for the aa index. Results using the aa index only for the ap index period of availability are also computed (and then correspond to the same mean solar activity) and plotted in Figure 13. The most important result of this paragraph is the similarity of the two intensity curves in Figure 13, showing that the method is stable against a change of data set.

Prevision

The model proposed here also can be used for forecasting the intensity of extreme storms for the current 24th solar cycle. Data needed are: -beginning date; -solar activity index (the X covariate); -peak date; -end date, and so, cycle length D.

Beginning date is known but the three other data must be estimated since the 24th cycle has not ended. NOAA predictions are used and presented in Table 3. Estimation (from the beginning to present (April 2013)) and prediction are represented in Figure 14 (plain line). Solar activity index of the current cycle is estimated at 87.9 and the mean solar activity index of the previous seven cycles is 146.7. Thus the risk factor of the 24th cycle (compared to a cycle with mean solar activity of 146.7) is exp( β(87.9 -146.7)) = 0.70. Thus, occurrence intensity for the current cycle is 0.70 times weaker than the one estimated for a cycle of mean solar activity (gray dashdotted line in Figure 14).

CONCLUSION

In this paper, an ionospheric storm occurrence risk estimation is performed. As extreme storms are rare events, classical methods for probability estimation can not be applied. Moreover, the occurrence probability of ionospheric storms varies over time and the estimation method should take into account this non-stationarity. Based on the Cox model, the construction of an innovative approach using a proportional hazard model has been motivated by the proportional hazard assumptions. In this model, the intensity function can depend on time and a covariate influence can be handled. Mathematically, the number of storms for a cycle is modeled by a non-homogeneous Poisson process of whom the intensity depends on time and one covariate (the solar activity index).

However, extreme events are too scarce to apply straightforward estimations. Then, adopting an Extreme Value Theory process, estimations are made using high events, less rare, and the occurrence intensity across a solar cycle is extrapolated to extreme events, as in the classical Peaks Over Threshold method.

The ap index is chosen to describe ionospheric perturbations. Advantages of this index are its global character, the large amount of data (7 complete solar cycles) and the completeness of time series.

Sensitivity of the model is investigated in two different ways:

-firstly, the model response to a change of low level parameter is examined (this level is used to define the high level storms). Analysis are performed for different low levels. First results of ap index are given for a low level of 111 and a comparison is made using three other low levels: 132, 154 and 236. The shape similarity between all the intensity curves confirms the stability of the method.

-Secondly, the data set used for estimation is changed. In the first approach, the maximal ap value reached during a storm is used to characterise its strength. A second approach considers gradients of storms to classify them as high or extreme events. In a third step, an analogous study is performed using the aa index. All the results obtained are coherent, showing a strong similarity between all the estimated intensity curves.

The main results of this study show a strong correlation between the occurrence intensity of magnetic storms and their positions on the solar cycle. This intensity is higher during the second half of the cycle.

This work also shows that solar activity has an influence on occurrence intensity of extreme storms. With the proportional hazard model, the occurrence intensity of a cycle can be expressed as a relative risk compared to a cycle with a mean solar activity index (146.7 for the ap data set and 119.7 for aa). The proposed model can also be used to forecast the intensity of occurrence of extreme storms for the current 24th solar cycle. level reached by a storm does not depend on its occurrence time. P 400 is then estimated using the empirical frequency. To verify this independence assumption, a Chi-square test is performed.

The warped time interval [-0.5, 0.5] is partitioned into two parts of low and high intensity using the final estimated instantaneous intensity λ 400 . The period of low intensity is defined as T l = {t | λ 400 (t) < 0.30} and the period of high intensity as T h = {t | λ 400 (t) ≥ 0.30}. A threshold of 0.30 is the lowest possible which makes the application of a Chi-square test possible (under the hypothesis of independence the expected number of storms in each cell should be greater or equal to 5).

The numbers of high and extreme storms for T h and T l are gathered in a 2 × 2 contingency table ( The Chi-square statistic is computed for this contingency table and the corresponding p-value (for a significance level of α = 5%) is 0.50. Thus, the independence hypothesis is not rejected. The same test is applied with different thresholds for the warped time partition into T h and T l (0.40, 0.50 and 0.60) and always leads to the same conclusion with p-values of 0.67, 0.39 and 0.61, respectively.

B Maximum likelihood estimator of β

The use of N j (the total number of high level storms in cycle j) instead of N j (t) for the estimation of β raises the question of sufficiency of this statistic. Consider only one cycle and the model: N (t) ∼ P (λ 0 (t)dt D exp(βX)) , for t ∈ [-0.5, 0.5].

Then, consider ∆ 1 , ∆ 2 , ..., ∆ n a partition of [-0.5, 0.5] into n sub-segments. For i = 1...n, note N (∆ i ) = ∆ i dN (t) the number of events in ∆ i . The {N (∆ i ), i = 1...n} are independent variables (by definition of a Poisson process) and N (∆ i ) ∼ P ∆ i λ 0 (s)ds D exp(βX) . Noting C i = ∆ i λ 0 (s)ds D, the Log-likelihood with respect to the counting measure (in which we integrate the weights 1/N (∆ i )! ) is:

exp(βX)

n i=1 C i + n i=1 [N (∆ i ) log(C i )] + βX n i=1 N (∆ i ).
β is linked to the N (∆ i ) only by the term n i=1 N (∆ i ). Hence there is no loss of information when using the total number of events per cycle for the estimation of β.

It is now possible to compute the maximum likelihood estimator. For the j th cycle, the likelihood with respect to the counting measure with weights 1/N j ! is, noting α = 1/2 -1/2 λ 0 (s)ds exp (-α D j exp(βX j )) (α D j exp(βX j )) N j , and the Log-likelihood for all the J cycles: This implicit equation can be resolved only numerically (by the secant method). The Fisher information matrix is also computable:

α -1 J j=1 D j exp(βX j ) J j=1 D j X j exp(βX j ) J j=1 D j X j exp(βX j ) α J j=1 D j X 2 j exp(βX j )

.

The (2,2) coefficient of the inverse matrix of the Fisher information matrix provides the variance of β, used for the construction of a confidence interval.

Figure 1 :

 1 Figure 1: ap index during the first three solar cycles (from the 17th to the 19th). Vertical plain lines represent the change of cycle and dotted lines correspond to the peak of each cycle.

  t)dt , with P(ξ) denoting the Poisson distribution of intensity ξ.

Figure 2 :

 2 Figure 2: Positions of observatories for the Dst (⋆) and the Kp/ap indices ( ).

Figure 3 :

 3 Figure 3: The ap index during the first cycle (from September 1933 to February 1944). The dotted vertical line represents the peak and the dash-dotted line the warped time.

Figure 4 :

 4 Figure 4: Mean relative errors of estimation for λ 0 (t) (top) and β (bottom) as functions of the number of cycles (in logarithmic scale). Mean values are computed over 700 error values.

Figure 5 :

 5 Figure 5: Estimation of λ 0 (t) with a sample of 10 cycles. The curve to estimate is plotted in gray.

Figure 6 :

 6 Figure 6: Estimated baseline intensity (years -1 ), with confidence interval, of the storms of level greater or equal to 111, for a mean solar activity of 146.7.

Figure 7 :

 7 Figure 7: Total number of storms per cycle for a low level of 111, 132 or 154 as a function of solar activity (centered).

Figure 8 :

 8 Figure 8: Estimated instantaneous intensity (years -1 ), with confidence interval, of the storms of level 400 obtained by extrapolation from the low level 111, for a mean solar activity of 146.7. The horizontal dash-dotted line represents the empirical frequency of storms of level 400 (stationary intensity).

Figure 9 :

 9 Figure 9: Estimated instantaneous intensity (years -1 ) of the storms of level 400 obtained by extrapolation from low levels 111 (plain line), 132 (dashdotted line), 154 (dashed line) and 236 (dotted line).

  β = 0.0053 with the 95% C.I. [0.0038; 0.0069].

Figure 10 :

 10 Figure 10: Estimated instantaneous intensity (years -1 ), with confidence interval, of the storms with extreme gradient (≥ 100) obtained by extrapolation from the low gradient level 35, for a mean solar activity of 146.7. The dashdotted line represents the empirical frequency of storms with an extreme gradient.

Figure 11 :

 11 Figure 11: Total number of storms per cycle for low levels of 130, 150 and 170 according to the solar activity (centered).

Figure 12 :

 12 Figure 12: Estimated instantaneous intensity (years -1 ) of the storms defined by aa ≥ 500 (plain line). Estimation is obtained by extrapolation from the low level 130, for a mean solar activity of 119.7. The dotted lines represent the 95% confidence interval. The gray horizontal line corresponds to the empirical frequency of extreme storms. The gray curve with stars is the intensity estimated with the ap data set (for a low level of 111).

Figure 13 :

 13 Figure 13: Similar to Figure 12 with the aa index considered on the 7 complete cycles of the ap data set.

Figure 14 :

 14 Figure 14: Estimation and prediction of instantaneous intensity (years -1 ) of the storms of level 400 for the 24th solar cycle, with confidence interval. For comparison, in dash dotted gray, same intensity for a cycle with a mean solar activity index of 146.7.

DNDD

  j exp(βX j ) + log(α) j log(D j ) + β J j=1 N j X j .The derivatives in α anb β respectively give: j X j exp(βX j ) = J j=1 N j X j .Replacing α by the solution of the first equation, we obtain: J j=1 D j X j exp(βX j ) j exp(βX j ) J j=1 N j X j .

Table 1 [

 1 19]:

	Ionosphere Condition Kp index	ap index
	Quiet	0-1	<7
	Unsettled	2	7 to <15
	Active	3	15 to <27
	Minor storm	4	27 to <48
	Major storm	5	48 to <80
	Severe storm	6	80 to <132
	Large severe	7	132 to <236
	Extreme	8	236 to <400
	Extreme	9	400

Table 1 :

 1 Relation between Kp, ap and ionosphere activity.

Table 2 :

 2 . Number of occurrences and frequency of storms by level.

	Level Number of storms Freq.(×10 4 ) Freq.(year -1 )
	111	182	7.99	2.33
	132	158	6.93	2.02
	154	103	4.52	1.32
	179	84	3.69	1.08
	207	51	2.24	0.65
	236	57	2.50	0.73
	300	44	1.93	0.56
	400	23	1.01	0.29

Table 3 :

 3 . Values of end and peak dates as well as length and solar activity index of this 24th cycle are predictions (in gray), since it is not ended. End of the 24th cycle is estimated around December 2019. The NOAA values are used[4].

	Cycle Beginning	Peak	X	D	Extreme
						storms
	17	Sep 1933 Apr 1937 119.2 10.4	5
	18	Feb 1944 May 1947 151.8 10.2	2
	19	Apr 1954 Mar 1958 201.3 10.5	7
	20	Oct 1964 Nov 1968 110.6 11.7	3
	21	Jun 1976 Dec 1979 164.5 10.3	2
	22	Sep 1986	Jul 1989 158.5 9.7	1
	23	May 1996 Mar 2000 120.8 12.6	3
	24	Dec 2008 Nov 2013 87.9	11	-

Beginning and peak dates, solar activity index X (in sunspots), length D (in years) and number of extreme storms for each cycle from the 17th to the 24th (predicted values are written in gray).

Table 4 .

 4

	Low level	111	132	154
	P in %	3.14	4.19	5.93
	95 % C.I. [1.85 ; 4.43] [2.48 ; 5.90] [3.53 ; 8.33]

Table 4 :

 4 P 400 (probability for a high storm to grow into a storm of level 400) and 95 % confidence interval for each low level.

Table 5 :

 5 . Sources of aa index for different time intervals and corresponding weights.

	Period	British obs.	Australian obs.
		(w Eng = corr. wt.) (w Aus = corr. wt.)
	1868-1919	Greenwich	Melbourne
		(1.007)	(0.967)
	1920-1925	Greenwich	Toolangi
		(1.007)	(1.033)
	1926-1956	Albinger	Toolangi
		(0.934)	(1.033)
	1957-1979	Hartland	Toolangi
		(0.934)	(1.033)
	1980-present	Hartland	Canberra
		(1.059)	(1.084)

  6, 87.9, 64.2, 105.4, 78.1, 119.2, 151.8, 201.3, 110.6, 164.5, 158.5, 120.8) with a mean solar activity of about 119.7. The estimations for β and P 500 are: β = 0.0068 with the 95% C.I. [0.0055; 0.0081], P 500 = 2.49% with the 95% C.I. [1.65; 3.33].

Table 6

 6 ).

		Extreme storms High storms
	T h	19	538
	T l	4	164

Table 6 :

 6 Number of extreme and high storms for low and high intensity periods, respectively.

400 is the value of the ap index characterizing extreme storms. More details about high and extreme levels are given in Sections 3.2 and 3.3

Acknowledgements

The authors wish to thank the editor, the associate editor and two referees of Statistics and Its Interface (SII) for providing helpful comments and suggestions. This PhD candidate work at IMT lab of UPS is founded by CNES and TASF.

Appendix A Chi-square test

The utilization of the parameter P 400 (the probability for a high level storm to grow into an extreme level storm) for the extrapolation supposes that the


Marco Pizzolato
email: marco.pizzolato@inria.fr

Demian Wassermann
email: demian.wassermann@inria.fr

Timothé Boutelier
email: timothe.boutelier@olea-medical.com

Rachid Deriche
email: rachid.deriche@inria.fr

Exploiting the Phase in Diffusion MRI for Microstructure Recovery: Towards Axonal Tortuosity via Asymmetric Diffusion Processes

Keywords: diffusion MRI, phase, EAP, asymmetry, tortuosity, axon

Microstructure recovery procedures via Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) usually discard the signal's phase, assuming symmetry in the underlying diffusion process. We propose to recover the Ensemble Average Propagator (EAP) directly from the complex DW signal in order to describe also eventual diffusional asymmetry, thus obtaining an asymmetric EAP. The asymmetry of the EAP is then related to tortuosity of undulated white matter axons, which are found in pathological scenarios associated with axonal elongation or compression. We derive a model of the EAP for this geometry and quantify its asymmetry. Results show that the EAP obtained when accounting for the DW signal's phase provides useful microstructural information in such pathological scenarios. Furthermore, we validate these results insilico through 3D Monte-Carlo simulations of white matter tissue that has experienced different degrees of elongation/compression.

Introduction

A rather unexplored field in Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) is the use of the phase component of the signal to quantify white matter microstructural characteristics. Despite the DW signal's complex nature, current approaches discard the phase and characterize microstructure using just the signal's magnitude [START_REF] Assaf | AxCaliber: a method for measuring axon diameter distribution from diffusion MRI[END_REF][START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF][START_REF] Alexander | A general framework for experiment design in diffusion MRI and its application in measuring direct tissue-microstructure features[END_REF]. In this paper, we show that the signal's phase component plays an important role in microstructure characterization. It provides a novel approach towards the quantification of diffusion properties. Particularly, we exploit the signal's phase to quantify different levels of tortuosity in undulated axons within a voxel, which can be associated with pathology through, for instance, axonal compression.

Marco Pizzolato and Demian Wassermann contributed equally to this work. The author expresses his thanks to Olea Medical and the Provence-Alpes-Côte d'Azur (PACA) Regional Council for providing grant and support for this work. Fig. 1. From left to right: elongated and compressed WM tissue, reproduced with permission from [START_REF] Shacklock | Biomechanics of the Nervous System: Breig Revisited[END_REF]; straight/undulated axonal models for different tortuosity rates α.

White matter (WM) tissue may present different microstructural organizations depending on its nature and pathological condition. For instance, axons composing the tissue can be straight as shown in fig. 1 left. However axons may appear sinusoidal, and pathological conditions such as those caused by cervical cord injury can result in changes of undulation amplitude and tortuosity due to longitudinal elongation or compression [START_REF] Nilsson | The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study[END_REF][START_REF] Shacklock | Biomechanics of the Nervous System: Breig Revisited[END_REF], as shown in central fig. 1. The movement of spins within the WM tissue is influenced by the type of organization. Thus, a profound knowledge of the diffusion process can give insights into its microstructural properties.

DW-MRI allows to probe diffusion in a tissue sample by measuring a complexvalued signal related to the displacement probability of spins. The average displacement probability of the ensemble of spins contained within a voxel is known as Ensemble Average Propagator (EAP), and is related to the signal attenuation via an inverse Fourier transform [START_REF] Tanner | Restricted Self-Diffusion of Protons in Colloidal Systems by the Pulsed-Gradient, Spin-Echo Method[END_REF]. Tissue microstructure influences spins displacements. Hence, EAP and DW signal can be exploited to measure microstructural properties. Recent work [START_REF] Nilsson | The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study[END_REF] proves that undulated axons, modeled as sinusoids, are a confound on axonal diameter estimation from DW images.

We show that by exploiting the DW signal's phase it is possible to characterize axonal tortuosity, which can be linked to elongation/compression phenomena in pathological scenarios. Compressed axons are partially convoluted and show irregular undulation with a non-uniform tortuosity along the longitudinal direction (see fig. 1 right). Hence, we model the compressed axon as a sinusoid with varying tortuosity, inducing asymmetry in its shape. Moreover, as the degree of compression increases, the tortuosity variation rate along the axon also increases. Spins diffusing within these axons undergo an asymmetric diffusion process resulting in an informative DW signal's phase, sensible to increasing tortuosity rates. The amount of diffusion asymmetry thus gives important insights into the underlying degree of compression.

Diffusion process asymmetries have already been investigated [START_REF] Özarslan | Remarks on q-space MR propagator in partially restricted, axially-symmetric, and isotropic environments[END_REF]. However, current approaches to obtain the EAP only use the DW magnitude signal and assume axial symmetry along each diffusion direction. We show that, by exploiting the DW signal's phase component, we recover an EAP describing also asymmetric diffusion which asymmetry is linked to the axonal tortuosity rate, and thus compression. To prove this, we first derive an analytic model of the EAP corresponding to the proposed axonal geometry, and obtain the correspondent DW complex signal. Using the complex signal we recompute the EAPs for increasing tortuosity rates, corresponding to scenarios of increasing axonal compression. We propose to quantify asymmetry as the Hellinger distance [START_REF] Hellinger | Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen[END_REF] between each EAP and its axially reflected version. We finally present experiments elucidating the link between EAP asymmetry and axonal tortuosity rates, thus compression degrees. Results are shown for EAPs generated analytically in 2D and validated by means of 3D Monte-Carlo simulations performed with Camino [START_REF] Hall | Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI[END_REF].

Theory and Methods

In this section, we introduce a model describing the undulated axon geometry with varying tortuosity. We derive the corresponding equations for the 2D Ensemble Average Propagator (EAP) and discuss its relation with the magnitude and phase of the DW-MRI signal.

Tortuous Axon Model

To analyze the DW signal in axonal compression cases, we model straight axons, in agreement with the usual hypothesis in the area, and undulated axons showing tortuosity. Our model for undulated axons is an extension of that proposed in [START_REF] Nilsson | The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study[END_REF]. We consider a sinusoidal axon that lies in the xz-plane and evolves along z with oscillations along x. Our axon model has infinitesimal thickness and goes through two full periods within the interval [-Z, Z] with Z > 0:

x = f (z) = A sin 2π(z -Z) L(z, α) (1)
where A and L(z, α) are respectively the amplitude and the wavelength of the sinusoid, with 0 < A ≤ X and L(z, α) > 0. We model the varying wavelength as

L(z, α) = α(z + Z) + l, f or α > 0 ( 2
)
where α is the spatial rate regulating the wavelength increment along the zaxis. The wavelength increases progressively along the z-axis starting from a minimum value l in z = -Z. The spatial rate α is related to the degree of intrinsic asymmetry induced in the undulated shape along the z-axis. We consider the case where α = 0 as the straight axon, that is when A = 0, as shown in fig. 1 right. To derive the EAP in the following section it is useful to define the curve length for eq. ( 1). In particular the length between z and z + ∆z is given by

s(∆z; z) = z+∆z z 1 + A cos 2π(ζ -Z) L(ζ, α) 2π(l + 2αZ) L(ζ, α) 2 2 dζ (3)
where ∆z is a displacement along the z-axis with respect to the coordinate z.

Ensemble Average Propagator for Tortuous Axons

To obtain the complex DW signal in the case of tortuous axons model in section 2.1, we derive the correspondent Ensemble Average Propagator (EAP). The EAP for a spin displacement (∆x, ∆z) subject to a diffusion time

t d is [6] EAP(∆x, ∆z; t d ) = ∞ -∞ ρ 0 (x, z)P (∆x, ∆z|x, z, t d ) dxdz = ∞ -∞ ρ 0 (z)P (∆x, ∆z|z, t d ) dz (4)
where ρ 0 is the initial spin density and P is the displacement pdf. The second equality in eq. ( 4) is due to our model parametrization in eq. ( 1). We compute ρ 0 by assuming an initial uniform spin density along our axon model

ρ 0 (z) = cos (arctan [f (z)]) -1 (5)
where f (z) is the derivative of f in eq. ( 1). To produce an approximation for the displacement probability P in eq. ( 4), we decompose it using:

P (∆x, ∆z|z, t d ) = P (∆x|∆z, z)P (∆z|z, t d ). (6)
Using this decomposition, we assume that spins diffuse freely along the curve, i.e. they follow a normal distribution with variance 2Dt d where D is the selfdiffusion coefficient of water. Then, the probability of a spin experiencing a net displacement along the z-axis ∆z can be approximated by [START_REF] Nilsson | The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study[END_REF] P (∆z|z,

t d ) = 1 √ 4πDt d e - s(∆z;z) 2 4Dt d (7)
where s, the arc length of the axon between z and z + ∆ z , is in eq. ( 3). As the coordinate x on the axon is univocally determined by z, as specified in eq. ( 1), the displacement density along x provided that we know ∆z and z is [4]

P (∆x|∆z, z) = δ D (f (z) + ∆x -f (z + ∆z)) (8)
where δ D is the Dirac delta function. Given eqs. ( 5) to [START_REF] Hellinger | Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen[END_REF], the EAP of a tortuous axon, according to our model in section 2.1, can be obtained by simple numerical integration of eq. ( 4). We proceed to express the relationship between EAP and signal, then we quintify EAP asymmetry by exploiting the phase information.

Diffusion Signal and EAP Asymmetry

The signal acquired in DW-MRI can be represented as function of the diffusion time t d and frequency vector q = γδG/2π where γ is the gyromagnetic ratio, δ is the diffusion pulse duration and G is the diffusion gradient vector. The measured q-dependent signal S(q) has a complex nature and its attenuation E(q) = S(q)/S(0) is related to the EAP via a Fourier relationship [START_REF] Tanner | Restricted Self-Diffusion of Protons in Colloidal Systems by the Pulsed-Gradient, Spin-Echo Method[END_REF]:

E(q x , q z ; t d ) = ∞ -∞ EAP(∆x, ∆z|t d )e -j2π(qx∆x+qz∆z) d∆z d∆x (9)
which is valid under narrow pulse approximation and where q x and q z are the components of the frequency vector q. The EAP usually calculated from the magnitude attenuation |E| via inverse Fourier transform is symmetric due to the zero phase assumption [START_REF] Bracewell | The Fourier transform and its applications[END_REF]. We consider instead the complex attenuation and characterize the asymmetry of the diffusion process by measuring the distance between the corresponding EAP and its axially reflected version. Despite different measures are possible, after defining the displacement vector r = (∆x, ∆z), we adopt the Hellinger distance [START_REF] Hellinger | Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen[END_REF] because it constitutes a proper metric between probability measures

H 2 = 1 2 EAP(r|t d ) -EAP(-r|t d ) 2 dr ( 10
)
where 0 ≤ H ≤ 1, 0 corresponding to equality and 1 to maximum inequality.

Monte-Carlo Simulation of DW-MRI in Tortuous Axons

To validate analytic results, we generate the DW complex signal by performing 3D Monte-Carlo simulations using an in-house customized version of Camino [START_REF] Hall | Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI[END_REF] that allows to obtain both Magnitude and Phase. This is done by considering both real and imaginary parts of the signal, thus replacing eq. ( 7) in [START_REF] Hall | Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI[END_REF] with

S(G, ∆, δ) = S(0, ∆, δ) ∞ -∞ P (φ) (cos(φ) + j sin(φ)) dφ ( 11
)
where φ is the phase accumulation, P (φ) the spin phase distribution and ∆ the time separation between diffusion gradients.

Experiments

In this section we present experiments showing the relationship between EAP asymmetry and axonal tortuosity for the EAP obtained as described in section 2.2 and with 3D Monte-Carlo simulations.

In the first case we compute the EAP of straight and undulated axons via eq. ( 4) for increasing tortuosity rates α ∈ [START_REF] Assaf | AxCaliber: a method for measuring axon diameter distribution from diffusion MRI[END_REF][START_REF] Hellinger | Neue Begründung der Theorie quadratischer Formen von unendlichvielen Veränderlichen[END_REF]. Axons present undulation amplitude A = 4 µm with a basal wavelength l = 50 µm [START_REF] Nilsson | The importance of axonal undulation in diffusion MR measurements: a Monte Carlo simulation study[END_REF]. We set D = 2×10 -9 m 2 s -1 [START_REF] Hall | Convergence and parameter choice for Monte-Carlo simulations of diffusion MRI[END_REF] and t d = 28.6 ms according to a realistic scenario. We fix the maximum observable displacement in both directions to X = Z = 50 µm and we define the observation frame as the xz-plane limited to the interval [-X, X] along x and [-Z, Z] along z. As α increases along the z-axis, eq. ( 2), the undulated axonal trajectory becomes compressed towards the frame's negative bound

α =0 α =1 α =2 α =3 α =4 α =5 α =6 α =7
Usual EAP from magnitude signal α =8 -Z (see fig. 1 right). The observation frame is discretized over a 257 × 257 displacement grid, and the EAP is calculated for each sample via numerical evaluation of eq. ( 4). The complex signal attenuation is obtained, for each considered α, as the Fourier transform of the EAP over the displacement grid according to eq. ( 9). The EAPs are then calculated as the inverse Fourier transform of both the complex and magnitude signal attenuations. Images in fig. 2 show the log-transformed EAPs calculated for the straight axon (α = 0) and for different degrees of axonal compression (α > 0). The first row contains the EAPs as usually calculated from the magnitude signal. The second row shows the EAPs obtained from the complex signal. We note that no differences exist for the straight axon α = 0. Despite the shape of the magnitude-derived EAP changes as α increases, we note that it always shows axial symmetry. On the other end, evident asymmetry is present in the case of the complex-derived EAPs for α > 0. When looking at the EAPs in the second row of fig. 2, the difference between the upper and lower lobes of each EAP clearly appears. We further note that the degree of asymmetry between the lobes is more marked at high α, for instance the difference between the lower and upper lobes for α = 8 is more accentuated than in the case of α = 1.

We quantify asymmetry in the EAPs by computing the proposed Hellinger distance, eq. ( 10), as function of α (fig. 3 left). We note the incremental trend between the complex-derived EAP asymmetry and the tortuosity rate, in agreement with images in the second row of fig. 2. To validate the results for the more realistic case in which axons have finite diameter, we reproduce experiments using 3D Monte-Carlo simulations. Complex signals are obtained via PGSE sequences with ∆ = 24.3 ms, δ = 12.9 ms and maximum gradient strength |G max | = 0.8 T m -1 . Gradient directions are selected over the undersampled

Discussion and Conclusion

The overall objective of this work is to assess the use of the DW signal's phase as microstructural marker. We show that the phase becomes informative whenever the underlying diffusion process manifests asymmetry along specific directions as a consequence of, for instance, a variation in the axonal tortuosity.

EAP asymmetry increases with larger tortuosity rates in our 2D model. This is observable in fig. 3 left and fig. 2 (second row), where values along radial lines passing through the origin are asymmetric for α > 0. The explanation of this phenomenon lies in the spin displacements asymmetry along the axes. Indeed the maximum distance a spin can travel along an axon evolving along the z-axis is defined by the diffusion coefficient and time. This distance is fixed irrespective of the fact that the spin displaces, along the axon, towards the positive or the negative values of the z-axis. This is the case of the straight axon (α = 0) where diffusion occurs equally in each direction (fig. 2). With undulation the maximum displacement along z, resulting from a movement along the axon, depends on the movement direction of the spin, and precisely on the average tortuosity the axon shows in each direction: the ratio between the length of the trajectory and the corresponding length along the z-axis. If the tortuosity in the two directions is different (α > 0), then the displacement probability of the spin along the z-axis, thus the x-axis, is asymmetric with respect to the origin. Since the EAP accounts for all of the contributions of the spin displacement probabilities along the axon, the EAP is symmetric when axonal tortuosity is constant (α = 0), and is asymmetric otherwise (α > 0). EAP asymmetry is sensible to axonal tortuosity in realistic scenarios. In 3D experiments spins diffuse transversely to the axon as well as along it. Hence, the trajectory a spin performs is convoluted and leads to smaller displacements along the axon, causing the asymmetry to have a milder effect on the overall EAP. This directly implies that a good characterization of the EAP asymmetry may be better achieved at very high |q| values as shown in fig. 3 central and right. Lower gradient strengths, hence displacement resolutions, affect the measured asymmetry at very short wavelengths as for the case of α = 4 and |G max | = 0.1 T m -1 . However, our proposed measure of asymmetry using the complex signal combined with the Hellinger distance, eq. (10), as microstructure marker for axonal tortuosity/compression, remains effective. As a future work we suggest to explore higher diffusion times, and thus ∆. While a very high |q| is hard to achieve, longer diffusion times allow spins to diffuse more thus sensing the asymmetry at greater displacements, which are observable also at lower |q|.

Real data acquisition for this scenario is challenging, specially due to technical difficulties such as the extreme sensibility of phase to noise. However, although other techniques assume straight axons, simulations show that undulation combined with tortuosity significantly affect the DW signal and can be detected through diffusion asymmetry by exploiting the phase. We strongly believe that this pioneering work opens new roads and perspectives towards measuring WM tissue properties via DW-MRI.

Fig. 2 .

 2 Fig. 2. Log-scaled EAPs obtained for straight and undulated axons with different tortuosity rates α. The first row shows EAPs obtained via inverse Fourier transform of the magnitude signal attenuation: EAPs show axial symmetry with respect to the axis passing through the origin. The second row shows EAPs obtained from the complex signal: axial asymmetry is present and increases with the tortuosity rate α.

Fig. 3 .

 3 Fig. 3. Hellinger distances as function of tortuosity rate α obtained for the analytic model (left) and via Monte-Carlo simulation (central, right). Right figure shows the behaviour for state-of-the-art human MRI scanners such as the HCP one.