B Cooper Boniece

Herwig Wendt

Gustavo Didier

Patrice Abry

WAVELET-BASED DETECTION AND ESTIMATION OF FRACTIONAL L ÉVY SIGNALS IN HIGH DIMENSIONS

Keywords: multivariate self-similarity, operator fractional Lévy motion, wavelets, random matrices, high dimensions

In the modern world, systems are routinely monitored by multiple sensors, generating "Big Data" in the form of a large collection of time series. However, dynamic signals are often low-dimensional and characterized by joint scale-free dynamics (self-similarity) and non-Gaussianity. In this paper, we put forward a statistical methodology for identifying the number of multivariate self-similar, Lévydriven components immersed in high-dimensional noise, as well as for estimating the underlying scaling exponents. It relies on the analysis of the evolution over scales of the eigenvalues of random wavelet matrices. Monte Carlo simulations show that the proposed methodology is accurate for realistic sample sizes. This holds even at low signal-to-noise ratios and for a large number of observed mixed and noisy time series. The mathematical framework further allows us to analyze the impact of the tails of the Lévy noise marginal distribution on the estimation performance.

INTRODUCTION

Context: many sensors, few sources. In the modern world of "Big Data," a plethora of sensors monitor natural and artificial systems, generating large data sets in the form of several joint time series. This can be seen in many fields of application. In neuroscience, for example, the number of macroscopic brain activity time series ranges from hundreds (MEG data) to several tens of thousands (fMRI data) [START_REF] Ciuciu | Scale-free and multifractal time dynamics of fMRI signals during rest and task[END_REF]. Likewise, in climate studies, dealing with large numbers of measured components has become standard [START_REF] Isotta | The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data[END_REF]. Nevertheless, the number of physical mechanisms, or sources, driving the spatio-temporal dynamics of a system is usually far smaller than the number of recorded time series (sparsity in high dimensions). This has created a strong demand for data analytic methods that can statistically identify sources amidst noisy measurements and extract their informational content [START_REF] Wainwright | High-Dimensional Statistics: a Non-Asymptotic Viewpoint[END_REF].

(Covariance) self-similarity and non-Gaussian behavior.

A signal is called scale invariant when its temporal dynamics lack a characteristic scale. In the presence of scale invariance, the analyst's main interest is in identifying mechanisms, often in the form of scaling exponents, that relate the continuum of scales [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. Scale Work supported by Grant ANR-16-CE33-0020 MultiFracs. G.D. was partially supported by the ARO grant W911NF-14-1-0475. G.D.'s long term visits to France were supported by ENS de Lyon, CNRS and the Carol Lavin Bernick faculty grant.

invariance often manifests itself in physics and engineering in the form of self-similarity. A signal X is called self-similar when its finite-dimensional distributions (fdd) are invariant under suitable time scaling, i.e., {X(t)} t∈R fdd = {a H X(t/a)} t∈R , a > 0, where 0 < H < 1 is called the Hurst parameter. A cornerstone model of scale invariance, fractional Brownian motion (fBm) is the only Gaussian, self-similar process with stationary increments [START_REF] Taqqu | Fractional Brownian motion and long range dependence[END_REF]. Estimation of the Hurst parameter H plays a key role in signal processing tasks such as characterization, diagnosis, classification and detection. The wavelet transform underpins efficient multiscale representations of self-similar or fractional systems, and it provides the analytical basis for well-established methodologies for the estimation of H [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF].

On the other hand, non-Gaussian behavior is pervasive in a myriad of natural and artificial systems displaying scale invariance (e.g., turbulence [START_REF] Boniece | On fractional Lévy processes: tempering, sample path properties and stochastic integration[END_REF] or Internet traffic [START_REF] Willinger | Scaling phenomena in the Internet: critically examining criticality[END_REF]). Due to its breadth and flexibility, the family of fractional Lévy processes [START_REF] Benassi | Identification and properties of real harmonizable fractional Lévy motions[END_REF][START_REF] Marquardt | Fractional Lévy processes with an application to long memory moving average processes[END_REF][START_REF] Lacaux | Hurst exponent estimation of fractional Lévy motion[END_REF] has become popular in physical applications involving non-Gaussianity [START_REF] Zhang | A least squares estimator for Lévy-driven moving averages based on discrete time observations[END_REF][START_REF] Xu | The switch in a genetic toggle system with Lévy noise[END_REF]. Nevertheless, while of great importance in practice, the modeling of multivariate fractional Lévy signals is a research topic that has been relatively little explored [START_REF] Barndorff-Nielsen | Multivariate supOU processes[END_REF][START_REF] Moser | Functional regular variation of Lévy-driven multivariate mixed moving average processes[END_REF]. In [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF], operator fractional Lévy motion (ofLm), a generalization of the univariate fBm and its multivariate counterpart, operator fBm (ofBm) [START_REF] Mason | Sample path properties of operatorself-similar Gaussian random fields[END_REF][START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF], is put forward as a model for non-Gaussian multivariate (covariance) self-similarity. In particular, the ofLm family comprises instances where there are multiple correlated fractional Lévy coordinate processes with possibly distinct self-similarity exponents Hm, m = 1, . . . , M , occurring in a non-canonical set of coordinates (mixing).

Related work.

