Davide Faranda
email: davide.faranda@lsce.ipsl.fr

Yuzuru Sato

Gabriele Messori

Nicholas R Moloney

Pascal Yiou

Minimal dynamical systems model of the northern hemisphere jet stream via embedding of climate data

We derive a minimal dynamical systems model for the northern hemisphere mid-latitude jet dynamics by embedding atmospheric data, and investigate its properties (bifurcation structure, stability, local dimensions) for different atmospheric flow regimes. The derivation is a three-step process: first, we obtain a 1-D description of the mid-latitude jet-stream by computing the position of the jet at each longitude using the ERA-Interim reanalysis. Next, we use the embedding procedure to derive a map of the local jet position dynamics. Finally, we introduce the coupling and stochastic effects deriving from both atmospheric turbulence and topographic disturbances to the jet. We then analyze the dynamical properties of the model in different regimes: one that gives the closest representation of the properties extracted from real data; one featuring a stronger jet (strong coupling); one featuring a weaker jet (weak coupling); and one with modified topography. Our model, notwithstanding its simplicity, provides an instructive description of the dynamical properties of the atmospheric jet.

Introduction

Jet streams are narrow, fast-flowing westerly air currents near the tropopause. They are a major feature of the large-scale atmospheric circulation and modulate the frequency, severity and persistence of weather events across the extra-tropics (e.g. [START_REF] Röthlisberger | Regional-scale jet waviness modulates the occurrence of midlatitude weather extremes[END_REF]). Their location and intensity also affects commercial aviation and shipping [START_REF] Reiter | Jet-stream structure and clear-air turbulence (CAT)[END_REF][START_REF] Hadlock | The Experiment on Rapidly Intensifying Cyclones over the Atlantic (ERICA) field study: Objectives and plans[END_REF][START_REF] Williams | Intensification of winter transatlantic aviation turbulence in response to climate change[END_REF]. Two types of atmospheric jets can be identified: thermally-driven subtropical jets, and eddy-driven jets associated with baroclinic instability at the polar front. In the Northern Hemisphere (NH), the two are not always clearly separated [START_REF] Lee | The dynamical relationship between subtropical and eddy-driven jets[END_REF], and when considering monthly or longer time-averages a single, 1 spiral-shaped jet structure emerges (e.g. [START_REF] Archer | Historical trends in the jet streams[END_REF]. In this paper we consider a single NH jet (NHJ), rather than attempting to separate the subtropical and eddy-driven jets (e.g. [START_REF] Belmecheri | Northern Hemisphere jet stream position indices as diagnostic tools for climate and ecosystem dynamics[END_REF].

Even though the climatological NHJ is a westerly flow, it can present large meanders on synoptic timescales (e.g. [START_REF] Koch | An event-based jet-stream climatology and typology[END_REF][START_REF] Röthlisberger | Regional-scale jet waviness modulates the occurrence of midlatitude weather extremes[END_REF]. These can cause the local flow to become predominantly meridional, or can even determine a splitting or breaking of the jet [START_REF] Haines | Isolated anomalies in westerly jet streams: A unified approach[END_REF]. The occurrence of these large meanders in the jet is often associated with events such as temperature and precipitation extremes (e.g. [START_REF] Dole | Was there a basis for anticipating the 2010 Russian heat wave?[END_REF][START_REF] Screen | Amplified mid-latitude planetary waves favour particular regional weather extremes[END_REF].

Although jet dynamics are well understood in a climatological sense, our insights into dynamical features such as jet splitting or meandering are still limited.

The dynamics of meanders and split jets has often been framed in terms of transitions between zonal and blocked flows since the seminal work by [START_REF] Charney | Multiple flow equilibria in the atmosphere and blocking[END_REF]. [START_REF] Legras | Persistent anomalies, blocking and variations in atmospheric predictability[END_REF] and [START_REF] Ghil | Dynamics, statistics and predictability of planetary flow regimes[END_REF] used an intermediate complexity barotropic model with dissipation forcing and topography, and observed two distinct equilibria associated with the zonal and blocked flows. Similar mechanisms have been proposed by [START_REF] Mo | Cluster analysis of multiple planetary flow regimes[END_REF] and then performed in experimental facilities [START_REF] Weeks | Three-dimensional direct imaging of structural relaxation near the colloidal glass transition[END_REF]. However, there is no consensus about the nature of flow multistability, and a wide range of theoretical explanations and models have been proposed (e.g. [START_REF] Tung | A theory of stationary long waves. I-A simple theory of blocking[END_REF][START_REF] Simmons | Barotropic wave propagation and instability, and atmospheric teleconnection patterns[END_REF][START_REF] Frederiksen | A unified three-dimensional instability theory of the onset of blocking and cyclogenesis[END_REF]Faranda et al., 2016b). Moreover, jet dynamics have been described as a manifestation of multiple equilibria in asymmetrically forced flows [START_REF] Hansen | Observational characteristics of atmospheric planetary waves with bimodal amplitude distributions[END_REF] or as a result of soliton-modon structures [START_REF] Mcwilliams | Numerical studies of barotropic modons[END_REF].

In order to advance our understanding of the jet dynamics, we employ a low-dimensional dynamical systems model derived from reanalysis data. The best-known example of a low dimensional model for atmospheric phenomena is Lorenz' simple threedimensional system representing some features of Rayleigh-Bénard convection [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. Thereafter, simple dynamical systems models have been devised to study El Niño [START_REF] Penland | A balance condition for stochastic numerical models with application to the El Nino-Southern Oscillation[END_REF], ocean-atmosphere interactions [START_REF] Dijkstra | Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach[END_REF], climate tipping points [START_REF] Stommel | Thermohaline convection with two stable regimes of flow[END_REF][START_REF] Benzi | Stochastic resonance in climatic change[END_REF], large scale atmospheric motions [START_REF] Lorenz | Irregularity: A fundamental property of the atmosphere[END_REF][START_REF] Lorenz | Predictability: A problem partly solved[END_REF] and many other phenomena. The goal of these investigations was not to provide the most realistic representation of the relevant systems, but rather to capture key emerging behaviours (such as chaos, intermittency, multistability). The main drawback of those investigations was the weakness of the connection between models and real data due to the scarcity of observations as well as theoretical limitations. Until very recently, there was a strong case against the use of embedding techniques to derive low dimensional models from experimental data [START_REF] Letellier | How the choice of the observable may influence the analysis of nonlinear dynamical systems[END_REF]. This opposition was motivated by a long sequel of papers that appeared between 1984 and 1991. The initial claim that low dimensional models for complex phenomena could be derived using a very small numbers of variables (see e.g. [START_REF] Nicolis | Is there a climatic attractor?[END_REF]; [START_REF] Fraedrich | Estimating the dimensions of weather and climate attractors[END_REF]) was disproved by rigorous numerical computations by [START_REF] Grassberger | Do climatic attractors exist?[END_REF] and [START_REF] Lorenz | Dimension of weather and climate attractors[END_REF].

Progress in data quality and availability and the advent of stochastic dynamical systems have renewed the attention for data embedding. Recently, Faranda et al. (2017c) have shown that embedding techniques can yield effective low-dimensional dynamics provided that the chosen observables reflect the symmetries of the system and that small-scale (subgrid) dynamics is represented as stochastic perturbations. Here, we use these results to develop a minimal model of the effective dynamics of the mid-latitude jet. This is useful to explore a range of possible behaviors beyond those displayed in the available data, that could have appeared in past climates and could appear again in future climates. In analogy to the model derived by Faranda et al. (2017c) for the von Karman turbulent flow, the jet model is based on a coupled map lattice (CML, see Appendix A).

Each element of the lattice reflects the dynamics of the jet at each longitude. Such a model does not require physical sub-grid terms a priori, but only if they are found to be essential to capture the large-scale phenomenology -which we show is not the case. We then evaluate how this model represents key dynamical features of the jet, namely its stability, the statistics of splitting/breaking and the response to topographical features, and relate the results back to the original ERA-Interim data.

First, we provide the details of the ERA-Interim data and of the jet detection algorithm (Section 2). We then present the stochastic coupled lattice map model and compute its bifurcation structure (Section 3). Next, we introduce some instantaneous dynamical indicators (Section 4) and use them to relate the conceptual model to more complex climate models and reanalysis data (Section 5). Finally, we highlight the open questions our results can answer and the new questions they pose (Section 6).

Data and Methods

ERA-Interim data and jet position algorithm

The analysis is based on the European Centre for Medium Range Weather Forecasts' ERA-Interim reanalysis [START_REF] Dee | The ERA-Interim reanalysis: Configuration and performance of the data assimilation system[END_REF].

We consider daily data with a 1 • horizontal resolution over the period 1979-2016.

The jet position is diagnosed through a modified version of the approach by [START_REF] Woollings | Variability of the North Atlantic eddy-driven jet stream[END_REF]. We take daily mean wind-speed averaged over 200-400 hPa and apply a 10-day low-pass Lanczos filter [START_REF] Duchon | Lanczos filtering in one and two dimensions[END_REF]. We then identify the latitudinal position of the jet at every longitude as the location of the strongest wind, over the band 15 • -75 • N. This approach is intended to provide a "raw" measure of the jet variability. We then consider the longitude and time dependence of the latitude of the jet to monitor its waviness.

We define an index of large jet meanders, or breaks (Breaking Index, BRI), as the daily number of meridional variations in jet position of more than 10 • of latitude across adjacent longitude gridpoints, except at longitude 0. The analysis has been repeated for BRI thresholds between 5 • and 15 • , with no significant qualitative differences.

Figure 1 shows a snapshot of the jet position on February 4th 1979, together with the time series of the daily jet position recorded in 1979 at longitude 120 • W. An animation of the jet location for the year 1980 is provided as supplementary video.

Both the time series and the snapshot show large jumps in the jet position. A qualitative analysis of the jet position data suggests that the jet fluctuates around a central latitude (Central Jet, CJ) and seldom shifts to more northerly (NJ) or southerly (SJ) latitudes.

In order to embed the data and derive the effective maps of the dynamics, we remove the seasonal cycle from the data by subtracting, longitude by longitude, the average meridional position for each calendar day and dividing by the standard deviation. For the deseasonalized data, the dimensionless threshold for the computation of BRI corresponding to about 10 • latitude is |x| > 1.

Local dynamical systems metrics

Our analysis leverages two recently developed dynamical systems metrics, namely: the local dimension of the attractor d and the stability or persistence of phase-space trajectories θ -1 . Instantaneity in time corresponds to locality in phase-space, such that a value of d and θ -1 can be computed for a given variable (in our case the jet position data) at every timestep. d is a proxy for the system's active number of degrees of freedom. It provides information on how the system can reach a given state, and how it can evolve from such state. θ -1 describes the persistence of a state in time, thus providing complementary information to d.

Local Dimension

The local dimension is estimated by making use of extreme value statistics applied to Poincaré recurrences. The Freitas et al.

(2010) theorem, modified by [START_REF] Lucarini | Extreme value theory for singular measures[END_REF], states that the probability of entering a hyperball with a small radius centred on a state ζ on a chaotic attractor obeys a generalized Pareto distribution [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF]. In order to compute this probability empirically, we first calculate the series of distances dist(x(t), ζ) between the point on the attractor ζ and all other points x(t) on the trajectory. This series is transformed via the distance function:

g(x(t)) = -log(dist(x(t), ζ)), (1)
such that close recurrences of ζ correspond to large values of g(x(t)) [START_REF] Collet | Iterated maps on the interval as dynamical systems[END_REF]. Thus, the probability of entering a small hyperball around ζ is transformed into the probability of exceeding a high threshold s(q), where q is a percentile of the series g(x(t)) itself. In the limit of an infinitely long trajectory, it can be shown that the choice g(x(t)) in Eq. ( 1) locks this probability into the exponential member of the generalized Pareto distribution:

Pr(z > s(q)) ≃ exp [ -ϑ(ζ) ( z -µ(ζ) β(ζ) )] (2)
where where z = g(x(t)), and µ and β (obtained via fitting) depend on the point ζ. These are the location and the scale parameters of the distribution. Remarkably,

β(ζ) = 1/d(ζ), where d(ζ)
is the local dimension around the point ζ. This result has been recently applied to a range of atmospheric and oceanic fields (e.g. Faranda et al., 2017b, a;[START_REF] Messori | A dynamical systems approach to studying midlatitude weather extremes[END_REF]Faranda et al., 2019a, b). In this paper, we use the quantile 0.975 of the series g(x(t)) to determine q. Our results are robust with respect to reasonable changes in this quantile.

Local Persistence

Extreme value statistics also allow estimating the persistence of a given state ζ, by inspecting the temporal evolution of the dynamics around it. A measure of persistence around ζ can be obtained from the mean residence time of the trajectory within the neighborhood of ζ. To quantify this, we employ the so-called extremal index ϑ [START_REF] Freitas | The extremal index, hitting time statistics and periodicity[END_REF]Faranda et al., 2016a):

an adimensional parameter 0 < ϑ(ζ) < 1 which can be interpreted as the inverse of the mean residence time. We can then

compute θ -1 (ζ) = dt/ϑ(ζ)
, where dt is the timestep of our data. Heuristically, if the ith visit to the neighbourhood of ζ lasts τ i consecutive time steps, and N such visits are made in total, then θ -1 ≈ (1/N ) ∑ i τ i . In practice, instead of this naive estimator, we compute the extremal index using the likelihood estimator of [START_REF] Süveges | Likelihood estimation of the extremal index[END_REF]. θ = 0 corresponds to a stable fixed point of the dynamics, so that the trajectory resides an infinite amount of time in the neighbourhood of ζ. θ = 1 corresponds to a trajectory residing in the neighbourhood of ζ for only one time step per visit. The estimate of θ is thus sensitive to the dt used. If dt is too large, the time dependence structure is unresolved and θ will be close to 1. Conversely, if dt is too small, θ will be close to 0. Faranda et al. (2017b) observed that θ varies between 0.3 and 0.5, when dt = 1 day, for sea-level pressure fields over the North Atlantic. In this work, we use the same dt.

3 Derivation of the lattice jet model

Model Framework

While not directly issued from the Navier-Stokes equations, our framework builds on concrete physical hypotheses, namely that: (i) the physics of the jet is the same at every longitude and it is only slightly modified by the presence of topographical constraints, (ii) the jet can experience sudden breaks and shifts from its central position (CJ) to northern (NJ) or southern latitudes (SJ), (iii) the jet must propagate to the west, (iv) smaller scale phenomena, such as turbulence and baroclinic waves, will be introduced in the model only if necessary to reproduce the effective dynamics in the data. This latter point is fundamentally different from the philosophy of direct numerical simulations.

We construct our model starting from the local time series of the non-dimensionalised jet position x measured at each longitude i and time t. We use the simplest possible embedding procedure (see Appendix B), which consists of plotting the return map x (i)

t vs. x (i)
t+1 (an example is shown in Figure 2) and searching for a function f (i) such that x

(i) t+1 = f (i) (x (i) t ).
The first thing to verify is that the same functional form f (i) may be used at all longitudes i. This is equivalent to asking that there is only one dynamics driving the jet independently of the location. With the choice:

f (i) (x) =        -A(A+x) A-c + r (i) , x < -c, sinh(βx) + r (i) , -c ≤ x ≤ c, A(A-x) A-c + r (i) , c < x.
(

) 3
where we have dropped the dependecies of x for clarity, the parameters can be fixed at all longitudes as: β = 0.75, A = 3, and c = sinh -1 (A)/β ≈ 2.4246. Even though the functional form of f (i) is independent of longitude, a dependence on i remains in the form of the parameter r (i) , which represents the effects of topography in terms of spatial inhomogeneities of the local dynamics. As first order approximation, we consider only the difference between land and ocean, and assign one of two discrete values to each r (i) . The choice of the function f (i) is not unique; however the one we propose here is a suitable option that satisfies hypotheses (i) and (ii) above. In order to reproduce the eastward propagation of the jet (hypothesis iii), we introduce the coupled map lattice (CML, see [START_REF] Kaneko | Transition from Torus to Chaos Accompanied by Frequency Lockings with Symmetry Breaking: In Connection with the Coupled-Logistic Map[END_REF] and Appendix A):

x (i) t+1 = (1 -ϵ)f (i) (x (i) t ) + ϵf (i-1) (x (i-1) t
), (i = 1, 2, . . . , N ; t = 1, 2, . . .).

(4)

With this geometry, the dynamics is divided into N = 360 cells. Periodic boundary conditions are applied at N = 360. The dynamics in each cell i is time-independent, but perturbed by the cell i-1 (i.e. its neighbour to the west) with intensity ϵ, which we estimate and scale based on the observed data. This further implies that our reference length-scale in the model corresponds to that of 1 • longitude in the mid latitudes, namely order 100km.

Subgrid feedbacks to jet dynamics

If we perform a numerical simulation of Eq. ( 4), the dynamics is fixed to one of the three states (CJ, SJ, NJ), depending on the value of ϵ. This means that the role of small scales perturbations in triggering the transitions between the states is fundamental.

We therefore have to include an additive noise term ξ (i) t in Eq. ( 4):

x (i) t+1 = (1 -ϵ)f (i) (x (i) t ) + ϵf (i-1) (x (i-1) t ) + ξ (i) t , (i = 1, 2, . . . , N ; t = 1, 2, . . .). ( 5
)
The noise is a fundamental ingredient for the breaking of the jet and the transition between zonal and blocked states, as shown in tank experiments and numerical simulations [START_REF] Jacoby | Generation of inertia-gravity waves in the rotating thermal annulus by a localised boundary layer instability[END_REF]. Physically, noise arises from key sub-grid processes that affect the jet dynamics, such as convection or the interaction between the jet stream and gravity waves [START_REF] Williams | Intensification of winter transatlantic aviation turbulence in response to climate change[END_REF][START_REF] Williams | On the generation mechanisms of short-scale unbalanced modes in rotating two-layer flows with vertical shear[END_REF]. Translated to our model with reference spatial scale of order 100km these phenomena, ranging from a few meters to a few kilometers, imply a perturbation in the range 10 -4 < ν < 10 -3 . Several subgrid parametrization of turbulence exists: the seminal works of [START_REF] Kraichnan | Dynamics of nonlinear stochastic systems[END_REF] and [START_REF] Thomson | Criteria for the selection of stochastic models of particle trajectories in turbulent flows[END_REF] showed that if large scales are represented by a deterministic term, a single random variable can drive the turbulence term. This means that Langevin model representations are appropriate to describe turbulent eddies [START_REF] Mccomb | The physics of fluid turbulence[END_REF][START_REF] Frederiksen | Eddy viscosity and stochastic backscatter parameterizations on the sphere for atmospheric circulation models[END_REF]. Following these ideas, we model ν

(i)
t ∈ [-δ, δ] as an uniform random variable.

However, considering small scale turbulent disturbances to the jet dynamics is not sufficient to reproduce the blocking and breaking of the jet. Even if the introduction of ν as a stochastic term can account for the direct Kolmogorov turbulent cascade [START_REF] Kolmogorov | AN Kolmogorov[END_REF], the jet dynamics is also driven by the effect of an inverse cascade transferring energy to large scales via baroclinic waves [START_REF] Held | A scaling theory for horizontally homogeneous, baroclinically unstable flow on a beta plane[END_REF].

Baroclinic activity is associated with extra-tropical cyclones and anticyclones, on scales of order 10 3 km. These can affect the jet position by several degrees of latitude. Again, there is not a unique way to model baroclinic waves in our framework.

We follow the rationale of multiscale parametrizations as they can be theoretically justified (e.g. [START_REF] Wouters | Multi-level dynamical systems: Connecting the Ruelle response theory and the Mori-Zwanzig approach[END_REF]; [START_REF] Kitsios | Subgrid Parameterizations of the Eddy-Eddy, Eddy-Mean Field, Eddy-Topographic, Mean Field-Mean Field, and Mean Field-Topographic Interactions in Atmospheric Models[END_REF]) and are numerically efficient [START_REF] Faranda | Modelling and analysis of turbulent datasets using Auto Regressive Moving Average processes[END_REF]. The simple introduction of another source of noise η (i) , acting at intermidiate scales (i.e. between the scale of the jet and the scale of turbulence), is enough to obtain a reliable jet breaking dynamics (see Sect. 4.1). The simplest choice for η (i)

t ∈ [-µ, µ] is a block noise taking the same value over bl cells (the one-dimensionalized size of cyclones/anticyclones, see Appendix B), and obeying the uniform distribution.

Another choice for modelling baroclinic disturbances to the jet could be to introduce a second deterministic equation, weakly coupled with the jet position. However, this choice requires additional hypotheses and parameters and does not emerge naturally from the embedding procedure used to derive the dynamics of x.

The minimal sub-grid parametrization can thus be written in the form:

ξ (i) t = ν (i) t + η (i) t , ( 6
)
where ν

(i) t
and η

(i)

t model the effects of small turbulent disturbances and baroclinic eddies, respectively (Figure 3). The noise terms are discussed further in Appendix B.

Local dynamics

Owing to the uni-directional coupling in our model and to the large N , the local dynamics can be approximated by a nonautonomous dynamical system:

x (i) t+1 ≃ f (i) (x (i) t ) + p (i) t , (7)
where a non-autonomous external force p

(i)
t is given by:

p (i) t = ϵ [ f (i-1) (x (i-1) t ) -f (i) (x (i) t ) ] + ν (i) t + η (i) t . ( 8
)
When |p

(i)
t | → 0, the local dynamics has a stable fixed point at x ≈ 0, and two unstable chaotic sets near x ≈ ±2. When p (i) t > 0, the resulting perturbed dynamics may exhibit escape behaviour from the fixed point to the chaotic regions with positive Lyapunov exponents. Figure 4 shows the bifurcation diagrams as a function of κ over land (r (i) = -0.02), and ocean (r (i) = 0, see Appendix B) obtained by approximating the external perturbation p (i) t as a random variable obeying the uniform distribution in [-κ, κ]. They both indicate a bifurcation to chaotic and partially chaotic behaviour (Sato et al., submitted). The different value of r (i) over land gives rise to an asymmetry in the invariant sets, namely the sets delimiting the accessible region of the dynamics with respect to all possible external perturbations p (i) t . In Figure 4, these dynamically reachable regions are depicted in grey, while a realization of the dynamics is depicted by the black dots. With r (i) = -0.02 over land and 0.1574 < κ < 0.1985, there is a small chance to reach SJ positions and no chance to reach NJ positions. This is reflected in the skewed distribution of x (i) t . For the sake of conciseness, we do not report the detailed bifurcation analysis of the local dynamics here. A brief analysis for the global dynamics is presented in Sect. 5.

In this section we compare the ERA-Interim deseasonalized jet position data with numerical simulations of our model. In order to have the same statistical sample as as for the reanalysis, we simulate 37 years of daily snapshots of the jet position. The best fit of our model to the data is obtained, by a trial and error procedure, for the parameters: η = 1.2, bl = 15, ϵ = 0.33, δ = 10 -4 .

We further compare the results of model runs containing all noise terms with runs where individual terms are suppressed in turn: the coupling (ϵ = 0), topography (r = 0), and baroclinic waves (bl = 1).

Spatiotemporal Dynamics

We first consider the latitudinal distribution of the yearly median jet positions at each longitude (dots in Figure 5) and their interannual mean (solid lines in Figure 5). The ERA-Interim dataset (Figure 5-a) presents a negative interannual mean jet position at almost all longitudes, with noticeable zonal asymmetries and a marked interannual variability. The best model run (Figure 5-b) captures both the interannual variability and, thanks to the term r, the longitudinal variations in average location.

A run without coupling (ϵ = 0) is shown in Figure 5-c). In this simulation the dynamics is local, except for the presence of block noise, resulting in a discontinuous jet profile. Unlike the ERA-Interim data, the run with no geography (Figure 5-d) has median values which are roughly symmetric around zero. Finally, the run with suppressed baroclinc activity (Figure 5-e) has a smaller interannual variability than the ERA-Interim data and sharp changes in the median values of x following the geographic constraints.

We next consider the NHJ's shifts from CJ towards NJ or SJ positions. We binarize the dynamics by the detecting all the events such that |x| > 1. Note that this corresponds to breaks in the jet position with the same threshold defined by the breaking index BRI, although there is not exact correspondence. We then assign '0' to all the observations with |x| < 1 (CJ) and '1' to all the others (NJ or SJ). This procedure, known as coding, is widely used in dynamical systems analysis to identify different dynamical phases in complex systems [START_REF] Kaneko | Clustering, coding, switching, hierarchical ordering, and control in a network of chaotic elements[END_REF]. The so-obtained binary spatio-temporal dynamics are shown in Figure 6a-e for all the previously described runs. In the ERA-Interim data (Figure 6a), the switch from CJ to NJ and SJ phases occurs in clusters displaying a characteristic longitudinal extent and temporal persistence. There is also some indication of westerly propagation of the clusters. The best model fit captures the qualitative aspects of this behavior, although the longitudinal coherence is weaker (see section 4.2 below). In the remaining model simulations, the suppression of different noise terms alters either the cluster size or the westerly propagation of the clusters (6c-e). A quantitative analysis of the cluster size spectra is presented in Figure -6f for space clusters and Figure-6g for time clusters. There is a clear power law behavior, reminiscent of a multiscale structure [START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF]. This is coherent with the claim that the underlying jet dynamics is turbulent, with energy at all scales. Despite its simplicity, our model is reproduces this power law behavior. The theoretical reasons are nontrivial and can be related to the possibility of building turbulent cascades starting from simple Langevin equations [START_REF] Wouters | Multi-level dynamical systems: Connecting the Ruelle response theory and the Mori-Zwanzig approach[END_REF][START_REF] Faranda | Modelling and analysis of turbulent datasets using Auto Regressive Moving Average processes[END_REF]. We underline here the necessity of adding ϵ and having bl > 1. Indeed when ϵ = 0 the spatial cluster spectrum consist of discrete peaks with the energy concentrated at precises scales. These are a resonance of the block noise size. When instead we impose bl = 1, we still recover a power-law behavior, but the slope for the temporal clustering strongly deviates from that observed in the ERA-Interim data.

Dynamical indicators

We further assess our model by means of the d and θ metrics described in Sect. 2.2, computed on both ERA-Interim data and the coupled map lattice. We also compare here the statistics of the spatial breaks of the jet, detected via the indicator BRI.

Figure 7 show the box-plots of d (Figure 7-a), and θ (Figure 7-b) for each day in the data set and the breaking index BRI (Figure 7-c). The ranges of values of d and θ for the ERA-Interim data resemble closely those found for sea-level pressure fields over the Northern Hemisphere (Faranda et al., 2017a). This supports the claim that the position of the jet is indicative of large-scale features of the NH atmospheric circulation. Similar claims about the relevance of low dimensional projections in describing the mid-latitude atmospheric circulation are presented by [START_REF] Madonna | The link between eddy-driven jet variability and weather regimes in the North Atlantic-European sector[END_REF]. The model runs can produce average dimensions comparable to those observed in the ERA-Interim data, except for the bl = 1 case. There, the fragmented dynamics leads to a much higher dimension. This is consistent with the spatio-temporal diagrams shown in Figure 6. The models' inverse persistence θ are slightly larger than those observed in reanalysis, but still of order 2 days. Here we can clearly see the effect of the noise suppression (ϵ = 0 and bl = 1) in modifying the dynamical properties by leading to lower persistence.

Finally, we remark that the number of breaks is correlated to the local dimension. This result is consistent with Faranda et al. (2017b), who found that high d match blocking-like atmospheric configurations in the North Atlantic region. For the limiting bl = 1 case, BRI is also correlated with θ: the more breaks, the lower the persistence of the flow.

Bifurcation diagram and jet regimes

The bifurcation diagram in Figure 8 is constructed by plotting the empirical density ρ(x) of the jet position at all longitudes as a function of ϵ. The vertical gray line corresponds to the value of ϵ that best fits the ERA-Interim data. The diagram would look symmetric with respect to x = 0 if r = 0 everywhere, but the addition of geography via r (i) alters the relative proportions of time spent in SJ versus NJ. Specifically, our asymmetric land/ocean distribution implies a southward shift of the average CJ position with increasing coupling. This is reminiscent of the behaviour in the stochastic bifurcation obtained from the approximated local dynamics (Figure 4). By analysing the bifurcation structure of the conceptual model as a function of the coupling coefficient -which mimics the coherence of the jet -we identify three behaviours: (i) a strong and uniform jet where large meridional excursions in the jet location are relatively rare events (ϵ < 0.35), which is close to the jet dynamics as inferred from the ERA-Interim data; (ii) a state with sharp meridional excursions in which the jet is very unstable and on average shifted far to the South (0.6 < ϵ < 0.9); and (iii) an intermediate state of transition between the two. These jet regimes are broadly consistent with those obtained in idealised atmospheric simulations [START_REF] Lachmy | Wave and Jet Maintenance in Different Flow Regimes[END_REF][START_REF] Son | The response of westerly jets to thermal driving in a primitive equation model[END_REF], although here we do not delve into the physical mechanisms underlying the different behaviours. It is also noteworthy that our model qualitatively reproduces a southern jet configuration, even though we provide it with a single NHJ and do not distinguish between eddy-driven and subtropical jets.

We have derived a minimal model of the jet stream position dynamics, based on a stochastic coupled map lattice, by embedding data extracted from the ERA-Interim reanalysis. This procedure innovates over earlier studies e.g. (Faranda et al. (2017c)) by making use of a coupled map lattice derived from a local embedding of the data, and could be adapted to systems with several degrees of freedom. Instead of embedding the data of a global observable in a high-dimensional space, we have constructed the return map for the local position of the jet and then added, via coupling and noise, the physical ingredients identified in previous studies as drivers of the jet dynamics. The conceptual model is then validated and tuned using dynamical indicators of the jets dimension and persistence in the reanalysis data.

Future analyses could apply this approach to the southern hemisphere, where the role of topography is less important than in its northern counterpart. This would allow us to better constrain the influence of topography on the dimension-persistence diagrams. Another possibility would be to use the low-dimensional model to build a surrogate data set of the jet positions and then apply this to atmospheric analogues, so as to construct realistic atmospheric dynamics. Finally, it would be interesting to study whether further projections of the atmospheric dynamics to a lower dimensional space are possible, beyond the model developed here, and to test possible relations between different atmospheric blocking indices and the breaking index BRI defined here.

The analysis we have conducted can however already answer some of the questions left open in Faranda et al. (2017a) and [START_REF] Madonna | The link between eddy-driven jet variability and weather regimes in the North Atlantic-European sector[END_REF] concerning the possibility of reducing the complex mid-latitude circulation dynamics to lowdimensional representations given by blocking indices or conceptual models. The fact that the dimension-persistence diagram of our minimal model qualitatively matches many features obtained for the ERA-Interim jet position and sea-level pressure fields shows that a substantial part of the dynamics projects along a single line (the jet position). This may explain why previous investigations observed relatively low dimensions when considering the full sea-level pressure fields (Faranda et al., 2017b, a).

It also suggests that breaks in the jet are responsible for higher dimensions.

Appendix A: Coupled map lattice

A coupled map lattice (CML, [START_REF] Kaneko | Transition from Torus to Chaos Accompanied by Frequency Lockings with Symmetry Breaking: In Connection with the Coupled-Logistic Map[END_REF]) is given by:

x (i) t+1 = (1 -ϵ)f (x (i) t ) + ϵ 2 [ f (x (i-1) t ) + f (x (i+1) t ) ] , (i = 1, 2, . . . , N ; t = 1, 2, . . .). ( 9
)
where ϵ ∈ [0, 1], x (i) t ∈ R, and f (x) : R -→ R. For our jet dynamics, we adopt the open flow model, which is class of CML with uni-directional coupling [START_REF] Kaneko | Spatial period-doubling in open flow[END_REF]:

x (i) t+1 = (1 -ϵ)f (x (i) t ) + ϵf (x (i-1) t
), (i = 1, 2, . . . , N ; t = 1, 2, . . .).

(

) 10
The CML is a phenomenological model to study complex spatio-temporal dynamics in systems with large numbers of degrees of freedom. The idea is to discretize the dynamics in space and time, while capturing the global phenomenology of the system. CMLs have been successfully applied to processes such as turbulence in thermal convection [START_REF] Yanagita | Coupled map lattice model for convection[END_REF] and turbulent puff in pipe flow [START_REF] Avila | Nature of laminar-turbulence intermittency in shear flows[END_REF]. It is convenient for us to model the jet dynamics leveraging the CML approach because we can extract a local one-dimensional map from the observed time series.

Appendix B: Average return map and noise

To extract the local jet dynamics, we construct an average return map. We first coarse-grain the state space into M partitions

L (i) k (k = 1, 2, .
. . , M ) and let x(i,k) be the midpoint of

L (i) k . Then, we construct a set Y (i,k) = {x (i) t |x (i) t-1 ∈ L (i)
k } (t = 2, 3, . . .) and a return map f (i) via the return plot of (x (i,k) , ⟨Y (i,k) ⟩), where ⟨•⟩ is the average over the elements of Y (i,k) at each longitude i, and at each partition k: 2, . . . , , N, k = 1, 2, . . . , M, (11) where |x| ≤ c. In the region |x| > c, we assume linear reflection effects. As a result, we have the return map f (i) in Eq. (3).

⟨Y (i,k) ⟩ = f (i) (x (i,k) ), i = 1,
Figure 2 illustrates the construction for i = 1 and M = 500.

An important ingredient of the jet dynamics is the presence of topographic obstacles to the mid-latitude zonal flow. Mountain ranges and land-sea boundaries cause meridional deviations in the mean jet location [START_REF] Tibaldi | Orographically induced cyclogenesis: Analysis of numerical experiments[END_REF]. This inhomogeneity can be modeled via a parameter r (i) that mimics this "spatial noise". Since the topography is at most a few kilometers in height, this translates to a perturbation of the order of 10 -3 in the model. Reasonable geographical constraints are therefore: As discussed in Sect. 3.2, noise is a fundamental ingredient in the jet dynamics. The "turbulent noise" term ν relates to physical phenomena in the range from a few meters to a few kilometers, implying a perturbation in the range 10 -4 < ν < 10 -3 , where ν is a random variable obeying the uniform distribution. The second noise term, η, relates to baroclinic activity and we model it as a block noise taking the same value over bl blocks (the one-dimensionalized size of cyclones/anticyclones in our model) with an amplitude of order 1. The latter value is determined empirically as it is an indicative magnitude of the large shifts mid-latitude baroclinic systems can induce in the jet. To determine a realistic length for bl, we reason as follows: given that our model has a reference scale of about 100km, and assuming a typical scale for extra-tropical cyclones of about 3000 km, we then have that bl ≈ 30 blocks. However, the perturbations are associated with the cyclone radius rather than diameter: upstream of the cyclone, the jet will mostly be deviated southwards, while downstream of the cyclone, the jet will mostly be deviated northwards. We therefore take the block perturbation to be of size bl = 15 blocks.

r (i) = -0.02 (i ∈
Owing to the uni-directional coupling in our lattice jet model and to the large N , the local dynamics can be approximated by a non-autonomous dynamical system x

(i) t+1 ≃ f (i) (x (i) t )+p (i)
t , where the external force p

(i) t = ϵ [ f (i-1) (x (i-1) t ) -f (i) (x (i) t ) ] + ν (i) t + η (i) t .
Assuming that the time averages ⟨f (i-1) (x

(i-1) t ) -f (i) (x (i) t )⟩, ⟨ν (i)
t ⟩, and ⟨η (i)

t ⟩ are all 0 by symmetry, we have ⟨p (i) t ⟩ ≈ 0. Thus we recover the average return map given in Eq. (11). For r (i) = -0.02 and 0.1574 < κ < 0.1985, there is a small chance to reach SJ positions and no chance to reach NJ positions.

  land) and r (i) = 0.0 (i ∈ ocean), where land spans the ranges 0 ≤ i < 161 and 239 ≤ i < 301, and ocean spans the ranges 161 ≤ i < 239 and 301 ≤ i < 360. The negative sign for the jet shifts over land is justified by the negative median values of the ERA-Interim jet position anomalies over land (Compare Figure 5-a,b with Figure 5-c where no topography is present).

5

  All the authors were supported by the ERC Grant No. 338965-A2C2. DF and YS were supported by the CNRS PICS Grant No. 74774 and by London Mathematical Laboratory External Fellowships. DF and PY were supported by a INSU-CNRS-LEFE-MANU grant (project DINCLIC). YS was supported by the Grant in Aid for Scientific Research (C) No. 18K03441, JSPS, Japan. GM was supported by grant No. 2016-03724 from the Swedish Research Council VetenskapsrÃědet, and the Department of Meteorology of Stockholm University. YS, GM and NRM thank the LSCE for the hospitality.

Figure 1 .

 1 Figure 1. Snapshot of the jet position extracted from the ERA-Interim dataset on February 4th 1979 and time series of the jet position for the year 1979, recorded at longitude 120 • W.

Figure 2 .Figure 3 .

 23 Figure 2. The average return map extracted from the data at longitude i = 1. This is constructed by coarse-graining the state space at i = 1 into M partitions L (1) k (k = 1, 2, . . . , M ). We then define x(1,k) as the midpoint of the partition L (1) k , andY (1,k) = {x (1) t |x (1) t-1 ∈ L (1)k } (t = 2, 3, . . .). The black dots represent (x(1,k) , ⟨Y (1,k) ⟩) for k = 1, 2, . . . , 500, where ⟨•⟩ is the average over the elements of Y(1,k) , computed based on the observed data. The red line represents the approximated average return map ⟨Y (1,k) ⟩ = f (1) (x(1,k) ), when |x| ≤ c. In the region |x| > c, we assume linear reflection effects. As a result, we have the return map f (1) in Eq. (3).

Figure 4 .

 4 Figure 4. Bifurcation diagrams as a function of κ for (a) land (r (i) = -0.02) and (b) ocean (r (i) = 0.0). The grey regions delimit the accessible region of the dynamics with respect to all possible external forcings. A realization of the dynamics is depicted by the black dots.

Figure 5 .Figure 7 .Figure 8 .

 578 Figure 5. Single-year median location (dotted points) and multi-year average (solid lines) of the meridional jet position for ERA-Interim (a) and model (b-e) data. b) Best-fit model, obtained with η = 1.2, bl = 15, ϵ = 0.33, δ = 10 -4 , c) as in b) but with ϵ = 0, d) as in b) but with r = 0, e) as in b) but with bl = 1. The simulations consist of 37 years of daily jet positions.