In the statistical signal processing literature, the problem of statistically identifying the number and properties of sources in multivariate (fixed dimensions) or high-dimensional noisy signals has been studied for at least three decades now [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF][START_REF] Stoica | Model-order selection: a review of information criterion rules[END_REF][START_REF] Liavas | On the behavior of information theoretic criteria for model order selection[END_REF][START_REF] Da Costa | Comparison of model order selection techniques for highresolution parameter estimation algorithms[END_REF]. Examples of the proposed techniques include principal component analysis, factor analysis and sparse graphical Gaussian models [START_REF] Johnson | Applied Multivariate Statistical Analysis[END_REF]. Nevertheless, there has been a paucity of estimation methodologies for both high-dimensional and scale invariant signals [START_REF] Zhang | Identifying cointegration by eigenanalysis[END_REF]. To address this issue, the so-named wavelet eigenanalysis estimation methodology was recently constructed [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF][START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. It was shown to lead to efficient and robust estimation of self-similarity exponents in both Gaussian multivariate/high-dimensional [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF][START_REF] Abry | Detecting and estimating multivariate self-similar sources in high-dimensional mixtures[END_REF] and non-Gaussian multivariate [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF] settings.

Goal, contributions and outline.

In this paper, we tackle the problem of high-dimensional inference for non-Gaussian covariance self-similar signals. The observations are assumed to appear in the form of an ofLm plus high-dimensional noise. The concept of multi-variate (covariance) self-similarity is briefly laid out at the beginning of Section 2. The definition and properties of ofLm are recapped in Section 2.1. The recently proposed wavelet eigenanalysis methodology is detailed in Section 2.2. The high-dimensional model and model order selection procedure are described in Section 3. Performance with respect to non-Gaussianity and increasing dimension of both model order selection and Hurst exponent estimation are assessed by means of Monte Carlo studies and discussed in Section 4.

OPERATOR FRACTIONAL L ÉVY MOTION AND WAVELET EIGENANALYSIS-BASED ESTIMATION

Covariance operator self-similarity. Let H = P diag(H)P -1 be the so-named Hurst matrix parameter, where the vector H contains the Hurst exponents (eigenvalues). We say a stochastic process X = (X1, . . . , XM ) is covariance (multivariate) operator selfsimilar (c.o.s.s.) with Hurst matrix parameter H when

EX(s)X(t) * = a H EX(s/a)X(t/a) * a H * (1)
for all a > 0, where a H := +∞ k=0 log k (a)H k /k!. When the mixing matrix P is diagonal, namely, when we can set P ≡ I, the covariance self-similarity relation (1) takes the simple form of componentwise covariance self-similarity relations

EX (s)X (t) = a H +H EX (s/a)X (t/a) (2)
for a > 0 and , = 1, 2, . . . , M .

Operator fractional Lévy motion

OfLm is a class of non-Gaussian multivariate fractional processes displaying the same covariance structure of ofBm [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF]. In particular, ofLm is c.o.s.s. For clarity of exposition and without loss of generality, we construct ofLm in the bivariate case. Let {L(t) = (L1(t), L2(t))} t∈R be a two-sided symmetric Lévy process in R 2 with EL(1)L(1) * =: ΣL, |ΣL| < ∞ (e.g., [START_REF] Sato | Lévy Processes and Infinitely Divisible Distributions[END_REF]). Let 0 < H1 ≤ H2 < 1, and define the (pre-mixed) process X by the componentwise convolution X (t) = (gt * L)(t), where g is the diagonalmatrix-valued fractional kernel gt(u)

:= u D + -(u -t) D + , D = diag(H) -1 2 I.
Then, the entrywise processes XH ( = X), = 1, 2, are generally correlated fractional Lévy processes with corresponding Hurst parameters H ∈ (0, 1) (e.g., [START_REF] Benassi | Identification and properties of real harmonizable fractional Lévy motions[END_REF][START_REF] Marquardt | Fractional Lévy processes with an application to long memory moving average processes[END_REF]). In particular, each has stationary increments and its covariance function is identical to that of fBm, i.e.,

EXH (t)XH (s) = {|t| 2H + |s| 2H -|t -s| 2H }σ 2 /2. (3)
Note that, if non-Gaussian, the stochastic behavior of the process XH is not characterized by its covariance function. Let P be a 2 × 2, real-valued, invertible matrix. We define the (bivariate) ofLm X H,L,P by

{X H,L,P 1 (t), X H,L,P 2 (t)} t∈R = P {XH 1 (t), XH 2 (t)} t∈R
(in short, X H,L,P = P X ). Note that the parameters H and ρ := Corr(XH 1 (1), XH 2 (1)) cannot be selected independently (see [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF] for more details).

Wavelet-based scaling exponent estimation

Multivariate wavelet transform. Let ψ0 be a mother wavelet, i.e., an oscillating function that is well localized in both time and frequency [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. A mother wavelet is parametrized by the so-named number of vanishing moments N ψ , i.e., a positive integer such that,

for κ = 0, . . . , N ψ -1, R t κ ψ0(t)dt ≡ 0 and R t N ψ ψ0(t)dt = 0. Let {ψ j,k (t) = 2 -j/2 ψ0(2 -j t -k)} (j,k)∈Z 2 be
the collection of dilated and translated templates of ψ0 that forms an orthonormal basis of L 2 (R). For the M -variate process X = {X(t)} t∈R and its components, let DX m (2 j , k)

:= 2 -j/2 ψ(2 -j t -k)|Xm(t) , m ∈ {1, . . . , M }. The multivariate DWT of X is defined as DX (2 j , k) := (DX 1 (2 j , k), . . . , DX M (2 j , k)) * , ∀k ∈ Z, ∀j ∈ Z.
Covariance operator self-similarity in the wavelet domain. It can be shown that the wavelet coefficients {DX (2 j , k)} k∈Z satisfy the covariance (operator) self-similarity relation [START_REF] Boniece | On Scale Invariance and Wavelet Analysis: Transience, Operator Fractional Lévy Motion, and High-Dimensional Inference[END_REF]

EDX (2 j , k)DX (2 j , k) * = 2 j(H+ 1 2 I) EDX (1, k)DX (1, k) * 2 j(H * + 1 2 I) , (4)
for every fixed octave j and any shift k. When P = I (no mixing), similar to (2), M -variate covariance (operator) self-similarity in the wavelet domain boils down to the entrywise covariance selfsimilarity relations, for , = 1, . . . , M ,

EDX (2 j , k)DX (2 j , k) * = 2 j(H +H +1) EDX (1, k)DX (1, k) * . (5)
Estimation of H. The sample wavelet spectrum (variance) is given by the M × M random matrices

SX (2 j ) = 1 nj n j k=1 DX (2 j , k)DX (2 j , k) * , j = j1, . . . , j2, (6)
where n is the time series sample size and nj = n/2 j . When there is no mixing (P = I), i.e., when relation (5) holds, the scaling exponents may be estimated by means of entrywise multiscale logregressions based on the system (6) (see [START_REF] Amblard | Identification of the multivariate fractional Brownian motion[END_REF]). Otherwise, i.e., starting from (4) (mixing, or non-diagonal P ), this procedure does not yield meaningful results (see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF][START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] for further discussions). An estimation methodology for general (non-diagonal) P can be constructed as follows. Consider the eigenvalues ΛX (2 j ) = (λ1(2 j ), . . . , λM (2 j )) * of the random matrix SX (2 j ). The estimators H = ( H1, . . . , HM ) for (H1, . . . , HM ) are defined by means of the multiscale weighted linear regressions

Hm = j 2 j=j 1 wj log 2 λm(2 j ) 2 - 1 2 , ∀m. (7)
It was shown mathematically or based on Monte Carlo experiments that H has good performance over large and finite samples in multivariate/high-dimensional Gaussian frameworks [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF][START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF][START_REF] Abry | Detecting and estimating multivariate self-similar sources in high-dimensional mixtures[END_REF], as well as in multivariate non-Gaussian (Lévy-driven) frameworks [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF].

ESTIMATION IN HIGH DIMENSIONS

High-dimensional model

Let W be a full rank L × M matrix, where L M . Also let X = {X(t)} t∈R be a M -variate ofLm with Hurst matrix H = P diag(H)P -1 , as described in Section 2. We assume the L observed time series follow the model

Y (t) := W X(t) + σN N (t). (8)
In [START_REF] Willinger | Scaling phenomena in the Internet: critically examining criticality[END_REF], the noise term N (t) = N1(t), . . . , NL(t) * t∈R is made up of L independent vectors, each consisting of i.i.d. standardized Gaussian variables, and σ 2 N is a control parameter for the noise variance. Assume, without loss of generality, that W P has rows with unit norm, and further suppose σX := σ1 = . . . = σM in (3). The Signal-to-Noise Ratio (SNR) is defined as

SNR := σX σN . ( 9
)
This definition is natural for scaling processes because, for a given value of the SNR, it sets the size of eigenvalues at fine scales j → 0, regardless of the values of M and L; see Fig. 1. Note that it is much more penalizing than the conventional definition of SNR because [START_REF] Benassi | Identification and properties of real harmonizable fractional Lévy motions[END_REF] implies that, for a fixed SNR value, the overall energy of noise in the observations increases linearly with L.

We assume the Lévy prrocess driving the ofLm has a symmetric exponentially tempered α-stable distribution (e.g. [START_REF] Grabchak | Tempered Stable Distributions: Stochastic Models for Multiscale Processes[END_REF]). Such distributions have moments of all orders and are parametrized by a stability index α ∈ (0, 2) and a tempering parameter γ > 0. A smaller tempering parameter γ > 0 corresponds to heavier-tailed behavior. Furthermore, it can be shown [START_REF] Boniece | On Scale Invariance and Wavelet Analysis: Transience, Operator Fractional Lévy Motion, and High-Dimensional Inference[END_REF] that at fine timescales the corresponding ofLm is approximately operator stable (heavy-tailed), whereas on coarse timescales it is close to an ofBm (i.e., Gaussian), with transition controlled by γ. Accordingly, for this class of ofLm, we treat γ as a proxy for measuring departures from Gaussianity.

Model order selection procedure

The proposed procedure relies on the eigenstructure of the random matrices SY (2 j ) over several octaves j, i.e., on the sample wavelet spectra of the L-variate observations Y . In the absence of noise (SNR = +∞), SY (2 j ) has rank M for all scales 2 j (as long as nj > M ). Conversely, when noise is added (SNR < +∞), SY (2 j ) has rank L. Hence, assuming the SNR is sufficiently large, heuristically speaking each random matrix SY (j) possesses L -M small eigenvalues corresponding to noise and M large eigenvalues corresponding to the hidden ofLm components. The former set of eigenvalues manifests itself in the form of log-eigenvalue curves with negative slopes, whereas the latter leads to log-eigenvalue curves with positive slopes driven by the exponents Hm > 0. This is illustrated in Fig. 1, which shows averages and error bars for log-eigenvalues generated based on 1000 independent realizations of ofLm (see Sec. 4 for details on numerical simulations). So, let ΛY (2 j ) ∈ R L be the vector of eigenvalues of SY (2 j ). We propose to estimate the number of components M as the number of estimated slopes that are found to be positive, where each slope estimate is obtained as in [START_REF] Boniece | On fractional Lévy processes: tempering, sample path properties and stochastic integration[END_REF] by log-linear regressions on the entries of ΛY (2 j ).

Once the number of ofLm components M has been estimated, a vector of estimated scaling exponents ( H1, . . . , H M ) stems naturally from the procedure given by ( 7), applied to the M largest values of ΛY (2 j ).

PERFORMANCE

Monte Carlo simulation.

We evaluate the performance of the proposed model order detection and high-dimensional estimation procedure for 1000 independent realizations of size n = 2 14 of ofLm with M = 4 components and exponents H = (0.4, 0.6, 0.75, 0.9). We study different degrees of deviation from Gaussian behavior (tempering parameter γ = (10 -7 , 10 -6 , . . . , 10 -1 )) with α = 1.1 fixed and L = (10, 20, 50, 100) observed time series, for various SNR values. The mixing matrix is set to

P = ((1 1 1 -1) T , (1 -1 1 1) T , (1 1 1 1) T , (1 -1 -1 1) T ).
The matrix W is picked randomly and is normalized such that the rows of W P have unit norm. Estimation of the scaling exponents Hm, m = 1, . . . , M , is performed as in ( 7) by means of weighted linear regressions for j ∈ (4, j (m) max ), where j (m) max is the largest scale for which λm(j) > 0, ∀m. The ofLm paths were generated via an adaptation of the FFT-based algorithm described in [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF]. In the notation of [START_REF] Stoev | Simulation methods for linear fractional stable motion and FARIMA using the Fast Fourier Transform[END_REF], the mesh size was chosen m = 16 with kernel cutoff parameter M = 2 14 . The symmetric tempered stable processes driving the ofLm paths were simulated via an accept-reject procedure (see [START_REF] Baeumer | Tempered stable Lévy motion and transient super-diffusion[END_REF]). Model order selection performance.

For different values of γ (namely, levels of departure from Gaussianity) and of the SNR, and for different dimensions L, Fig. 2 plots histograms of the estimated model orders M generated by means of the proposed model order selection procedure. The results demonstrate that for γ = 10 -1 , 10 -3 (weak departure from Gaussianity), the procedure selects the correct model order M = 4 with high probability as long as SNR ≥ 0dB.

Note that the procedure's good performance is robust with respect to the value of L. Strikingly, we also observe a phenomenon that may be called the blessing of dimensionality: for lower SNR, larger dimensional observations turn out to be beneficial. For example, for SNR= -6dB, the selected model orders are correct with probability close to one for L ≥ 50, yet the correct model order is only detected one out of two times for L = 10. This phenomenon has been observed in other frameworks (cf. [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]), in particular in Gaussian systems [START_REF] Abry | Detecting and estimating multivariate self-similar sources in high-dimensional mixtures[END_REF], and may be related to the fecund eigenvalue repulsion effect [START_REF] Tao | Topics in Random Matrix Theory[END_REF]. The fact that model order selection becomes easier in high dimensions can also be seen in Fig. 1, where the separation between the curves with negative and positive slopes becomes sharper as the dimension L increases. Turning back to Fig. 2, as expected, the stronger the departure from Gaussianity, the more difficult the model order selection problem. Nevertheless, the model order selection procedure remains quite reliable. For L ≥ 20 and SNR=0dB, the probability of correct detection is about 80% for γ = 10 -5 , 10 -7 (instead of ≈100% in lighther-tailed instances). For low SNR values, it is generally true that the majority of detection attempts are often too conservative ( M < M ). Note that model order selection performance could be optimized by searching the best threshold value, which in our simulations is set a priori to 0.

Estimation performance in high dimensions. For instances involving low noise levels (30dB SNR), Fig. 3 plots bias, standard deviation (std) and root-mean-squared error (rms) for estimates Hm as a function of the dimension L = 10, 20, 50, 100 and γ. The plots reveal a number of properties. First, as expected, it can be seen that the estimation performance tends to deteriorate for stronger departures from Gaussianity (smaller γ). This is essentially due to a reduction in the standard deviation as γ increases, whereas, remarkably, the bias remains unaffected. These results corroborate similar findings for the bivariate, bi-dimensional case L = M = 2 without noise [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF]. Second, and interestingly, the results indicate that there is an inversion in the relative estimation difficulty between the smallest exponent H1 than the largest exponent H4 as γ increases. In other words, for small values for γ (heavier tails), the smallest exponent H1 has lower std and rms than the largest exponent H4; by contrast, the opposite tends to be observed for large values for γ. Finally, estimation performance generally improves for larger L, both in terms of bias and std. This resonates with the blessing of dimensionality phenomenon (improved model order selection performance for large L) observed in Figs. 1 and2.

Gaussianity of estimates Hm. In Fig. 4, we investigate the influence of L and γ on the finite-sample distribution of Hm. We report qq-plots for estimates Hm against the standard Normal distribution for instances with very little noise (30dB SNR). We can observe that, regardless of L, the finite-sample distribution of Hm becomes more and more Gaussian as the regularization value γ increases, as expected and reported in [START_REF] Boniece | On multivariate non-Gaussian scale invariance: fractional Lévy processes and wavelet estimation[END_REF] for M = L = 2. The results also show that there is no clear difference between the histograms for the smaller exponent H1 and for the largest exponent H4. Moreover, the dimension L seems to bear little influence on the distribution of Hm.

CONCLUSIONS

In this paper, we construct statistical methodology for the identification of the number of multivariate self-similar, Lévy-driven components immersed in high-dimensional noise, as well as for the estimation of the associated scaling parameters. The technique draws upon the analysis of the evolution over scales of the eigenvalues of random wavelet matrices (wavelet eigenanalysis). Monte Carlo simulations reveal the effect of high dimensions on the estimation procedure and show that the proposed methodology is accurate for realistic sample sizes. This holds even at relatively low signal-to-noise ratios and for a large number of observed mixed and noisy time series. By tuning the tail parameter of the Lévy noise distribution, we further quantify the impact of non-Gaussianity on the estimation performance, where strongly non-Gaussian instances are shown to be typically more statistically challenging.

Fig. 1 .

 1 Fig. 1. Log-eigenvalues of SY (2 j ). For γ = 10 -5 and SNR=0dB, plots of log 2 ΛY (2 j ) for M = 4 and L = 10, 100 (left and right plots, respectively).

Fig. 2 .

 2 Fig. 2. Model order selection performance. Histograms of selected model orders M for different SNR (increasing from left to right columns) and γ (decreasing departure from Gaussian from top to bottom row), respectively, and M = 4. The red/green/blue/grey bars correspond to L = 10/20/50/100, respectively.

Fig. 3 .

 3 Fig. 3. Estimation performance. Bias, std and rms (left to right column) for L = 10/20/50/100(top to bottom row) and fixed SNR=30dB as a function of γ.

Fig. 4 .

 4 Fig.[START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. Gaussianity of Hm. For SNR=30dB, γ = {10 -7 , 10 -3 } (first two and last two columns, respectively) and L = 10, 100, qqplots of the estimates H4 and H1 (largest and smallest exponents, top and bottom rows of subplots) against the standard Normal distribution.


Nicolas Riteau

Jean-François Bernaudin

Christine Huppertz

Benedikt Jäger

Grazyna Wieczorek

Peggy Engelhard

Stephen J Oliver

Franz-Georg Bauernfeind

Amanda Littlewood-Evans

Tobias Welte

In addition to mTOR and JAK/STAT, NLRP3 inflammasome is another key pathway activated in sarcoidosis

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In addition to mTOR and JAK/STAT, NLRP3 inflammasome is another key pathway activated in sarcoidosis Nicolas Riteau 1 and Jean-François Bernaudin 2,3,4 Affiliations: Nicolas.riteau@cnrs-orleans.fr Jean-francois.bernaudin@aphp.fr Sarcoidosis is a rare "old disease" first described over a century ago, known as a multifaceted systemic multiorgan inflammatory granulomatosis characterized by the presence of non-necrotizing epithelioid granulomas [START_REF] Valeyre | Sarcoidosis[END_REF][START_REF] Grunewald | Sarcoidosis. Nature reviews Disease primers[END_REF]. In many ways sarcoidosis remains a mysterious disease with numerous unsolved knowledge gaps. The still unknown etiology/etiologies of sarcoidosis is one of these gaps, i.e. the disease develops in genetically predisposed individuals following as-yet-unknown antigens exposure. Various triggers or causes have been suspected, including bacterial agents (such as Propionibacterium acnes or mycobacterial antigens) or environmental particulates (such as crystalline silica), but their true identity and mechanisms of action are still non-elucidated [START_REF] Valeyre | Sarcoidosis[END_REF][START_REF] Grunewald | Sarcoidosis. Nature reviews Disease primers[END_REF][START_REF] Jonsson | Silica dust and sarcoidosis in Swedish construction workers[END_REF]. The second main gap regards the nature of the granuloma foundation, the hallmark of the disease, as its mechanisms related to complex immunopathogenesis are only partially understood [START_REF] Valeyre | Sarcoidosis[END_REF][START_REF] Grunewald | Sarcoidosis. Nature reviews Disease primers[END_REF][START_REF] Jonsson | Silica dust and sarcoidosis in Swedish construction workers[END_REF].

Furthermore, if the highly variable course of the disease is well documented, disease evolution is still quite unpredictable for a given patient. In many cases tissue granulomas resolve spontaneously while in others they persist or reoccur requiring treatment, mostly prednisone alone or other drugs as immunosuppressants [START_REF] Valeyre | Sarcoidosis[END_REF][START_REF] Grunewald | Sarcoidosis. Nature reviews Disease primers[END_REF][START_REF] Schupp | Phenotypes of organ involvement in sarcoidosis[END_REF].

The current therapeutic monitoring is hence mainly based on empirical knowledge and not on well-identified pathogenesis mechanisms. Therefore, a better knowledge of the intracellular pathways potentially activated is needed.

As emphasized by Baughman and Valeyre in the preface of the sarcoidosis book they have coordinated "The world of sarcoidosis is changing on many fronts" [START_REF] Baughman | Sarcoidosis: A Clinician's Guide[END_REF]. The disease phenotypic variability has been recently challenged at both individual and epidemiological levels. Data obtained from the GenPhenReSa (Genotype-phenotype relationship in sarcoidosis) European project allowed a smart subgrouping of sarcoidosis patients according to main organ clusters [START_REF] Schupp | Phenotypes of organ involvement in sarcoidosis[END_REF]. According to their clinical patterns, 2163 European patients were distributed in five distinct subgroups according to predominant organ involvement. In addition, three geoepidemiological factors i.e. geolocation, ethnicity and personal environment from 128 955 patients were compared on a basis of published studies reporting at least 100 patients each and showed different patterns among people of White, Black and Asian origin [START_REF] Brito-Zeron | Geoepidemiological big data approach to sarcoidosis: geographical and ethnic determinants[END_REF].

The authors concluded that "Sarcoidosis is probably one of the systemic diseases with the greatest influence of geoepidemiological factors on the frequency and phenotypic expression". These results emphasize that sarcoidosis is not a univocal disease. Its variability may reflect genetic background and environmental exposure differences as well as socioeconomic status of patients. Therefore, one can postulate that sarcoidosis pathogenesis implies the activation of different intracellular pathways.

The recent identification of significant key cellular pathways can be considered as the beginning of an answer related to the complexity of sarcoidosis. Indeed, the Mammalian Target of Rapamycin (mTOR) pathway [START_REF] Linke | Chronic signaling via the metabolic checkpoint kinase mTORC1 induces macrophage granuloma formation and marks sarcoidosis progression[END_REF], the Janus Kinase/signal transducers and activators of transcription (JAK/STAT) pathway [START_REF] Damsky | Tofacitinib Treatment and Molecular Analysis of Cutaneous Sarcoidosis[END_REF][START_REF] Rotenberg | Dramatic response of refractory sarcoidosis under ruxolitinib in a patient with associated JAK2-mutated polycythemia[END_REF] and the NLRP3 inflammasome pathway in this issue of the journal [START_REF] Huppertz | The NLRP3 inflammasome pathway is activated in sarcoidosis and involved in granuloma formation[END_REF] detected in most sarcoid granulomas [START_REF] Chen | Serum amyloid A regulates granulomatous inflammation in sarcoidosis through Toll-like receptor-2[END_REF], potentially locally produced by macrophage-derived cells [START_REF] Jumeau | Expression of SAA1, SAA2 and SAA4 genes in human primary monocytes and monocyte-derived macrophages[END_REF] is known to activate the NLRP3 inflammasome by involving the purinergic receptor P2X7 and cathepsin B [START_REF] Niemi | Serum amyloid A activates the NLRP3 inflammasome via P2X7 receptor and a cathepsin B-sensitive pathway[END_REF].

Multiple pathways have been described in the pathogenesis of sarcoidosis including autophagy/mTOR [START_REF] Linke | Chronic signaling via the metabolic checkpoint kinase mTORC1 induces macrophage granuloma formation and marks sarcoidosis progression[END_REF][START_REF] Calender | Exome sequencing and pathogenicity-network analysis of five French families implicate mTOR signalling and autophagy in familial sarcoidosis[END_REF] and JAK/STAT [START_REF] Damsky | Tofacitinib Treatment and Molecular Analysis of Cutaneous Sarcoidosis[END_REF][START_REF] Rotenberg | Dramatic response of refractory sarcoidosis under ruxolitinib in a patient with associated JAK2-mutated polycythemia[END_REF][START_REF] Rosenbach | Jak-inhibitors offer promise for a new era of targeted treatment for granulomatous disorders[END_REF] signaling pathways and future investigations are required to provide important leads toward the development of effective therapeutic strategies. Targeting the inflammasome using MCC950 [START_REF] Coll | A small-molecule inhibitor of the NLRP3 inflammasome for the treatment of inflammatory diseases[END_REF] or other molecules is currently under deep investigation in multiple conditions [START_REF] Swanson | The NLRP3 inflammasome: molecular activation and regulation to therapeutics[END_REF] such as cardiac sarcoidosis (MAGiC-ARTClinicalTrials.gov Identifier: NCT04017936). The intersection of granulomatosis and autoinflammatory diseases is a rare occurrence.

The field of autoinflammatory disorders expanded significantly in the last decade as IL-1-related autoinflammatory syndromes heterogeneous group of disorders share responsiveness to IL-1 blockers [START_REF] Pathak | Autoinflammatory diseases: update on classification diagnosis and management[END_REF]. The results obtained by Huppertz et al.

published in this issue of European Respiratory Journal reinforces the possibility to consider sarcoidosis as one of the polygenic autoinflammatory diseases [START_REF] Pathak | Autoinflammatory diseases: update on classification diagnosis and management[END_REF].

  Editorial for: The NLRP3 inflammasome pathway is activated in sarcoidosis and involved in granuloma formation by Huppertz C et al. Eur Respir J. 2020 Jan 16. pii: 1900119. doi: 10.1183/13993003.00119-2019.

  have been brought to light within the last 3 years. Linke et al. showed that activation of the metabolic checkpoint kinase mTORC1 in macrophages by Tsc2 deletion resulted in a model of sarcoid-like granulomatosis in mice[START_REF] Linke | Chronic signaling via the metabolic checkpoint kinase mTORC1 induces macrophage granuloma formation and marks sarcoidosis progression[END_REF]. In addition, the authors identified activated mTORC1 signaling in one third of granulomas obtained from 27 biopsies from patients. These results were comforted by Calender et al. genetic studies on five French families which led to the identification of genes enriched in the regulation of mTOR signaling pathway and autophagy[START_REF] Calender | Exome sequencing and pathogenicity-network analysis of five French families implicate mTOR signalling and autophagy in familial sarcoidosis[END_REF]. As a consequence, a clinical trial targeting mTOR in a pilot study of sirolimus inhibitors in sarcoidosis with cutaneous affection is ongoing (EudraCT Number: 2017-004930-27). A second breakthrough came from Damsky et al. which successfully treated a patient who had refractory cutaneous sarcoidosis with JAK2 oral inhibitors[START_REF] Damsky | Tofacitinib Treatment and Molecular Analysis of Cutaneous Sarcoidosis[END_REF]. The authors also assessed the role of activated JAK-STAT signaling in these granulomatous lesions from consecutive biopsies of this patient and 21 archived tissue samples. Concomitantly, Rotenberg et al. reported a dramatic improvement in a long standing sarcoidosis with multivisceral flares refractory to corticosteroids associated with JAK2-mutated polycythemia using the JAK inhibitor ruxolitinib[START_REF] Rotenberg | Dramatic response of refractory sarcoidosis under ruxolitinib in a patient with associated JAK2-mutated polycythemia[END_REF]. Since, JAK inhibitors are considered promise for a new era of targeted treatment for granulomatous disorders[START_REF] Rosenbach | Jak-inhibitors offer promise for a new era of targeted treatment for granulomatous disorders[END_REF], and clinical trials are ongoing mostly on refractory cutaneous sarcoidosis (ClinicalTrials.gov Identifier: NCT03910543).In this vein, the Huppertz et al. report on the role of NLRP3 inflammasome in sarcoidosis adds a new stone[START_REF] Huppertz | The NLRP3 inflammasome pathway is activated in sarcoidosis and involved in granuloma formation[END_REF]. A role for IL-1β secretion by macrophages in sarcoidosis has been suggested by multiple studies almost 30 years ago but quite overlooked since[START_REF] Wewers | Evaluation of alveolar macrophages in normals and individuals with active pulmonary sarcoidosis for the spontaneous expression of the interleukin-1 beta gene[END_REF][START_REF] Pueringer | The relationship between alveolar macrophage TNF, IL-1, and PGE2 release, alveolitis, and disease severity in sarcoidosis[END_REF][START_REF] Rolfe | Interleukin-1 receptor antagonist expression in sarcoidosis[END_REF]. However, the IL-1β family has recently regained interest in part due to the molecular delineation of the inflammasome pathway. Since its discovery in 2002, the inflammasome is regarded as a key inflammatory pathway with major implications in various conditions and as a potential target for treatments[START_REF] Lamkanfi | Inflammasomes and their roles in health and disease[END_REF][START_REF] Swanson | The NLRP3 inflammasome: molecular activation and regulation to therapeutics[END_REF]. The inflammasome was first coined in 2002 in Molecular Cell by the group of Jürg Tschopp [18]. Employing 293T cell line or cell extracts from human monocytesderived cell line (THP1), Martinon et al. discovered that IL-1β maturation requires caspase-1 (Casp-1) activity through NLRP1-based caspase-activating complex that they call the inflammasome. The NLRP3 inflammasome was discovered by the same group in a study showing that patients suffering from Muckle-Wells syndrome exhibit gain-of-function mutations in NLRP3 gene leading to constitutive macrophage inflammasome activation and mature IL-1β production [19]. Uric acid crystal was the first reported activator of the NLRP3 inflammasome, leading to gout disease [20]. In this manuscript, Christine Huppertz et al. provided evidence of inflammasome activation and IL-1β expression in sarcoid granulomas by showing consecutive sections of transbronchial biopsy from a sarcoid patient in which CD68 expressing cells (used as a monocyte marker) coexpressed active Casp-1 and IL-1β. In addition, IL-1β staining score positively correlates with disease severity. The further show that while untreated BAL cells from healthy volunteers or patients with scleroderma secrete basal levels of IL-1β, cells from severe sarcoidosis patients display increased IL-1β production. Augmented IL-1β and cleaved Casp-1 productions in BAL cells from sarcoid patients were further confirmed using inflammasome activating conditions in vitro. The authors also report that mRNA levels of miR-223, a micro RNA known to downregulate NLRP3, were decreased in alveolar macrophages sorted from BAL cells while NLRP3 mRNA expression levels were indeed increased. The study of noncoding RNA in sarcoidosis is a new field of investigation, with potential critical significance in terms of diagnosis and treatment [21]. The authors then performed experiments employing a mouse model of granuloma to mimic some features of human sarcoidosis by injecting trehalose 6,6 dimycolate (TDM) mixed in incomplete Freund's Adjuvant (IFA). TDM (also known as Cord Factor) is a glycolipid molecule extracted from the cell wall of Mycobacterium tuberculosis or related species. As compared to untreated conditions, BAL cells from TDM treated mice displayed increased IL-1β release following inflammasometriggering conditions in vitro. Immunohistochemical staining of consecutive lung sections showed that macrophages in the granulomas (Iba-1-positive areas) express IL-1β. The authors investigated the role of NLRP3 as well as its downregulator miR-223 by using deficient animals. NLPR3-and miR-223 deficient mice displayed decreased or increased lung granuloma counts, respectively, confirming the pathological role of the inflammasome using this in vivo model. Finally, the authors treated WT animals with the NLRP3 inhibitor MCC950 [22] or with anti-IL-1β antibody to assess the effects of pharmacological interventions targeting the inflammasome pathway. Both MCC950 or anti-IL-1β treatments resulted in reduced numbers of pulmonary granuloma in TDM-treated WT mice. While these data strongly suggest that the inflammasome is active in sarcoidosis, additional studies are required to confirm the exact contribution of IL-1β in the disease, especially at early stages. Interestingly, in this manuscript Christine Huppertz et al. showed that scleroderma patients with severe lung function impairment display no constitutive NLRP3 inflammasome activation. In contrast, numerous NLRP3 inflammasome activation triggers are potential involved in sarcoidosis, from viral or bacterial PAMPS to endogenous molecules resulting from cell stress/death. For instance, the acute phase protein serum amyloid A (SAA)

  1. CNRS, INEM-UMR7355 and University of Orleans, France. 2. Sorbonne Université, Paris, France. 3. INSERM UMR 1272 Université Paris 13, Bobigny, France. 4. Pneumology Department, Hôpital Avicenne APHP Bobigny, France

	Correspondence:


Lauric Thiault
email: lauric.thiault@gmail.com

Damian Weekers

Matt Curnock

Nadine Marshall

Petina Pert

Roger Beeden

Michelle Dyer

Joachim Claudet

Keywords: compliance, enforcement, recreational fishing, risky facilities, wildlife crime

Marine Protected Areas (MPAs) are effective resource management and conservation measures, but their effectiveness is often hindered by non-compliant activities such as poaching. Better predicting poaching risk and spatial patterns is crucial for efficient law enforcement and to ensure MPAs are delivering the outcomes they were established for.

Here, we predicted poaching risk from recreational fishers within fully protected MPAs of Australia's Great Barrier Reef Marine Park (GBRMP). Combining patrol effort data, observed distribution of reported incidents, and geographically-referenced predictor variables, we modeled the occurrence probability of poaching incidents using boosted regression trees and mapped poaching risk at fine-scale. Our findings (i) reinforce the key role of fishing attractiveness, accessibility and capacity in shaping the spatial patterns of illegal recreational fishing; (ii) identify key interactions among XXX and tipping points beyond which poaching is more likely to occur; and (iii) highlight gaps in patrol effort that could be filled for improved resource allocation. The methods developed through this study provide a novel approach to quantify the relative influence of multiple interacting factors in shaping poaching occurrence, and hold promises for replication across a broad range of marine or terrestrial settings.

Introduction

Marine Protected Areas (MPAs) are widely promoted as a tool for resource management and conservation (REFs). While various elements of MPA design and implementation are essential [START_REF] Claudet | Marine reserves: size and age do matter[END_REF][START_REF] Jupiter | Ecosystem-based management in Fiji: Successes and challenges after five years of implementation[END_REF][START_REF] Green | Larval dispersal and movement patterns of coral reef fishes, and implications for marine reserve network design[END_REF][START_REF] Gill | Capacity shortfalls hinder the performance of marine protected areas globally[END_REF], the effectiveness of an MPA is ultimately reliant on its users' compliance with regulations. Yet, ensuring compliance, or infringements, remains a persistent problem, and numerous noncompliant activities (e.g. harvest, waste disposal, dampening, or illegal constructions) continue to occur within many MPAs worldwide (REF?). Of these, illegal fishing inside MPAs (i.e., poaching) is particularly prevalent, and can easily render MPAs ineffective [START_REF] Samoilys | Effectiveness of five small Philippines' coral reef reserves for fish populations depends on site-specific factors, particularly enforcement history[END_REF][START_REF] Guidetti | Italian marine reserve effectiveness: Does enforcement matter?[END_REF][START_REF] Campbell | Weak Compliance Undermines the Success of No-Take Zones in a Large Government-Controlled Marine Protected Area[END_REF]) and erode trust in management (REF).

As with all regulatory frameworks, individual reasons for not complying with rules vary between negligent, opportunistic and intentional offending. To be effective, compliance management should address each level of offending through appropriate strategies including, education, engagement and enforcement, respectively [START_REF] Ivec | Applications of responsive regulatory theory in Australia and overseas: Update[END_REF].

Enforcement is often associated to resource-demanding management activities in MPAs (i.e. due to vessel, personnel and legal costs), resources need to be strategically allocated by effectively targeting patrols to locations and times at which poaching is more likely to be prevalent. This is particularly critical in large MPAs .

It has recently been shown that wildlife crime tend to be concentrated in space [START_REF] Brill | A decade of illegal fishing in Table Mountain National Park (2000-2009): trends in the illicit harvest of abalone Haliotis midae and West Coast rock lobster Jasus lalandii[END_REF][START_REF] Maingi | Spatiotemporal patterns of elephant poaching in south-eastern Kenya[END_REF][START_REF] Arias | Optimizing enforcement and compliance in offshore marine protected areas: a case study from Cocos Island, Costa Rica[END_REF], time [START_REF] Diogo | Catch me if you can: Non-compliance of limpet protection in the Azores[END_REF], on specific ecosystem components [START_REF] Pires | A CRAVED Analysis of Multiple Illicit Parrot Markets in Peru and Bolivia[END_REF][START_REF] Kurland | Forest Policy and Economics The spatial pattern of redwood burl poaching and implications for prevention[END_REF]) and were concerning a restricted set of offenders [START_REF] Weekers | Conservation Criminology: Modelling Offender Target Selection for Illegal Fishing in Marine Protected Areas[END_REF]. Understanding what prevails/drives these patterns can assist in the development of cost-effective surveillance strategies.

Surveillance activities are more effective when adequately and sustainably funded and staffed, and targeted to the right places at the right times [START_REF] Jachmann | Monitoring law-enforcement performance in nine protected areas in Ghana[END_REF][START_REF] Critchlow | Spatiotemporal trends of illegal activities from ranger-collected data in a Ugandan national park[END_REF][START_REF] Critchlow | Improving law enforcement effectiveness and efficiency in protected areas using ranger-collected monitoring data[END_REF][START_REF] Petrossian | Preventing illegal, unreported and unregulated (IUU) fishing: a situational approach[END_REF]. In line with the 'law of crime concentration' [START_REF] Weisburd | The law of crime concentration and the criminology of place[END_REF], poaching is far from being a random activity. Poaching is a highly structured activity driven by the convergence of willing offenders and vulnerable targets at suitable places [START_REF] Moreto | Wildlife crime: An environmental criminology and crime science perspective[END_REF]. 'Poaching hotspots' are formed when these points of convergence are repeated over space and time, revealing an underlying opportunity structure for illegal activities.

Poaching hotspots tend to share a set of common characteristics, or risk factors (REF). In coastal and marine environments, they include the harvest target availability and attractiveness, accessibility (e.g. travel costs and travel time), opportunism (e.g. along MPA boundaries), guardianship effectiveness, and the perceived likelihood and consequence of getting caught [START_REF] Arias | Optimizing enforcement and compliance in offshore marine protected areas: a case study from Cocos Island, Costa Rica[END_REF]Bergseth et al. 2017;[START_REF] Weekers | Risky facilities : Analysis of illegal recreational fishing in the Great Barrier Reef[END_REF][START_REF] Weekers | Conservation Criminology: Modelling Offender Target Selection for Illegal Fishing in Marine Protected Areas[END_REF].