Author contributions. DF and YS performed the analysis and derived the conceptual model. GM computed the jet position data. All the authors participated in the writing and the discussions.

Competing interests. The authors declare no competing interest.


Antoine Guillon

Deborah Brea

Emilie Luczka

Virginie Herve

Soujoud Hasanat

Camille Thorey

Magdiel Pérez-Cruz

Juliette Hordeaux

Julie Mankikian

Philippe Gosset

Virginie Hervé

Christelle Coraux

Mustapha Si- Tahar

Dr Mustapha Si-Tahar

Inactivation of the interleukin-22 pathway in the airways of cystic fibrosis patients

Keywords: IL-22, Pseudomonas aeruginosa, cystic fibrosis, respiratory infection, neutrophil proteases

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Cystic fibrosis (CF) is characterized by recurrent bronchial obstruction due to mucus accumulation, bacterial airways infections and persistent inflammation. Pseudomonas aeruginosa (PA) is the predominant pathogen responsible of chronic lung infection, resulting in lung tissue damage [START_REF] Broquet | Interleukin-22 level is negatively correlated with neutrophil recruitment in the lungs in a Pseudomonas aeruginosa pneumonia model[END_REF]. Also, human airway surface epithelial regeneration is abnormal in CF [START_REF] Hajj | Human airway surface epithelial regeneration is delayed and abnormal in cystic fibrosis[END_REF]. This CF feature is still not perfectly understood and could be a contributing factor to recurrent infections. IL-22, a member of the IL-10 family, is an important component of immune-epithelial cell cross-talk and plays a critical role in regulating the maintenance of the mucosal barrier [START_REF] Pociask | IL-22 is essential for lung epithelial repair following influenza infection[END_REF]. IL-22 targets respiratory epithelial cells through its receptor IL-22R to control defense mechanisms and repair of respiratory epithelia damaged by infection [START_REF] Pociask | IL-22 is essential for lung epithelial repair following influenza infection[END_REF]. The role of IL-22 in controlling PA infections is not completely clear. A protective role against PA lung infection likely is debated [START_REF] Broquet | Interleukin-22 level is negatively correlated with neutrophil recruitment in the lungs in a Pseudomonas aeruginosa pneumonia model[END_REF][START_REF] Mear | Candida albicans airway exposure primes the lung innate immune response against Pseudomonas aeruginosa infection through innate lymphoid cell recruitment and interleukin-22-associated mucosal response[END_REF][START_REF] Bayes | IL-22 exacerbates weight loss in a murine model of chronic pulmonary Pseudomonas aeruginosa infection[END_REF]. CF lung disease is characterized by non-resolving inflammation, driven by continuous recruitment of immune cells into airways [START_REF] Hartl | Innate immunity in cystic fibrosis lung disease[END_REF]. In particular, neutrophil infiltration is a classic feature of CF and has been shown to be crucial to fight PA [START_REF] Rada | Interactions between Neutrophils and Pseudomonas aeruginosa in Cystic Fibrosis[END_REF]. Neutrophils are a major source of serine proteases (elastase, cathepsin G, protease-3) [START_REF] Twigg | The Role of Serine Proteases and Antiproteases in the Cystic Fibrosis Lung[END_REF]. Neutrophil serine proteases have beneficial roles in microbial killing but could also be harmful by destructing host tissue components [START_REF] Hartl | Innate immunity in cystic fibrosis lung disease[END_REF]. We previously demonstrated a regulation of IL-22 signaling in chronic obstructive pulmonary disease by the proteolytic action of neutrophil enzymes [START_REF] Guillon | Neutrophil proteases alter the interleukin-22-receptor-dependent lung antimicrobial defence[END_REF]. In the context of CF lung disease, the protease burden overwhelms the anti-protease shields and causes a variety of damaging effects, for instance by cleaving immune proteins [START_REF] Twigg | The Role of Serine Proteases and Antiproteases in the Cystic Fibrosis Lung[END_REF]. IL-22 production is evident in CF patients infected by PA [START_REF] Bayes | T helper cell subsets specific for Pseudomonas aeruginosa in healthy individuals and patients with cystic fibrosis[END_REF], but the functionality of the IL-22-dependent pathway has not been examined in the CF context. Here, we aimed to specifically investigate whether IL-22 could be altered in the protease-rich microenvironment of CF airways.