Determining the relative importance of risk factors and patterns of poaching is nonetheless a challenging and often context-specific tasks. While approaches assessing patterns of illegal activities based on raw patterns of incidents reported by patrols are useful tools to identify where crimes are occurring more frequently [START_REF] Hilborn | Effective Enforcement in a Conservation Area[END_REF], they are typically biased towards areas that are routinely surveyed. Failure to account for spatiotemporal variation in surveillance effort runs the risk of systematically over-or underestimating non-compliant activities [START_REF] Keane | The sleeping policeman: understanding issues of enforcement and compliance in conservation[END_REF]. Other approaches have been proposed to explicitly account for detection biases [START_REF] Critchlow | Spatiotemporal trends of illegal activities from ranger-collected data in a Ugandan national park[END_REF][START_REF] Critchlow | Improving law enforcement effectiveness and efficiency in protected areas using ranger-collected monitoring data[END_REF], but they tend to rely on assumptions that may not necessarily hold in marine systems (e.g. on the form of the relationship between predictors and incident occurrence) and fail to account for interactions between risk factors.

Here, we aim at improving patrol efficiency by better understanding and predicting how [potentially complex] interactions among [various] risk factors shape poaching risk. We use spatially-explicit environmental and human predictors combined with patrol-collected incident and monitoring data commonly available across various settings to (1) quantify the relative influence of various risk factors in shaping poaching risk; (2) identify main interactions and critical tipping points; (3) predict poaching risk in 44 no-take MPAs within Australia's Great Barrier Reef Marine Park (GBRMP); and (4) identify potential gaps in patrol surveillance effort.