Materials and methods.

Ethical considerations

Animal experiments were conducted according to protocols approved by the Comité d'Ethique en Expérimentation Animale Nord Pas-De-Calais (agreement number (L'Espace de Réflexion Ethique Région Centre) and was conducted in accordance with the ethical standards of the Helsinki Declaration. Informed written consent was obtained from each participant.

Mice and P. aeruginosa-induced lung Infection

Il-22-/-C57BL/6 mice of both sexes, 6-8 weeks old and wild-type (WT) were challenged intranasally with 50 µL of PBS containing 5 x 10 7 CFU of PAKΔpscF strain or PBS alone.

Bronchoalveolar lavage (BAL) fluid and whole lungs were collected 72 hours post-infection.

The lungs were perfused and immersed in 10 % formalin, dehydrated, paraffin-embedded.

Sections (5 µm) were stained with H&E. A multiparameter cumulative scoring system was developed by the pathologist and validated in mice infected with various pathogens (grading details are provided in Online Table1). Bacterial burden was evaluated by plating serially diluted samples of lung homogenate or BAL.

Human airway tissues

Non-infected human airway tissues were collected after nasal polypectomy from 11 CF subjects (median age: 22; range: 12-26; male gender: 9/11; genotypes: F508del/F508del: 7; F508del/G85E: 1; F508del/I507del: 1; F508del/N1303K: 1; F508del/Y1092X: 1) and 7 non-CF subjects who did not suffer from any other disease at the time of experiments were carried out (median age: 58; range: 40-71; male gender: 6/7).