Methods

Study site

The Cairns Management Region (CMR) is located within Australia's Great Barrier Reef Marine Park (GBRMP) (Fig. 1). The current zoning plan was established in 2003 and came into force on 1 July 2004. It consists of various types of multiple-use areas, but we focused here on Marine National Park (Green) Zones, a key zone type for the managing authority's strategy (Great Barrier Reef Marine Park Authority 2018). Marine National Park (Green) Zones are 'no-take/regulated access' areas [START_REF] Horta E Costa | A regulation-based classification system for Marine Protected Areas (MPAs)[END_REF] where all extractive activities like fishing or collecting are not allowed. They fall within the broader type of Fully Protected Areas (Horta e Costa et al. 2016) (from now on, we refer to Marine National Park (Green) Zones as no-take zones). No-take zones represent about a third of the GBRMP total area (11.7% of the CRM) and thus require extensive surveillance effort. The primary form of offence in Australia's Great Barrier Reef Marine Park (GBRMP) is poaching by recreational fishers in no-take zones (Great Barrier Reef Marine Park Authority 2018), and a number of studies have suggested that such activity may be occurring significantly more frequently than previously thought [START_REF] Davis | Surveillance and poaching on inshore reefs of the Great Barrier Reef Marine Park[END_REF][START_REF] Castro-Sanguino C | Detecting conservation benefits of marine reserves on remote reefs of the northern GBR[END_REF]Bergseth et al. 2017).