Human airway epithelial cell dissociation and culture at the air-liquid interface (ALI)

Non-CF and CF human airway epithelial (HAE) cells were obtained after dissociation by 0.05% type XIV collagenase (Pronase E; Sigma Aldrich, St Louis, MO) in RPMI 1640 medium supplemented with 20 mM HEPES (Gibco, Life technologies, Paisley, UK) and antibiotics (200 UI/ml penicillin, 200 µg/ml streptomycin; Gibco) overnight at 4°C. HAE cells from 7 non-CF and 7 CF donors (genotypes: F508del/F508del: 5; F508del/G85E: 1; F508del/N1303K: 1) were suspended in CnT17 medium (CELLnTEc, Bern, Switzerland) and seeded (5x10 4 cells per membrane) on type IV collagen (Sigma Aldrich)-coated porous 12 mm diameter Transwell polyester membranes (Transwell-Clear, 0.4 µm pores; Corning, Acton, MA). Cells were grown in liquid-liquid conditions in CnT17 medium until confluence.

The culture medium was then definitively removed from the upper compartment and the epithelium was allowed to differentiate at the ALI for 35 days in medium consisting in 1:1 DMEM/F12 (Gibco,) and BEGM (Lonza, Walkersville, MD) with the Lonza supplements for hEGF, epinephrine, BPE, hydrocortisone, insulin, triiodothyronine and transferrin, and supplemented with 200 UI/ml penicillin, 200 µg/ml streptomycin (Gibco), 0.1 nM retinoic acid (Sigma Aldrich) and 1.5 µg/ml bovin serum albumin (BSA; Sigma Aldrich). All the cultures were incubated at 37°C in a 100% humidified incubator in the presence of 5% CO2.

RNA isolation, cDNA preparation, and RT-PCR. [11]

Immunohistochemical analysis of airway tissues and ALI cultures.

Immunohistochemistry was performed on formalin-fixed paraffin-embedded non-CF and CF human nasal tissues and well-differentiated ALI cultures. Five µm sections were blocked with 10% BSA in PBS for 30 min at room temperature. The slides were then incubated overnight at 4°C with rabbit polyclonal anti-IL-22R1 (Abcam, Cambridge, United Kingdom; 10µg/ml in 10% BSA in PBS). Slides were washed with PBS and incubated with donkey anti-rabbit Alexa Fluor 488 secondary antibodies (Molecular Probe, Eugene, OR) and cell nuclei stained using DAPI (Sigma Aldrich). Negative controls were performed by omitting the primary antibodies. Images were taken by Confocal Zeiss LSM710 microscope (20X DIC/1.40 oil).

Human airway epithelial cell culture for wound-healing assay.

HAE cells from 4 CF donors (genotypes: F508del/F508del: 2; F508del/I507del: 1; F508del/ Y1092X: 1) and from 3 non-CF donors were suspended in BEGM medium and seeded (13x10 3 cells per cm2) on type IV collagen-coated 24-well plates. When confluence was reached, epithelial monolayers were injured mechanically with a P20 pipette tip. After injury, the monolayers were washed to remove detached injured cells, and fresh BEGM medium containing only 10% of each supplement was added. This medium was supplemented or not with human recombinant IL-22 (20 ng/ml; ImmunoTools GmbH; Friesoythe, Germany). A mark on the Petri dishes allowed us to photograph the wounds exactly at the same place at various times and to measure with ImageJ software, the wound area initially and after repair.