Incidents' presence and pseudo-absence

Drawing on spatially explicit occupancy models increasingly used in the predictive ecological community [START_REF] Marmion | Evaluation of consensus methods in predictive species distribution modelling[END_REF], we modelled the spatial distribution of poaching risk within the CMR's no-take zones using Boosted Regression Trees (BRT; [START_REF] Elith | A working guide to boosted regression trees[END_REF] based on observed distribution of reported incidents as a function of geographicallyreferenced predictor variables,. Gradient boosted regression tree approaches such as BRT are increasingly used over statistical approaches for prediction because they better handle interactions among predictor variables and non-linearity than regression-based approaches; both of which were expected to emerge in our case. BRTs also can prevent overfitting by providing regularization [START_REF] Elith | A working guide to boosted regression trees[END_REF].

Presence records (i.e. occurrence of poaching incidents) were obtained from the Field Management Compliance Unit (FMCU) at the Great Barrier Reef Marine Park Authority. The data represents all reported incidents of illegal recreational fishing in CMR's no-take zones for the period January 2015 to March 2019 (n=221, Fig. S1). About 95% of these happened during daylight hours, between 7am and 7pm, reflecting effort data from the patrol vessel.

The presence data used in this study represents reliable records at GPS resolution, with heterogeneous detectability due to heterogeneous monitoring effort across the study area.

To account for this, we assigned a weight to presence points based on monitoring effort, on the basis that incidents detected in highly monitored areas had lower weight than incidents detected in areas that are monitored more rarely (see Supporting Information).

In our case, confirmed absences of incidents (i.e. locations where poaching never occurred) are more difficult to obtain due to the diffuse nature of offenders and the impracticability to monitor the entire area constantly. Therefore, we created artificial absence data (herein pseudo-absence) following guidelines from [START_REF] Cerasoli | Comparing pseudo-absences generation techniques in Boosted Regression Trees models for conservation purposes : A case study on amphibians in a protected area[END_REF]. Specifically, we generated the pseudo-absences using geographically stratified random selection (i.e. based on density estimate of presences) so that the sum of the weights on the pseudo-absences (i.e.

proportional to monitoring effort) equals the sum of those on the presence points (i.e. inverse of monitoring effort). This process yielded a total of 498 pseudo-absence points (Fig. S1).

Predictors of poaching risk