Our results are reported as percentages of the initial wound area. Each experiment was realized in triplicate. Sputum from four CF patients were collected (3/4 colonized by P. aeruginosa, 2/4 colonized by S. aureus). Sputum were homogenized in three volumes (w/v) of PBS, incubated with 1 mM dithiothreitol for one hour, and then centrifuged. Purification of human blood neutrophils and assessment of serine proteinase activity were performed as described [START_REF] Dubois | Influence of DNA on the activities and inhibition of neutrophil serine proteases in cystic fibrosis sputum[END_REF]. Total serine proteinase activity in activated cells was inhibited with 1µM of α1Pi (Preparatis, Krakow, Poland) and Pr3 was specifically inhibited with 10 µM of Ac-PYDAP(O-C6-H4-4-Cl)2 [START_REF] Guarino | New selective peptidyl di(chlorophenyl) phosphonate esters for visualizing and blocking neutrophil proteinase 3 in human diseases[END_REF].

Human recombinant IL22 (50 ng) were incubated at 37°C with 20 µL of homogenized sputum or supernatant of purified neutrophils or 50 µL of different concentrations of purified elastase (Preparatis), cathepsin G (Preparatis) and proteinase 3 (Athens Research, Athens, USA).

Reaction was stopped by adding of 6X Lemli then freezing (-20°C) at different time points. IL-22 concentration was determined using DuoSet ELISA assay kits (R&D Systems). SDS-PAGE was performed using anti-IL22 antibody diluted 1/3000e (Abcam, Paris, France).