We used a set of ten spatially-explicit predictors relating to environmental and human dimensions to predict the probability of incident occurrence (Table 1). Distance-related predictors (i.e. Accessibility, Features, Islands, Reefs, and Boundary) were derived from the most up-to-date data available on each elements' locations using the cost distance tool in ESRI's ArcGIS 10.5. Bathymetry data (Depth) was obtained from the DeepReef database (https://www.deepreef.org/bathymetry/65-3dgbr-bathy.html). Slope and Aspect were derived from the bathymetry model, using the 'Slope' and 'Aspect' tools in QGIS, respectively. Coral was modeled as the sum of the surrounding living coral patches, described as the number of 15x15m cells dominated by a coral taxon within a 1km radius around each focal cell on the basis of the Benthic cover type map for Reef Top areas of the Cairns Management Area. Finally, Fishing capacity, defined as the overall ability of the recreational fishery to extract resources in a 50km radius, was modeled by summing the number of motorized recreational boats (all size classes) registered in postcodes located within a 50km radius around each cell. Hence, our predictive model provides a static picture of poaching risk, and therefore assumes that other potential drivers are held evenly distributed throughout the study area.

Building a predictive model of poaching risk

We fitted the BRT models with a weighted logistic regression for binary classification against the ten predictors (Table 1) using the {dismo} package [START_REF] Hijmans | dismo: Species Distribution Modeling[END_REF] in the R statistical software version 3.4.0 (R Core Team 2017). This technique requires the specification of three main parameters: the shrinkage parameter, limiting the contribution of the single trees added to the model through the boosting algorithm (tc), the minimum loss reduction required to make a split (lr), and the bag fraction (bf), proportion of data to be selected at each step. In order to identify the best set of parameters, we implemented a two-step tuning process that retained the set of parameters maximizing cross-validated Area

Under the Curve (AUC) (Supplementary Information). We also explored the possibility of eliminating non-informative predictor variables to select the most parsimonious model, which led to the exclusion of the variable Coral. The final model explained 61% of the crossvalidated variance and had an AUC score of 0.93, indicating strong explanatory and predictive performance, respectively.

We quantified the relative interaction strength and significance between the nine remaining predictor variables using 500 bootstrap resampling [START_REF] Pinsky | Fishing, fast growth and climate variability increase the risk of collapse[END_REF].

We calculated the relative influence of the predictor 95% confidence intervals from 1,000 bootstrap replicates of the original dataset. Based on the same bootstrap replicates, we obtained partial dependency plots with 95% confidence intervals to visualize the relationships between the most influential predictor variables and the response (occurrence probability), while keeping all other predictors constant. Maps of the predicted poaching risk (i.e. probability of incident occurrence) were generated from the optimal BRT model's projections over the whole study area at each 50m * 50m cell with a continuous scale 0-1 for each bootstrap replicate, allowing the median predicted poaching risk +/-95% confidence interval to be mapped. Detailed methods used for model building and bootstrapping are provided in Supplementary Information.

Because the model underlying this map accounts for heterogeneous detectability, we were able to overlap the predicted poaching risk with patrol effort and identify potential spatial mismatch. We visualized how predicted poaching risk with patrol effort overlap using a bivariate choropleth map.

Results

Almost 75% of the spatio-temporal variability of incident occurrence was described by four predictors (Fig. 2). Fishing capacity was the most important predictor variable, accounting for 25.7% of the explained variability in incident occurrence. Depth, Accessibility, and Slope explained a broadly similar portion of the variability in incident occurrence, ranging between 17.8% and 14.3%. Boundary, Islands, and Features had smaller contributions to the model prediction (6.1% each). Reefs and Aspect explained little variability of incident occurrence(s). The occurrence probability remained low at low levels of Fishing capacity (25.7% relative influence) and then steadily increased from 100 boats per 50km radius onwards (Fig. 2).

Similar patterns were observed for Depth (17.8%), with initially low levels of poaching likelihood below -40m increasing until reaching a plateau around -20m depth. Accessibility was the third most important predictor of poaching occurrence (15.6%), with a negative non-linear relationship displaying a threshold around 45km from the nearest boat ramp.

Fitted function for Slope (14.3%) displayed a positive asymptotic relationship that reached a plateau around 15° angle. Other less significant predictors with negative relationships were distance to: Boundary (6.1%), Features (6.1%), Islands (6.1%), and Reefs (3.6%).

Models considering interactions between predictors performed better than simpler ones.

The four strongest pairwise interactions were Fishing capacity x Accessibility (71.51; p-value<0.001), Fishing capacity x Depth (36.9; p-value<0.001), Slope x Accessibility (14.1; p-value<0.001) and Fishing capacity x Slope (9.9; p-value<0.01). Occurrence probability for incidents was higher in areas characterized by higher Fishing capacity, shallower Depths, shorter distances to boat ramps (i.e. Accessibility), and steeper sea bottom (Figs. 23). As a consequence of these interactions and thresholds, predicted poaching risk was highly heterogeneous across the study area (Fig. 4). Poaching risk was concentrated on inshore and mid-shelf reefs located near major towns (Port Douglas, Cairns and Innisfail). Conversely, notake zones located far off-shore and in the north of the Cairns Management Region were exposed to lower levels of predicted poaching risk. Poaching hotspots include sites such as Low Isles Reef, Upolu Reef, Green Island Reef, Wide Bay, and Sisters-Stephens Reef (Fig. 4). Patrol allocation was skewed towards a few MPAs around Cairns and, to a lower extent, Port Douglas (Fig. 4). Indeed, among the 44 MPAs located within the CRM, 75% of the total patrol monitoring conducted between 2015 and 2019 was allocated in only four of them (Fig. S2).

This great heterogeneity of patrol effort partially matched with the spatial distribution of predicted poaching risk (Fig. 5). Indeed, while the highest levels of patrol effort were found in the three major poaching hotspots (yellow in Fig. 5), other areas with comparable levels of poaching risk received much less monitoring. Such areas, which may represent enforcement gaps, are mostly located in the south of the CMR (maroon in Fig. 5). Areas predicted to be exposed to low poaching risk tended to be monitored less (blue in Fig. 5).

Discussion

Identifying the underlying drivers of poaching and understanding how they interact to structure poaching risk offer great value to managers for cost effective surveillance. This study presents the first attempt at quantifying the relative importance of, and interaction between multiple risk factors of recreational poaching in a large Marine Protected Area (MPA) using commonly available patrol-collected data. As such, it offers novel insights that can inform management strategies and planning, via a new predictive approach that can potentially be applied to other settings.

The overwhelming influence of only four predictors (Fishing capacity, Depth, Accessibility and Slope) in explaining poaching risk is striking. Our analysis confirms the assumption that poaching 'hotspots' are characterized by substantial Fishing capacity; however, our model identifies a threshold (i.e. >100 boats within a 50km radius; Fig. 2) beyond which poaching risk increases markedly. Similarly, poaching 'hotspots' were predicted by the model at depths shallower than 40m, short(er) distances to the nearest boat ramp (0 to 45km) and in areas of complex topography, defined by steep(er) sea bottom (Figs. 23). Overall, these findings emphasize the value of these simple yet critical features in the assessment of an MPA's likelihood to be exposed to poaching and provide insights into the mechanisms by which they can interact.

Our results highlight the critical role of benthic topography (Depth and Slope) in driving poaching risk. In coral reefs, areas with high slope and low depth -where poaching risk is the highest-typically reflect reef slopes and edges, which often harbor higher target fish abundance and biomass. This suggests that these two variables can be appealing to recreational fishers and therefore may broadly define environmental bounds of "attractiveness" to recreational poachers. Since travel time or cost remain constraining factors [START_REF] Maire | How accessible are coral reefs to people? A global assessment based on travel time[END_REF], Accessibility proved also to be an important driver. Nevertheless, it is worth noting that the tipping point beyond which poaching risk diminishes significantly (around 45km from the nearest boat ramp) remains substantially higher in the GBRMP than in other places globally [START_REF] Daw | Spatial distribution of effort by artisanal fishers: Exploring economic factors affecting the lobster fisheries of the Corn Islands, Nicaragua[END_REF][START_REF] Daw | The spatial behaviour of artisanal fishers: Implications for fisheries management and development[END_REF][START_REF] Metcalfe | Addressing Uncertainty in Marine Resource Management; Combining Community Engagement and Tracking Technology to Characterize Human Behavior[END_REF]. Such difference in maximum travel distance might reflect the relative wealth of Australians and the investment by recreational fishers in faster and more sea-worthy vessels as fishing platforms (CRC 2017(CRC , 2018)).

Poaching risk was better predicted when drivers related to fishers' spatial preference (i.e.

accessibility and attractiveness) interacted with Fishing capacity (Fig. 3). The combined effects of attractiveness, accessibility and fishing capacity in driving fishing pressure generally [START_REF] Daw | Spatial distribution of effort by artisanal fishers: Exploring economic factors affecting the lobster fisheries of the Corn Islands, Nicaragua[END_REF][START_REF] Castro-Sanguino C | Detecting conservation benefits of marine reserves on remote reefs of the northern GBR[END_REF][START_REF] Thiault | Combining participatory and socioeconomic approaches to map fishing effort in small-scale fisheries[END_REF][START_REF] Metcalfe | Addressing Uncertainty in Marine Resource Management; Combining Community Engagement and Tracking Technology to Characterize Human Behavior[END_REF][START_REF] Harborne | Modelling and mapping regional-scale patterns of fishing impact and fish stocks to support coral-reef management in Micronesia[END_REF], and poaching risk specifically [START_REF] Diogo | Catch me if you can: Non-compliance of limpet protection in the Azores[END_REF][START_REF] Weekers | Risky facilities : Analysis of illegal recreational fishing in the Great Barrier Reef[END_REF] have been showed elsewhere. However, the fact that poaching-specific predictors (e.g. distance to the nearest boundary) are more marginal compared to fishing capacity suggests that most incidents may be due to fishers' negligence rather than driven by opportunistic or intentional motivations. This is particularly XXXX given the fishers' positive perception of the current zoning [START_REF] Sutton | Recreational fishers' attitudes towards the 2004 rezoning of the Great Barrier Reef Marine Park[END_REF] Our findings provide the preliminary indication that the patrol strategies applied in this case may have been sub-optimal, as patrolling effort only partially matches with the identified spatial patterns of poaching risk. Designing more cost-effective enforcement strategies may require redistributing partly patrol effort where enforcement gaps are likely to occur (i.e.