Statistical analysis.

Results are expressed as means ± SEM for the indicated number of experiments performed independently. Statistical significance between the different values was analysed by the Man-Whitney test (two groups) or the Kruskall-Wallis test with Dunn's multiple comparison test (more than two groups). Statistical analyses was performed using GraphPad Prism5 ® software.

A p value of < 0.05 was considered to be significant.

Results and Discussion

We first sought to clarify the role of IL-22 in PA-associated pulmonary disease using Il-22-/-and wild-type (WT) mice. Lack of representative models of human lung disease due to PA infection is a significant hurdle to translational research. Administration of PA to murine lungs results in either rapid bacterial clearance or immediate animal death, both situations uninformative in regard to host-pathogen interaction studies. The PA bead model has been widely used to study pulmonary PA infection but has limitations in accurately mimicking human disease. These limitations led us to develop an inflammatory model of nonlethal PA infection in Il22-/-mice with PAK ΔpscF, a PA mutant lacking the type III secretion system (T3SS). Lungs from WT and IL-22-/-mice were examined three days post-infection with PAK ΔpscF. Infected Il-22-/-mice exhibited greater lung injury and more inflammatory cell recruitment (mainly neutrophils) than infected WT controls (Fig. 1A,B

, Online Fig 1).

The infected Il-22-/-mice tended to have a higher bacterial load than infected WT mice, but the difference was not statistically significant (Fig. 1C). These results led us to conclude that the IL-22 pathway is likely important to prevent PA-induced lung damage.

We next took into consideration that pathological IL-22R expression at epithelial surfaces has been previously demonstrated and could contribute to the pathogenesis of certain diseases [START_REF] Guillon | Interleukin-22 receptor is overexpressed in nonsmall cell lung cancer and portends a poor prognosis[END_REF][START_REF] Tohyama | IL-17 and IL-22 mediate IL-20 subfamily cytokine production in cultured keratinocytes via increased IL-22 receptor expression[END_REF]. The IL-22R complex consists of the receptor chains IL-22R1 and IL-10R2.

As the IL-10R2 subunit is ubiquitously expressed, the expression of IL-22R1 by epithelial cells determines cellular sensitivity towards IL-22. We sought to investigate whether IL-22R1 was expressed in CF airway epithelium. Human airway tissues were collected after nasal polypectomy from seven CF and seven non-CF subjects. To determine whether IL-22R1 subunit expression was regulated at the transcriptional level, epithelial cells of all samples were cultured at the air-liquid interface (ALI) until complete differentiation, then processed for qRT-PCR. As shown in Figure 1D, IL-22R1 mRNA expression was similar in CF and control patients. We next assessed the IL-22R1 protein expression at epithelial surface by immunofluorescence analysis of CF and non-CF airways samples, along with non-CF and CF differentiated ALI cultures (Figure 1E). Immunostaining revealed an expression of the IL-22R1 protein at the apical surface of ciliated cells that seemed to be comparable in CF versus control patients. We concluded that the CFTR mutation did not alter the IL-22R1 epithelial expression. In addition, we conducted wound-healing assays on nasal epithelial cells obtained from four CF patients and three non-CF patients. We assessed the effect of IL-22 on epithelial repair and observed that the initial wound area was significantly reduced over a 6-h period in IL-22-treated cultures but not in non-treated ones (Fig. 1F). This result suggested that the IL-22-dependent signaling was functional in ex-vivo airway epithelia, including CF epithelium.

However, since the IL-22 pathway seems to be particularly sensitive to proteases [START_REF] Guillon | Neutrophil proteases alter the interleukin-22-receptor-dependent lung antimicrobial defence[END_REF], we hypothesized that the serine protease-rich environment of CF airway may be hardly compatible with a functional IL-22 signaling. To verify this issue, we assessed the stability of IL-22 in sputum of CF patients (Fig. 2A) or in supernatant of neutrophils (Fig. 2B). Sputum and human neutrophils were collected and conditioned as described previously [START_REF] Dubois | Influence of DNA on the activities and inhibition of neutrophil serine proteases in cystic fibrosis sputum[END_REF][START_REF] Guillon | Pseudomonas aeruginosa proteolytically alters the interleukin 22-dependent lung mucosal defense[END_REF]. Recombinant IL-22 was incubated in CF sputum or supernatant of neutrophils for 5 minutes and its concentration was measured by ELISA. We observed a clear decrease of IL-22 concentrations in these two protease-rich fluids (Fig. 2A,B). Consistently, the degradation of IL-22 by the supernatant of neutrophils was confirmed by western-blot (Fig. 2C). The addition of a specific serine protease inhibitor (α1Pi) prevented IL-22 degradation in both fluids (Fig. 2A,B). We further focused on the three major neutrophil serine proteases and observed that elastase, cathepsin G and protease-3 have the ability to cleave IL-22 (Fig. 2D).

However, protease-3 was more prone to destroy IL-22 in a time-(30 seconds) and concentration (0.01µM)-dependent manner (Fig. 2D,E). Of note, it has been previously demonstrated that the activities of these proteases in CF sputum were in the 10 to 100 nM range, with protease-3 and cathepsin G being the most active (>0.1 µM) [START_REF] Dubois | Influence of DNA on the activities and inhibition of neutrophil serine proteases in cystic fibrosis sputum[END_REF]. Finally, we exposed IL-22 to a neutrophil supernatant in the presence of a specific protease-3 inhibitor and observed that the proteolytic destruction of IL-22 was prevented (Fig 2F).

In conclusion, our findings suggest that IL-22 signaling is likely not directly affected by CFTR deficiency per se, notably because IL-22R1 is normally expressed at the epithelial surfaces of CF airways. However, the inflammatory proteases released by activated neutrophils, in particular protease-3, have the ability to rapidly and completely degrade IL-22.

Consequently, we suspect an acquired deficiency of the IL-22 pathway due to the persistent presence of neutrophils in the lungs of CF patients. These observations may have broader implications for other neutrophilic airway diseases. When confluence was reached, epithelial monolayers were injured mechanically, and fresh medium containing only 10% of each supplement was added. This medium was supplemented or not with human recombinant IL-22 (20 ng/ml). The monolayers were photographed at time 0 and after 6h. The wound areas were measured initially and after repair using ImageJ software. Data are shown as median (C) or mean ± SEM (B,D) and reported as percentage of the initial wound area (F). Ctl: control; ns: non-significant. C57BL/6 wild type (WT) and Il-22-/-mice were instilled intranasally with 5x10 7 CFU PAK ΔpscF (PA mutant deleted for T3SS). Three days post-infection, body weight and survival rate were not statistically different among infected WT mice (n=15) and infected Il-22-/-mice (n=14). Results are from two independent experiments and expressed in mean ± SEM.

Figure 1 .

 1 Figure 1. IL-22 pathway is determinant to limit P. aeruginosa-induced lung damage and to heal wound in airway epithelia of CF patients.

Figure 2 .

 2 Figure 2. IL-22 is degraded by neutrophil serine proteases.

Acknowledgements

We are grateful to Michael Breen for the English language editing.

This work was supported by a grant from the Association Vaincre la Mucoviscidose (grant ID RF20150501363). This agency played no role in study design, data collection and analysis, the decision to publish or preparation of the manuscript.

Declaration of interest

None.