higher poaching risk and lower patrol effort). No-take MPAs that may benefit from increased effort are often adjacent to shore. Thus, land-based compliance officers might be deployed more frequently in these areas. Systematic resource allocation methods (e.g., Marxan) could also be used to achieve an optimum deployment of patrol vessels, aiding the design of costeffective enforcement strategies [START_REF] Plumptre | Efficiently targeting resources to deter illegal activities in protected areas[END_REF].

The benefit of our approach provides a nuanced and precise understanding of the interactions among various risk factors related to recreational poaching, allowing one to reliably and accurately predict poaching risk and examine where patrols may be preferentially allocated. The relatively low sample size, however, means that we were not able to incorporate the temporal dimension and identify when such patrols should be deployed. Future applications based on multiple times and long periods would provide more dynamic predictions. Weather effects would for instance enable an improved contextual understanding and help to better predict poacher behavior on a day-to-day basis [START_REF] Critchlow | Spatiotemporal trends of illegal activities from ranger-collected data in a Ugandan national park[END_REF], while longer temporal changes could help determining the deterrence effects of patrols [START_REF] Dobson | Detecting deterrence from patrol data[END_REF], both of which will aid patrol strategy decisions and improve patrol efficiency as a whole. It also important to never forget to invest in all activities that also help promote compliance, such as information or participation to decision-making. Areas in greatest risk of poaching exposure with low surveillance effort are shown in maroon; they correspond with areas potentially requiring more patrolling given their predicted level of poaching risk. These include 1-Sisters-Stephens Reef; 2-Gardens Beach; 3-Feather Reef; 8-MNP-17-1070.

Figure 1 :

 1 Figure 1: Location of the Cairns Management Region (CMR) within the Great Barrier Reef Marine Park, Australia.

Figure 2 :

 2 Figure 2: Predictors of recreational incident occurrence in no-take zones. The left panel shows the relative influence (and 95% confidence intervals) of the predictor variables. The right panel shows partial dependency plots with 95% confidence intervals for the four most influential variables. The graphs show the effect of a given predictor on the probability of incident occurrence while holding all other predictor variables constant at their mean.

Figure 3 :

 3 Figure 3: Pairwise interaction plots of the four strongest interactions between variables predicting poaching risk. The left panels indicate the median fitted function calculated on 1,000 bootstrap replicates. The right panel indicates error around this fitted function, measured as the difference between the 0.95 and 0.05 quantiles calculated on the same bootstrap replicates. All interactions were significant (p<0.01).

Figure 4 :

 4 Figure 4: Predicted relative probability of incident occurrence within the no-take MPAs of the Cairns Management Region of the Great Barrier Reef Marine Park, and associated level of patrol monitoring effort. Numbers indicate poaching hotspots within no-take zones: 1-Tongue Reef; 2-Low Isles Reef; 3-Michaelmas Reef; 4-Upolu Reef; 5-Green Island Reef; 6-Wide Bay; 7-Scott Reef; 8-MNP-17-1070; 9-Normanby-Mabel Reef; 10-Feather Reef; 11-Sisters-Stephens Reef; 12-Gardens Beach.

Figure 5 :

 5 Figure 5: Congruence and mismatch between patrol effort and poaching risk. Bi-plot represent MPA averages. Areas in yellow and blue respectively indicate where patrol distribution matches with poaching risk. Areas in greatest risk of poaching exposure with low surveillance effort are shown in maroon; they correspond with

FiguresFigure 1 :

 1 Figures Figure 1: Location of the Cairns Management Region (CMR) within the Great Barrier Reef Marine Park, Australia.

Figure 2 :

 2 Figure 2: Predictors of poaching risk in no-take zones according to the final BRT model. The top-left plot shows the relative influence of the predictor variables used to model the occurrence of poaching incidents. The other plots show partial dependency plots with 95% confidence intervals for the four most influential variables. The graphs show the effect of a given predictor on the probability of incident occurrence while holding all otherpredictor variables constant at their mean.

Figure 3 :

 3 Figure 3: Pairwise interaction plots of the four strongest interactions between variables predicting poaching risk. The left panels indicate the median fitted function calculated on 1,000 bootstrap replicates. The right panel indicates error around this fitted function, measured as the difference between the 0.95 and 0.05 quantiles calculated on the same bootstrap replicates. All interactions were significant (p<0.01).

Figure 4 :

 4 Figure 4: Predicted relative probability of incident occurrence within the no-take MPAs of the Cairns Management Region, and associated level of patrol monitoring effort. Numbers indicate poaching hotspots within no-take zones: 1-Tongue Reef; 2-Low Isles Reef; 3-Michaelmas Reef; 4-Upolu Reef; 5-Green Island Reef; 6-Wide Bay; 7-Scott Reef; 8-MNP-17-1070; 9-Normanby-Mabel Reef; 10-Feather Reef; 11-Sisters-Stephens Reef; 12-Gardens Beach.

Figure 5 :

 5 Figure 5: Identifying potential mismatch between patrol effort and poaching risk. Bi-plot represent MPA averages. Areas in yellow and blue respectively indicate where patrol distribution matches with poaching risk.

Table 1 |

 1 Predictors used to predict poaching in no-take zones of the Cairns Management Area.

	Name	Description

Using this set of predictors, we were able to capture some previously unexplored potential risk factors in the GBRMP, although we acknowledge that poaching risk can have other dimensions such as individual determinants owned by offenders, weather, or time of the day/week/year

[START_REF] Bergseth | Discerning the culture of compliance through recreational fi sher ' s perceptions of poaching[END_REF] Bergseth et al. 2017;[START_REF] Oyanedel | Illegal fishing and territorial user rights in Chile[END_REF][START_REF] Weekers | Risky facilities : Analysis of illegal recreational fishing in the Great Barrier Reef[END_REF][START_REF] Weekers | Conservation Criminology: Modelling Offender Target Selection for Illegal Fishing in Marine Protected Areas[END_REF]

) that we were not able to incorporate here.

  and their negative attitude towards poaching activities in general (Bergseth & Roscher 2017) XXXXX .

	The critical roles of fishing capacity (determined by the number and location of registered
	recreational boats) and accessibility (determined by the number and location of boat ramps)
	indicate potential benefits associated with increased integration of new and updated data,
	for example, in coordination with the Queensland Department of Transport and Main Roads
	(TMR). Although the number of boats registered cannot be capped, TMR registration data,
	regulatory conditions, and planning schemes (e.g. for development and maintenance of
	recreational access points) represent potentially valuable points of opportunity around
	which to foster collaborative monitoring and management. [It is also worth noting that the
	main predictor of poaching risk in our case (i.e. Fishing capacity) is largely influenced by boat
	ownership, which is in turn driven by the broader external economic context (e.g.
	Queensland's resources boom and associated Fly-In-Fly-Out workforce patterns; Queesnland
	Government, 2015).]

This work was funded by the Great Barrier Reef Marine Park Authority (GBRMPA), Australia and Agence Nationale de la Recherche (ANR XXXX), France.