O Coquand
email: coquand@lptmc.jussieu.fr

K Essafi
email: essafi@lptmc.jussieu.fr

J.-P Kownacki
email: kownacki@u-cergy.fr

D Mouhanna
email: mouhanna@lptmc.jussieu.fr

Universal behaviors in the wrinkling transition of disordered membranes

Universal behaviors in the wrinkling transition of disordered membranes

I. INTRODUCTION

A considerable activity has been devoted these last years to understanding both experimentally and theoretically the effects of disorder in membranes, mainly within the contexts of the current study of graphene and graphene-like materials on the one hand and, in a more distant past, of partially polymerized lipid membranes on the other hand. Indeed, the synthesis of graphene [1,2] followed by the discovery of its outstanding mechanical, electronic, optical and thermal properties [3][START_REF] Katsnelson | Graphene: Carbon in Two dimensions[END_REF][START_REF] Amorim | [END_REF][6] has stimulated intensive researches aiming at understanding how the unavoidable presence of defects, vacancies, or adatoms would alter the physical properties of pristine compounds. Also, beyond the mere presence of native imperfections, the introduction of artificial defects, e.g. foreign adatoms or substitutional impurities, with the help of various processes -particle (electrons or ions) irradiation, chemical methods like oxidation or crystal growth -has given rise to the emergence of a whole defect engineering industry aiming at achieving new functionalities for these topologically designed graphene and graphene-like materials [7][8][9][START_REF] Ni | Handbook of Graphene[END_REF]. Among the numerous effects observed one finds: variation (increase or decrease) of electronic conductivity according to the size of the defects, increase of elasticity for moderate density of vacancies and decrease at higher density, decrease of thermal conductance, of fracture strength, enhancement of reactivity, appearance of ferromagnetism and so on [8][9][START_REF] Ni | Handbook of Graphene[END_REF][START_REF] López-Polin | [END_REF][12]. As part of this defect engineering activity, a specific effort involving various experimental or numerical techniques -(low pressure) chemical vapor deposition [13], ion/electron irradiation [14][15][16][17][18][19] or molecular dynamic simulation [20][21][22] -has been made toward the design of defect-induced two-dimensional (2D) amorphous counterparts of graphene and graphene-like materials. A highlight of this activity is the achievement by electron irradiation of a step-by-step, atom-by-atom, crystal-to-glass transition giving rise to a vacancy-amorphized graphene structure [13][14][15] similar to the continuous random network proposed by Zachariasen [23]. Many characteristics of this transition have been determined: the onset of the defect-induced amorphization process, its temperature dependence, the structural response to vacancy insertion, the nature of the electronic density of states of the defective configurations [20], a transition in the fracture response from brittle to ductile when increasing vacancy concentration [24]; finally a careful analysis of the glassy-graphene structure in terms of a proliferation of nonhexagonal carbon rings has been performed [15,22]. However we emphasize that the very nature of this glass transition is still unclear. Moreover there has been, up to now, neither within this last context nor within the more general one of the investigation of defective graphene and graphene-like materials, no characterization of a quantitative change between -still putative -ordered and disordered phases and a fortiori no indication of universal behaviors associated with them.

In marked contrast with this situation, in a very different context, recent investigations of partially polymerized lipid membranes by Chaieb et al. [25], following the pioneering work of Mutz et al. [26], have led to identify a remarkable folding-transition while varying the degree of polymerization. More precisely these authors have shown that, upon cooling below the chain melting temperature, partially polymerized phospholipid vesicles undergo a transition from a relatively smooth structure, at high polymerization, to a wrinkled structure, at low polymerization, characterized by randomly frozen normals. This has led them to suggest that this transition would be the counterpart of the spin-glass transition occurring in disordered spin systems [26,27]. Chaieb et al. [25], by considering the height-height correlation functions, have been able to characterize quantitatively the various phases as well as the wrinkling transition separating them. However, despite the large amount of theo-retical work oriented towards understanding the physics of disordered membranes, no theoretical explanation has been given so far on the grounds of these results [27][28][29][30][31][32][33][34][35][36][37][38].

In this article, we show that a recent nonperturbative renormaliation group (NPRG) study, performed by the present authors [39], of the effective theory used to study both curvature and metric disorders perfectly accounts for this situation. In a first part, we recall the experimental status of wrinkled partially polymerized membranes. In a second part, we lay out the unusually unsettled state of the theoretical situation. Finally, in a third part, performing an analysis of the long-distance morphology of membranes at and in the vicinity of the wrinkling transition, we show how the NPRG approach reproduces the experimental outputs. Finally we conclude, stressing the consequences of our analysis for the physics of graphene and graphene-like materials and claiming, in particular, that the behaviors observed in partially polymerized lipid membranes should also be observed in these materials.

II. WRINKLING TRANSITION IN PARTIALLY POLYMERIZED MEMBRANES

The identification of a wrinkling transition in partially polymerized membranes goes back to the work of Sackman et al. [40] on mixture of diacetylenic phospholipids and dimyristoylphosphatidylcholine, followed by those of Mutz et al. [26] and Chaieb et al. [25] on diacetylenic phospholipids [1,2-bis(10,12-tricosadiynoyl)sn-glycero-3-phosphocholine], who have taken advantage of the fact that, upon a chemical or photochemical process, notably ultraviolet (UV) irradiation, these compounds polymerize. In the case considered in [26] the polymerizable phospholipids are first prepared as giant vesicles and then cooled below the chain melting temperature T m 40 • C where they form tubular structures that are then partially polymerized by UV irradiation. The membranes are then reheated above T m where they reform spherical vesicles provided the degree of polymerization does not exceed the percolation threshold located around 40 %. These vesicles, of typical size ranging from 0.3 to 40 µm, are then cooled down to T w 18 -22 • C where they undergo a spontaneous, reversible, phase transition from a relatively smooth structure to a wrinkled, highly convoluted, rigid one displaying locally high spontaneous curvature. This observation has led Mutz et al. [26] to conjecture that this state of affairs should be well described by a theory of polymerized membranes submitted to quenched curvature disorder. The outcomes of this experiment have been made more quantitative by Chaieb et al. [25,[41][42][43] who have studied the transition by various techniques. Small angle neutron scattering has been used to investigate the local structure, giving access to the fractal dimension while environmental scanning electron microscopy has been employed for the study of the surface topography at mesoscopic scale [25,[41][42][43]. Finally a tapping-mode atomic force microscope has provided information on the mean-square fluctuations of the surface height h(x) at a point x = (x 1 , x 2 ), relative to the mean surface height, (h(x) -h(0)) 2 , and its Fourier transform, the power spectrum P (k) [25]. This last quantity has been found to display a remarkable power-law behavior in the range 0.1 -100 µm -1 : P (k) ∼ k -γ where the power exponent γ [64] is directly related to the roughness exponent ζ by: γ = 1 + 2ζ. Three clear distinct regimes have been observed [25] as the degree of polymerization φ is varied, see Fig. 1. At low polymerization, typically for φ < 30%, the surface of the vesicles presents -at large scaleslarge deformations, creases, of order of the vesicle size (500 nm) typical of a wrinkled state. In this case one finds [25]: γ = 2.9 ± 0.1 corresponding to ζ = 0.95 ± 0.05. At high polymerization, typically between 32 and 40% the vesicles are regular at large scales and the creases are less pronounced (of order 20 nm) and one has [25] : γ = 2 ± 0.06 and ζ = 0.5 ± 0.03. Finally, for φ in the intermediary region 30% ≤ φ < 32% a transition occurs and the vesicles display the morphology of a crumpled elastic sheet with γ = 2.51 ± 0.03 and ζ = 0.75 ± 0.02 [25].

FIG. 1: The three scaling behaviors of the power spectrum P (k) as function of k for various degrees of polymerization φ of the membrane. From top to bottom: φ = 40%, φ = 30% and φ = 9% and the corresponding membrane configurations. From Chaieb et al. [25] with permission.

III. THEORETICAL CONTEXT

Early investigations of the wrinkling transition by Mutz et al. [26] have triggered an impressive series of theoretical works aiming to understand the effects of quenched disorder contributions in the seminal model of Nelson and Peliti [44] used to describe the flat phase of disorder-free polymerized membranes [START_REF]Proceedings of the Fifth Jerusalem Winter School for Theoretical Physics[END_REF][START_REF] Bowick | [END_REF][47][48][49][50][51][52][53][54][55]. This series has been initiated by Nelson and Radzihovsky [28,29] who have mainly investigated the effects of impurityinduced disorder in the preferred metric tensor. They have in particular shown that, below the dimension D = 4 of the membrane, the flat phase of disorderfree membranes remains stable at any finite temperature T but should be destabilized at vanishing T for any amount of disorder due to a softening of the bending rigidity, making possible the emergence of a spin-glass behavior. This scenario has been strengthened by the work of Radzihovsky and Le Doussal [30] who, studying the limit of large embedding dimension d of the model, have identified an instability of the flat phase toward a spin-glass-like phase characterized by a nonvanishing Edwards-Anderson order parameter [56]. At the same time Morse et al. [57,58], extending the work done in [28,29], have considered the role of both curvature and metric quenched disorders. Using a perturbative, weakcoupling, = 4 -D expansion, they have shown that the curvature disorder gives rise to a new fixed point at T = 0, stable with respect to randomness but unstable with respect to the temperature. These works have been followed by an intensive search for various kinds, i.e. flat or crumpled, of glassy phases by means of mean-field approximations involving either short-range [27,[30][31][32][33][34][35] or long-range disorders [36,37] (see also [38] for a review) that have led to predict that, at sufficiently strong disorder, the flat phase could be even unstable at any finite temperature toward a glassy phase that would correspond to the experimentally observed wrinkled phase. However, we would like here to emphasize several facts. The first one is that these last approaches, based on mean-field, large d, computations, should be considered with great caution when their conclusions are extended in the finite d case, particularly when they involve a breakdown of the symmetry between replica used to perform the average over disorder. Second, none of the approaches -involving 1/d or = 4 -D expansions -has been carried out at next-to-leading order, where new physics could emerge. Finally, none of them has produced quantitative predictions or explanations as for the results of Chaieb et al. [25].

IV. NPRG ANALYSIS

Recently however, following previous works on disorder-free polymerized membranes [50,52,55,[59][60][61] the present authors [39] have performed a NPRG approach of the model considered within a perturbative framework by Morse et al. [57,58] and whose action is given by

S[R] = d D x κ 2 ∂ 2 i R(x) 2 + λ 2 u ii (x) 2 + µ u ij (x) 2 -c(x).∂ 2 i R(x) -σ ij (x) u ij (x) .
(1) In this action R(x) is a d-dimensional vector field parametrizing, in the embedding space, the points x ≡ x i , i = 1 . . . D of D-dimensional membrane while u ij is the strain tensor that represents the fluctuations around a flat reference configuration R 0 :

u ij = 1 2 (∂ i R.∂ j R - ∂ i R 0 .∂ j R 0 ) with [65] R 0 = [ R(x) ] = x i e i (2)
where . . . and [. . . ] denote thermal and disorder averages respectively. In Eq.( 2) the e i form an orthonormal set of D vectors. The coupling constants κ and (λ, µ) represent respectively the bending rigidity and the Lamé coefficients. The action (1) includes curvature and metric disorder contributions induced by two random fields c(x) and σ ij (x) that couple to the curvature and strain tensor respectively. They are considered here as shortrange, gaussian fields with [57,58] [

c i (x) c j (x )] = ∆ κ δ ij δ (D) (x -x ) σ ij (x) σ kl (x ) = (∆ λ δ ij δ kl + 2∆ µ I ijkl ) δ (D) (x -x ) (3
) where I ijkl = 1 2 (δ ik δ jl + δ il δ jk ), with i, j, k, l = 1 . . . D where ∆ κ and (∆ µ , ∆ λ + (2/D)∆ µ ) are positive coupling constants associated with curvature and metric disorders. Note finally that the ansatz (1), albeit reduced to four powers of the fields and field derivatives, is expected to lead to predictions not altered by higher orders, as this happens quite remarkably in the disorder-free case [52,55,60] and as this is strongly suggested by the very weak sensitivity of our results with the changes of renormalization group (RG) process -see below.

The RG equations corresponding to action (1) have been derived first perturbatively in [57,58] and then within a NPRG approach in [39]. Within this latter approach the RG equations have revealed that there exist, in the space of coupling constants, not only two, as found by Morse et al. [57,58] but actually three nontrivial fixed points: the usual finite-T , vanishing-disorder fixed point P 4 associated with disorder-free membranes [47]; a vanishing-T , finite-disorder fixed point P 5 identified for the first time in [57,58]; finally a finite-T c , finite-disorder fixed point P c found in [39], missed within previous approaches, unstable with respect to T , thus associated with a second-order phase transition and making the T = 0 fixed point fully attractive provided T < T c . The consequences of these facts are twofold: 1) a whole "glassy phase" associated with the T = 0 fixed point is predicted in agreement with the wrinkled phase observed in [25,26] and 2) three distinct universal scaling behaviors are expected, in agreement with the observations of Chaieb et al. [25]. The subsequent analysis shows, moreover, the quantitative agreement between the scaling behaviors predicted and those observed. The quantity to consider is the roughness exponent ζ. Let us recall how this quantity is defined in a field-theoretical context. Let us decompose R(x) around the flat phase configuration R 0 (x) according to R(x) = R 0 (x)+u(x)+h(x) where u(x) and h(x) are respectively the in-plane -phonon -and out-of-plane -flexuron, -degrees of freedom parametrizing the fluctuations around R 0 (x). Writing δh(x) = h(x) -h(x) one defines the connected and disconnected correlation functions of the h field by:

(h(x) -h(0)) 2 = T χ(x) + C(x)
where T χ(x) = [ (δh(x)-δh(0)) 2 ] and C(x) = [ h(x)h(0) 2 that, as usual, respectively measure the thermal and disorder fluctuations. The long-distance behavior of these correlation functions is typically given by: T χ(x) ∼ |x| 2ζ and C(x) ∼ |x| 2ζ that define two roughness exponents ζ and ζ . In the same way correlation functions are defined for the phonon field with two roughness exponents ζ u and ζ u . They are related to the previous ones by: ζ u = 2ζ -1 and ζ u = 2ζ -1. Similarly, in momentum space, writing δh(q) = h(q) -h(q) , one defines:

G hh (q) = h(q)h(-q) = T χ(q) + C(q) (4)
where T χ(q) = δh(q)δh(-q) and C(q) = h(q) h(-q) that behave, at low momenta, as χ(q) ∼ q -(4-η) and C(q) ∼ q -(4-η ) where η and η are the anomalous dimensions evaluated at the fixed points of the RG equations. As a consequence of expression (4) the scaling behavior expected for the height-height correlation function G hh (q) is determined by the relative value of η and η together with the position of the fixed point: at finite or at vanishing T [66]. These exponents are related to the roughness exponents by: ζ = 1 2 (4-D-η) and ζ = 1 2 (4-D-η ) and to the power exponent γ = 5-D-η or γ = 5 -D -η depending on the exponent -η or ηthat controls the long-distance behavior.

At the fully attractive, vanishing-T , finite-disorder fixed point P 5 we find, by improving the results of [39], η 5 = 0.448(2) and η 5 = 0.275 (2) [67]. This is in contrast with both [57,58] and [38] where η 5 = η 5 so that P 5 was found to correspond to a marginal -in fact marginally unstable -fixed point. As a consequence we find a roughness exponent ζ 5 = 0.862(1) that, according to (4) and the scaling laws of χ(q) and C(q), controls the long-distance behavior of the height-height correlation function G hh (q). This corresponds to a power exponent γ 5 = 2.725(2). This value is very close to that found in [25] at low polymerization -for φ in [10%; 30%] -and lies in the range [2.80, 2.92], see Fig. 2. As done in [25] we exclude the data point at lowest polymerization (corresponding to γ = 3 in Fig. 2), which does not belong to the plateau identified for φ in [10%; 30%]. Note that such a value would correspond to the expected value for a fluid membrane, which could be a hint that below φ = 10%, partially polymerized lipid membranes do not behave as polymerized membranes anymore. At the stable, finite-T , finite-disorder fixed point P 4 , we find η 4 = η 4 /2 = 0.849(3) (see [39,50,60] and also [START_REF] Amorim | [END_REF] for a review of other approaches) that corresponds to ζ 4 = 0.575(2) and γ 4 = 2.151(3). As seen in Fig. 2 this is again in good agreement with the results found in [25] at high polymerization with a value of γ that saturates at 2. This value corresponds to the case of disorder-free polymerized membranes. It is in agreement with that obtained by Locatelli et al.

V. CONCLUSION

The conclusion of our work is fourfold. First, we have shown that the longstanding problem of the wrinkling transition taking place in partially polymerized lipid membranes is both qualitatively and quantitatively clarified by means of the NPRG approach used in [39]. Second, reciprocally, this agreement validates the NPRG approach to disordered polymerized and, in particular, consolidates the prediction of the existence of three nontrivial fixed points in the RG flow of the model (1), in contradiction with all previous works. At the formal level this situation raises the question of the origin of the mismatch between the NPRG approach and the previous ones: the weak-coupling approach around D = 4 performed by Morse et al. [57,58] and the self-consistent screening approximation used by Radzihovsky and Le Doussal [30,38]. Third, as the three different kinds of scaling behaviors predicted in [39] are associated with fixed points or RG flow, they are universal and should be observed in a large class of defective materials able to display curvature disorder [START_REF] Morse | [END_REF]. This is in particular the case of defective graphene, whose sp 2 -hybridized carbon structure can reorganize into a non-hexagonal structure displaying nonvanishing curvature. Fourth, and finally, the glassy graphene configurations observed during the vacancy-amorphization process have been shown to display a rough, static, wrinkled structure with reduced thermal fluctuations with respect to their purely crystalline counterpart and exhibit a root mean squared roughness increasing with vacancy concentration indicating a change in the macroscopic morphological/shape structure of defective graphene sheets [20,22]. It would be of considerable interest to see if this transition can be moved closer to the wrinkling transition observed in partially polymerized membranes.

  [62] who have found ζ = 0.54±0.02 by means of low energy electron diffraction on free-standing graphene sheet. Finally, at the finitedisorder, finite-T , critical, fixed point P c found in[39] we get η c = η c = 0.490(2) that corresponds to ζ c = 0.755(1) and γ c = 2.510(2) that is in very good agreement with the value γ = 2.51 found by Chaieb et al.[25] at critical polymerization, see Fig.2.

FIG. 2 :

 2 FIG. 2: Power exponent γ as a function of the polymerization rate φ. Adapted from Chaieb et al. [25]. Horizontal lines correspond to the data extracted from our NPRG computations: γ5 = 2.725(2), γc = 2.510(2) and γ4 = 2.151(3).

  ACKOWLEDGEMENTS D.M. thanks F. Banhart, D. Bensimon, J.-N. Fuchs, A. Locatelli, J. Meyer and O. Pierre-Louis for fruitful discussions.


Adrien Le Coënt

Laurent Fribourg

Jonathan Vacher

Rafael Wisniewski

Probabilistic reachability and control synthesis for stochastic switched systems using the tamed Euler method

Keywords: Stochastic systems, numerical simulation, control system synthesis, switched control systems, nonlinear control systems

Introduction

Symbolic methods for the verification and control synthesis of hybrid systems (and, particularly, "switched systems") have received significant attention in the past few years. However, control systems involving stochastic differential equations remain difficult to handle with symbolic methods, and few methods have been developed for these systems. One distinguishes two main classes of symbolic methods for hybrid systems: indirect methods and direct methods [START_REF] Asarin | Effective synthesis of switching controllers for linear systems[END_REF].

Indirect methods proceed by constructing a finite abstraction of the original system by discretization of the dense state space R d (where d is the dimension of the state space). Among the indirect methods, one of the most successful proceeds by approximate bisimulation [START_REF] Girard | Approximately bisimilar symbolic models for incrementally stable switched systems[END_REF]. This method originally designed for deterministic switched systems has been recently extended for stochastic switched systems [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF][START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF][START_REF] Zamani | Towards scalable synthesis of stochastic control systems[END_REF]. This approach relies on the hypothesis of incremental stability of the stochastic switched system (or existence of a common/multiple Lyapunov function). Another method is presented in [START_REF] Esmaeil | Formal abstractions for automated verification and synthesis of stochastic systems[END_REF]. The associated tool [START_REF] Soudjani | Faust ˆ2: Formal abstractions of uncountable-state stochastic processes[END_REF] allows to generate formal abstractions for discrete-time Markov processes defined over uncountable (continuous) state spaces.

Indirect methods have also been developed for the verification of safety. Safety verification has been studied using the concept of barrier certificates [START_REF] Prajna | Barrier certificates for nonlinear model validation[END_REF]. The original inquiry has been later modified to cope with stochastic systems [START_REF] Wisniewski | Stochastic safety analysis of stochastic hybrid systems[END_REF], and to the design of a controller that makes the closed-loop system safe [START_REF] Wieland | Constructive safety using control barrier functions[END_REF][START_REF] Zakiyullah | Stabilization with guaranteed safety using control lyapunov-barrier function[END_REF]. Also related is the subject of stability and safety verification is the reach avoidance problem. It has been studied using dynamic programming. This approach has been developed both for deterministic hybrid system in [START_REF] Margellos | Hamilton-Jacobi formulation for reach-avoid differential games[END_REF], and for discrete-time stochastic hybrid systems in [START_REF] Summers | Verification of discrete time stochastic hybrid systems: A stochastic reach-avoid decision problem[END_REF].

A direct method proceeds by working directly at the level of the dense state of the system R d ; it computes "trajectory tubes", which are over-approximations of the set of all the controlled trajectories starting at a given subregion of R d . In previous work, we have followed such a direct approach first in [START_REF] Fribourg | Finite controlled invariants for sampled switched systems[END_REF], then in [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF], using the Euler approximation scheme for calculating over-approximations of tubes of trajectories. We show here how to extend this direct method to stochastic switched systems. The method extends the deterministic approach by replacing the classical Euler approximation scheme with a variant of the stochastic Euler-Maruyama scheme [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]. The correctness of these Euler-based methods does not rely on the hypothesis of incremental stability as in [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF][START_REF] Zamani | Towards scalable synthesis of stochastic control systems[END_REF], but on the hypothesis of 'one-sided Lipschitz (OSL)' condition with constant λ P R d (also called 'monotonicity'/'dissipativity', see [START_REF] Von Renesse | Existence and uniqueness of solutions of stochastic functional differential equations[END_REF]). It can be seen that if a stochastic switched system satisfies an OSL condition with λ ă 0, then the function V px, x 1 q " }x ´x1 } 2 is a common incremental Lyapunov function in the sense of [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF], from which it follows that the switched system is incrementally stable, and can be analysed by approximate bisimulation. However, Euler-based methods also apply when the system is not incrementally stable, in which case the constant λ is necessarily positive. We thus consider a class of systems different from that of [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF].

The plan of the paper is as follows: In Section 3, we give an explicit upper bound on the mean square error of the tamed Euler method for SDEs under OSL condition (Theorem 2). We apply the result in order to ensure reachability properties of stochastic switched systems with guaranteed minimum probability (Section 4). We conclude in Section 5.

This paper is an extended version of a paper (A. Le Coënt, L. Fribourg and J. Vacher. Control Synthesis for Stochastic Switched Systems using the Tamed Euler Method. In ADHS'18, IFAC-PapersOnLine 51 [START_REF] Liberzon | Switching in systems and control[END_REF], pages 259-264. Elsevier Science Publishers, 2018). The main additional material of the extended version is:

• Section 3.3, which gives analytical upper bounds on expressions occurring in Theorem 2,

• Section 3.4, which gives a lower bound on the probability of reaching a given set R (Proposition 3), and

• Section 4.4, which proposes an optimized method of control synthesis.

Preliminaries

Notations

The symbols N, N ě0 , R, R ą0 , R ě0 denote the set of natural, nonnegative integer, real, positive, and nonnegative real numbers.

The symbol }¨} denotes the Euclidean norm on R d . The symbol x¨, ¨y denotes the scalar product of two vectors of R d . Given a point x P R d and a positive real r ą 0, the ball Bpx, rq of centre x and radius r is the set ty P R d | }x ´y} ď ru.

Throughout this paper, the set U denotes a finite set of switched modes. Given a sequence of k modes a " pa 1 ¨¨¨a k q P U k , and a sequence of p modes b " pb 1 ¨¨¨b p q P N p , we denote the concatenation a ˚b of a and b the sequence of length k `p: a ˚b " pa 1 ¨¨¨a k ¨b1 ¨¨¨b p q.

A sequence of non-negative integers is written in bold when it is used as a multi-index, e.g. a " pa 1 , . . . , a p q P N p ě0 , the sum of the components of a is denoted by |a| " a 1 `¨¨¨`a p , and `b a ˘" b! a1!a2!...ap! is the multinomial coefficient.

Systems considered and assumptions

Let τ P R ą0 be a fixed real number, let pΩ, F, Pq be a probability space with normal filtration pF t q tPr0,τ s , let d, m P N let W " pW p1q , . . . , W pmq q : r0, τ s ˆΩ Ñ R m be an m-dimensional standard pW t q tPr0,τ s -Brownian motion and let x 0 : Ω Ñ R d be an F 0 {BpR d q-measurable mapping with Er}x 0 } p s ă 8 for all p P r1, 8q. Moreover, let f : R d Ñ R d be a continuously differentiable function whose derivative grows at most polynomially. Formally, let us suppose the existence of constants D P R ě0 and q P N such that, for all x, y P R d }f pxq ´f pyq} 2 ď D}x ´y} 2 p1 `}x} q `}y} q q (H1)

Let g " pg i,j q iPt1,...,du,jPt1,...,mu : R d Ñ R dˆm be a globally Lipschitz continuous function: there exists L g P R ě0 such that, for all x, y P R d }gpxq ´gpyq} ď L g }x ´y} (H2)

Finally, let us suppose that f is globally one-sided Lipschitz with constant λ P R:

Dλ P R @x, y P R d : xf pyq ´f pxq, y ´xy ď λ }y ´x} 2 (H3)

Then consider the Stochastic Differential Equations (SDE):

dX t " f pX t qdt `gpX t qdW t , X 0 " x 0 (1)
for t P r0, τ s. The drift coefficient f is the infinitesimal mean of the process X and the diffusion coefficient g is the infinitesimal standard deviation of the process X. Under the above assumptions, the SDE (1) is known to have a unique strong solution. More formally, there exists an adapted stochastic process pX t,x0 q : r0, τ s ˆΩ Ñ R d with continuous sample paths fulfilling

X t,x0 " x 0 `ż t 0 f pX s qds `ż t 0 gpX s qdW s
for all t P r0, τ s P-a.s. (see, e.g., ([21])). We denote by X t,x0 the solution of Equation ( 1) at time t from initial condition X 0,x0 " x 0 P-a.s., in which x 0 is a random variable that is measurable in F 0 .

Remark 1. Constants λ, L g and D can be computed using (constrained) optimization algorithms (see ( [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF])).

Overview of the approach

We aim at synthesizing reachability controllers for the class of systems described above. Our approach consists in adapting reachability analysis based control systhesis methods to stochastic systems. In [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF], we successfully used the Euler method for guaranteeing controlled reachability of deterministic systems using the Euler method in association to a new error bound relying on the OSL property of the vector field.

Consider equation (1) with g " 0. Classically, one knows that, if the function f is Lipschitz continuous with Lipschitz constant L, the solution of the ODE starting at a given initial value exists and is unique. Besides, one has:

}X t,x0 ´Xt,x1 } ď e Lt }x 0 ´x1 }, (2)
with two initial values x i (i " 0, 1). This gives a rough growth bound, i.e. a function bounding the distance of neighboring trajectories as t evolves. In [START_REF] Donchev | Stability and euler approximation of one-sided lipschitz differential inclusions[END_REF], it is proven that, if f is continuous and OSL with OSL constant λ, then the solution of the ODE starting at a given initial value exists and is unique, and, for all x 0 , x 1 P R n :

}X t,x0 ´Xt,x1 } ď e λt }x 0 ´x1 }. (3)
This gives a more accurate growth bound because a Lipschitz function f is always OSL, and the associated OSL constant λ is always less than or equal to its Lipschitz counterpart L. Using the OSL constant λ, it is also possible to bound the error }X t,x0 ´X t,x1 } in function of }x 0 ´x1 }, where Xt,x1 denotes the Euler approximate trajectory of the solution X t,x1 . We can then conclude that any trajectory starting from x 0 close enough from a given initial condition x 1 (}x 0 ´x1 } ď δ 0 ) will remain within some bound (denoted by δ t,δ0 ) of the Euler trajectory Xt,x1 . Our objective is to use a similar approach for stochastic systems, while making as few additional hypotheses on the dynamics as possible.

Therefore, we first explain how we can use a variant of the Euler-Maruyama scheme called tamed Euler scheme in order to guarantee reachability of the expectation of the state ErX t,x0 s of the stochastic system (1) in a given set. It consists in computing an approximate Euler-based trajectory Xt,x1 , we then provide an error bound for Er}X t,x0 ´X t,x1 } 2 s (Theorem 2) which allows to write a result for the expectation of the non approximated trajectory ErX t,x0 s (Proposition 2). This result can be extended further as a probability result (the probability with which the state reaches the set, Proposition 3).

This reachability analysis method is then used in association with tiling based control synthesis algorithms to ensure that the system is stabilised in a given region of interest with a given probability (Section 4). The first algorithm we use is based on an adaptive state-space decomposition algorithm [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF][START_REF] Le Coënt | An improved algorithm for the control synthesis of nonlinear sampled switched systems[END_REF], and uses an exhaustive enumeration of all possible control sequences (up to a given length) in order to find stabilising sequences (see Section 4.3). We then propose a faster approach (in Section 4.4) which relies on state-space discretization and dynamic programming techniques in order to minimize the average distance of the state with some given objective after a given number of steps. The second algorithm is merely an enhancement of the control sequence search of the first one, using a static tiling. Note that with both algorithms, there is no need to compute a finite-state abstraction of the system, which means that computing a controller is done with the application of a single algorithm, unlike [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF] which requires the computation of a finite-state abstraction before computing controllers.

Example 1. Throughout the paper, we illustrate our approach on a 2-mode switched system borrowed from [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF]. This case-study is a stochastic version of a well-known system originally introduced as an illustrative example in [START_REF] Liberzon | Switching in systems and control[END_REF]. The dynamics is given by the system of equations dx 1 " p´0.25x 1 `px 2 `p´1q p 0.25qdt `wx 1 dW 1 t , dx 2 " ppp ´3qx 1 ´0.25x 2 `p´1q p p3 ´pqqdt `wx 2 dW 2 t , where p " 1, 2 is the mode, and we consider different values of noise w P R ą0 .

Bounding the Error of the tamed Euler method

Tamed Euler approximation

The standard way to extend the classical Euler method for ordinary differential equations to the SDE (1) is the Euler-Maruyama scheme ( [START_REF] Maruyama | Continuous Markov processes and stochastic equations[END_REF]). More precisely, given z : Ω Ñ R d an F 0 {BpR d q-measurable mapping with Er}z} p s ă 8 for all p P r1, 8q, the explicit Euler-Maruyama (EM) method for the SDE (1) is given by the mappings

Y N n,z : Ω Ñ R d , n P t0, 1, . . . , N u, which satisfy Y N 0,z " z and Y N n`1,z " Y N n,z `τ N ¨f pY N n,z q `gpY N n,z qpW pn`1qτ {N ´Wnτ{N q
for all n P t0, 1, . . . , N ´1u and all N P N. See ( [START_REF] Maruyama | Continuous Markov processes and stochastic equations[END_REF]). Unfortunately, the convergence results for the EM scheme does not hold when the drift function f of the SDE (1) is polynomial (not an affine function). In order to consider non affine drift functions, we adopt a refined scheme, which has been proposed recently to overcome this difficulty ( [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]). Let

X N n,z : Ω Ñ R d , X N n`1,z " X N n,z `τ N ¨f pX N n,z q 1 `τ N ¨}f pX N n,z q} `gpX N n,z qpW pn`1qτ N ´W nτ N q (4)
for all n P t0, 1, . . . , N ´1u and all N P N. We refer to the numerical method (4) as a tamed Euler scheme ( [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]). In this method the drift term τ n ¨f pX N n,z q is "tamed" by the factor 1{p1 `τ N ¨}f pX N n,z q}q for n P t0, 1, . . . , N ´1u and N P N in (4). In order to establish results for the time continuous system, we use the following time continuous interpolation of the time discrete numerical approximations (4), also introduced in ([13]). Let XN z : r0, τ s ˆΩ Ñ R d , N P N, be a sequence of stochastic processes given by

XN t,z " XN n,z `pt ´nτ {N q ¨f p XN n,z q 1 `τ {N ¨}f p XN n,z q} `gp XN n,z qpW t ´W nτ N q (5)
for all t P r nτ N , pn`1qτ N s, n P t0, 1 . . . , N ´1u and all N P N. Note that XN t,z : r0, τ s ˆΩ Ñ R d is an adapted stochastic process with continuous sample paths for every N P N.

We finally introduce a piecewise constant interpolation X N t,z of the tamed Euler scheme as follows. It is mainly used in some constants when establishing the error bound and for proving Theorem 2.

X N t,z :" X N n,z for t P r nτ N , pn `1qτ N q. (6)
Note that XN t,z " X N t,z " X N n,z at time t " nτ N for n P t0, 1, . . . , N u. The following theorem is proven in ( [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]): Theorem 1. Let us suppose (H1)-(H2)-(H3). Let z : Ω Ñ R d be an F 0 {BpR d qmeasurable mapping with Er}z} p s ă 8 for all p P r1, 8q. Then, for all p P r1, 8q, the tamed Euler scheme (4) satisfies:

sup N PN sup nPt0,1,...,N u Er}X N n,z } p s ă 8
This theorem allows to ensure the strong convergence of the tamed Euler method. Any number N of subsampling steps can thus be used. This number is now left implicit for the sake of simplicity. From Theorem 1, it follows (cf. Lemma 4.3,[START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF])): Lemma 1. Let us suppose (H1)-(H2)-(H3). Let z : Ω Ñ R d be an F 0 {BpR d qmeasurable mapping with Er}z} p s ă 8 for all p P r1, 8q. Let us consider the time continuous interpolation of the tamed Euler scheme as defined in [START_REF] Davis | On the construction of a class of invariant polynomials in several matrices, extending the zonal polynomials[END_REF], and the piecewise constant interpolation as defined in [START_REF] Donchev | Stability and euler approximation of one-sided lipschitz differential inclusions[END_REF]. Then, for any even integer r ě 2, there exist two constants E r,z and F r,z such that

sup 0ďtďτ E}X t,z ´X t,z } r ď p∆ t q r 2 pE r,z p∆ t q r 2 `Fr,z dq. with ∆ t " τ {N and: E r,z " 2 r p}f p0q} r `D2 r`1 2 p1 `E sup 0ďtďτ }X t,z } qr q 1 2 pE sup 0ďtďτ }X t,z } 2r q 1 2 q, F r,z " 2 r p}gp0q} 2r `Lr g E sup 0ďtďτ }X t,z } r 2 q.
Remark 2. Constants E r,z and F r,z are computed using the constants λ and L g (see Remark 1), and the expected values of sup 0ďtďτ }X t,z } p for the required values of p at each time t " 0, ∆t, 2∆t, . . . , N ∆t. These expected values can either be computed numerically by using a Monte Carlo method, or by evaluating analytically the expectations E sup 0ďtďτ }X t,z } p for all required values of p. Proposition 2 explains how such expectations can be computed for the piecewise linear interpolation Xt,z , the computation for X p t,z is the same but much simpler since there is no time dependent term in [START_REF] Donchev | Stability and euler approximation of one-sided lipschitz differential inclusions[END_REF]. In this paper, we use the latter approach.

Mean square error bounding

The following Theorem holds for SDE [START_REF] Arnold | Stochastic differential equations[END_REF]. This corresponds to a stochastic version of Theorem 1 of ( [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]), showing that a similar result holds on average, using the tamed Euler method of ( [START_REF] Hutzenthaler | Strong convergence of an explicit numerical method for SDEs with nonglobally lipschitz continuous coefficients[END_REF]). It is an adaptation of Theorem 4.4 in ( [START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF]).

Theorem 2. Given the SDE system (1) satisfying (H1)-(H2)-(H3). Let δ 0 P R ě0 . Suppose that z is a random variable on R d such that Er}x 0 ´z} 2 s ď δ 2 0 .

Then, we have, for all τ ě 0:

Er sup 0ďtďτ }X t,x0 ´X t,z } 2 s ď δ 2 τ,δ0 ,
with δ 2 τ,δ0 :" βpτ qe γτ , where:

γ " 2p ? ∆ t `2λ `L2 g `128L 4
g q, and

βpτ q " 2δ 2 0 `2τ ∆ t L 2 g p1 `128L 2 g qpF 2,z d `E2,z ∆ t q `4τ a ∆ t DpF 4,z d `E4,z ∆ 2 t q 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 . ( 7
)
with ∆ t " τ {N .

Proof. The proof closely follows the proof of Theorem 4.4 in ( [START_REF] Higham | Strong convergence of Euler-type methods for nonlinear stochastic differential equations[END_REF]). Let e t " X t,x0 ´X t,z . We have, for all 0 ď t ď τ : de t " pf pX t,x0 q ´f pzqqdt `pgpX t,x0 q ´gpzqqdW t .

Then, by using Equation ( 8) and the integral version of Itô formula applied to function x Þ Ñ }x} 2 we obtain

}e t } 2 " }e 0 } 2 `ż t 0 2xe s , f pX s,x0 q ´f pX s,z qyds `ż t 0 }gpX s,x0 q ´gpX s,z q} 2 ds `M ptq, (9)
where e 0 " x 0 ´z, and M ptq "

ż t 0 2xe s , gpX s,x0 q ´gpX s,z qydW s .
So we have using (H2):

}e t } 2 ď }e 0 } 2 `ż t 0 2xe s , f pX s,x0 q ´f p Xs,z qyds `L2 g ż t 0 }X s,x0 ´Xs,z } 2 ds `ż t 0 2xe s , f p Xs,z q ´f pX s,z qyds `M ptq. (10)
So we have using (H3) and Young's inequality:

}e t } 2 ď }e 0 } 2 `ż t 0 p2λ}e s } 2 `L2 g }e s } 2 qds `L2 g ż t 0 }X s,z ´X s,z } 2 ds `ż t 0 p 1 ? ∆ t }f p Xs,z q ´f pX s,z q} 2 `a∆ t }e s } 2 qds `M ptq. (11)
So we have using (H1), for all 0 ď t ď τ :

}e t } 2 ď }e 0 } 2 `pa ∆ t `2λ `L2 g q ż t 0 }e s } 2 ds `L2 g ż t 0 }X s,z ´X s,z } 2 ds `D ? ∆ t ż t 0 p1 `}X s,z } q `} Xs,z } q q}X s,z ´X s,z } 2 ds `M ptq. (12)
It follows using Lemma 1 for r " 2, and Cauchy-Schwarz inequality:

Er sup 0ďsďt }e s } 2 s ď E}e 0 } 2 `pa ∆ t `2λ `L2 g q ż t 0 E}e s } 2 ds `L2 g τ ∆ t pE 2,z ∆ t `F2,z dq `D ? ∆ t ż t 0 pEp1 `}X s,z } q `} Xs,z } q q 2 q 1 2 pE}X s,z ´X s,z } 4 q 1 2 ds `mptq, (13)
where mptq " Er sup 0ďsďt }M psq}s.

Hence, using using Lemma 1 for r " 4, and inequality pa `bq p ď 2 p pa p `bp q:

Er sup 0ďsďt }e s } 2 s ď E}e 0 } 2 `pa ∆ t `2λ `L2 g qq ż t 0 E}e s } 2 ds `L2 g τ ∆ t pE 2,z ∆ t `F2,z dq `2Dτ a ∆ t pE 4,z ∆ 2 t `F4,z dq 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 `mptq. (14)
On the other hand, from the Burkholder-Davis-Gundy inequality, we get:

mptq ď 16Er ż t 0 }e s } 2 }gpX s,x0 q ´gpX s,z q} 2 dss 1 2
Hence, using (H2):

mptq ď 16L 2 g Er sup 0ďsďt }e s } 2 ż t 0 }X s,x0 ´Xs,z } 2 dss 1 2
Then, using Young's inequality (for any α ą 0):

mptq ď 8L 2 g pαEr sup 0ďsďt }e s } 2 s `1 α Er ż t 0 }X s,x0 ´Xs,z } 2 dssq.
Hence, by using Lemma 1 for r " 2:

mptq ď 8αL 2 g Er sup 0ďsďt }e s } 2 s `8L 2 g α ż t 0 Er sup 0ďrďs }e r } 2 sds `8L 2 g α τ ∆ t pE 2,z ∆ t `F2,z dq. (15)
Hence, letting α " 1 16L 2 g , we have by replacing in ( 14):

1 2 Er sup 0ďsďt }e s } 2 s ď δ 2 0 `pa ∆ t `2λ `L2 g `128L 4 g q ż t 0 Er sup 0ďrďs }e r } 2 sds `τ pL 2 g `128L 4 g q∆ t pE 2,z ∆ t `F2,z dq `τ 2D a ∆ t pE 4,z ∆ 2 t `F4,z dq 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 . ( 16
)
It results from Gronwall's inequality:

Er sup 0ďtďτ }e t } 2 s " βpτ qe γτ , with γ " 2p ? ∆ t `2λ `L2 g `128L 4
g q, and

βpτ q " 2δ 2 0 `2τ p∆ t L 2 g p1 `128L 2 g qpF 2,z d `E2,z ∆ t q `4τ a ∆ t DpF 4,z d `E4,z ∆ 2 t q 1 2 p1 `4E sup 0ďtďτ }X t,z } 2q `4E sup 0ďtďτ } Xt,z } 2q q 1 2 . (17)
It follows from Theorem 2 and Jensen's inequality:

Proposition 1. Consider two points x 0 and z in R d ,and a positive real number δ 0 . Suppose that x 0 P Bpz, δ 0 q. Then EX t,x0 P Bp Xt,z , δ t,δ0 q for all t P r0, τ s, where δ t,δ0 is defined in (7).

Numerical pre-computation of expectations

Theorem 2 requires the knowledge of E sup 0ďtďτ } Xt,z } 2q . In this section, we give a pre-computable upper bound of this quantity.

Proposition 2. Let Xt,z be the piecewise linear interpolation of the solution to equation (1) as defined in equation [START_REF] Davis | On the construction of a class of invariant polynomials in several matrices, extending the zonal polynomials[END_REF], written in the form Xt,z " αpzq βpzqt

`γpzqW t with αpzq " XN n,z ´pnτ{Nq¨fp XN n,z q 1`τ {N ¨}f p XN n,z q} , βpzq " f p XN n,z q 1`τ {N ¨}f p XN n,z q} , γpzq " gp XN n,z q. Then Er} Xt,z } 2q s " ÿ |k|"q ˆq k ˙Ck1,k2,k3 pzqµ k4,k5,k6 pz, 1qt k2`2k3`k4`k 4 `k5 2 `k6 , ( 18
)
where k " pk 1 , . . . , k 6 q. The previous equation can be written as follows

} Xt,z } 2q " ÿ |k|"q
ˆq k ˙Ck1,k2,k3 pzqt k2`2k3`k4 papzq T W t q k4 pbpzq T W t q k5 pW T t QpzqW t q k6

where apzq " 2βpzq T γpzq, bpzq " 2αpzq T γpzq, Qpzq " γpzq T γpzq and C k1,k2,k3 pzq " ´αpzq T αpzq ¯k1 ´2αpzq T βpzq ¯k2 ´βpzq T βpzq ¯k3 .

We denote µ κ pz, tq def.

" E " papzq T W t q k4 pbpzq T W t q k5 pW T t QpzqW t q k6 ‰ where κ " pk 4 , k 5 , k 6 q. The expectation µ κ pz, tq can always be written as the high-order moment of a product of three quadratic forms [START_REF]The moments of products of quadratic forms in normal variables[END_REF], therefore we have

µ κ pz, tq " $ ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' % λ p1q κ D k 5 ´k4 2 ,k4,k6 pbpzqbpzq T , P a,b pzq, Qpzqqt k 4 `k5 2 `k6 if k 4 ă k 5 and k 4 `k5 " 0 pmod 2q λ p2q κ D k 4 ´k5 2 ,k5,k6 papzqapzq T , P a,b pzq, Qpzqqt k 4 `k5 2 `k6 if k 4 ą k 5 and k 4 `k5 " 0 pmod 2q λ p3q κ D k5,k6 pP a,b pzq, Qpzqqt k5`k6
if k 4 " k 5 and k 5 " 0 pmod 2q 0 if k 4 `k5 " 1 pmod 2q where D n1,...,ni pA 1 , . . . , A i q denotes the normalized Davis polynomials [START_REF] Davis | On the construction of a class of invariant polynomials in several matrices, extending the zonal polynomials[END_REF][START_REF] Hillier | Generating functions and short recursions, with applications to the moments of quadratic forms in noncentral normal vectors[END_REF] and

P a,b pzq " apzqbpzq T `bpzqapzq T 2 , λ p1q κ " 2 k 4 `k5 2
`k6´1 k 4 k 6 pk 5 ´k4 q,

λ p2q κ " 2 k 4 `k5 2
`k6´1 k 5 k 6 pk 4 ´k5 q, λ p3q κ " 2 k5`k6 k 5 k 6 . Finally, the result holds.

The normalized Davis polynomials D n1,...,ni pA 1 , . . . , A i q can be computed efficiently using recursion formulas, see [START_REF] Hillier | Generating functions and short recursions, with applications to the moments of quadratic forms in noncentral normal vectors[END_REF].

From Proposition 2, we obtain an upper bound of E sup 0ďtďτ } Xt,z } 2q using the following lemma. Lemma 2. For every martingale Y t P L p (where p ą 1), we have:

Er sup 0ďtďτ }Y t } p s ď ˆp p ´1 ˙p Er}Y τ } p s.
Proof. See [START_REF] Arnold | Stochastic differential equations[END_REF].

We thus have:

E sup 0ďtďτ } Xt,z } 2q ď ˆ2q 2q ´1 ˙2q ÿ |k|"q ˆq k ˙Ck1,k2,k3 pzqµ k4,k5,k6 pz, 1qt k2`2k3`k4`k 4 `k5 2 `k6
where C k1,k2,k3 pzq and µ k4,k5,k6 pz, 1q are given in Proposition 2. The computation of this deterministic upper-bound is much faster than the direct Monte-Carlo approximation of the expectation. Furthermore, we get rid of the approximations inherent to the use of Monte-Carlo simulations.

Probabilistic reachability

Let us suppose for simplicity that R is a ball of the form BpO, ρq where O is the origin and ρ ą 0. Suppose that x 0 P Bpz, δ 0 q. We have:

P rpX τ,x0 R Rq " P rp}X τ,x0 } ą ρq. Furthermore: Proposition 3. P rpX τ,x0 P Rq ě 1 ´1 ρ 2 ˆδτ,δ0 `bEr} Xτ,z } 2 s ˙2 .
Proof. Using Chebyshev's inequality and triangular inequality, it follows

P rpX τ,x0 R Rq ď 1 ρ 2 Er}X τ,x0 } 2 s ď 1 ρ 2 ˆbEr} Xτ,z ´Xτ,x0 } 2 s `bEr} Xτ,z } 2 s ˙2 .
We thus get, using the inequality of the previous section (namely, Er} Xτ,z Xτ,x0 } 2 ď δ 2 τ,δ0 ):

P rpX τ,x0 R Rq ď 1 ρ 2 pδ τ,δ0 q 2 `bEr} Xτ,z } 2 s.
Alternatively, we have:

P rpX τ,x0 P Rq ě 1 ´1 ρ 2 ˆδτ,δ0 `bEr} Xτ,z } 2 s ˙2 .
NB: If we have R " BpO 1 , ρq where O 1 is not the origin, we have:

P rpX τ,x0 P Rq ě 1 ´1 ρ 2 ˆδτ,δ0 `bEr} Xτ,z ´O1 } 2 s ˙2 .

Implementation

This method has been implemented in the interpreted language Octave, and the experiments performed on a 2.80 GHz Intel Core i7-4810MQ CPU with 8 GB of memory. The implementation is an adaptation of the program described in ( [START_REF] Le Coënt | Control synthesis of nonlinear sampled switched systems using Euler's method[END_REF]) for controlling deterministic switched systems, but makes use of the tamed Euler scheme for SDEs (with the error function δ t,δ0 given in Theorem 2) instead of the classical Euler scheme.

Example 2. Consider the system of Example 1, for mode u " 1 and w " 0.05:

dx 1 " p´0.25x 1 `x2 `0.25qdt `0.05x 1 dW 1 t dx 2 " p´2x 1 ´0.25x 2 ´2qdt `0.05x 2 dW 2 t
The program gives (for τ " 1, ∆ t " τ {10 4 ): q " 0, D " 1.36, L g " 0.05, λ " 0.25; and for z " p´4, ´3.8q: E 2,z " 893.3, E 4,z " 2.14 ¨10 5 , F 2,z " 0.002, F 4,z " 4.9 ¨10 ´6. Consider now the system corresponding to Example 1 for mode u " 2 and w " 0.05:

dx 1 " p´0.25x 1 `2x 2 ´0.25qdt `0.05x 1 dW 1 t dx 2 " p´x 1 ´0.25x 2 `1qdt `0.05x 2 dW 2 t
The program gives (for τ " 1, ∆ t " τ {10 4 ): q " 0, D " 1.36, L g " 0.05, λ " 0.25, and, for z " p0, 3q: E 2,z " 543.2, E 4,z " 7.94 ¨10 4 , F 2,z " 0.0442, F 4,z " 0.00178. Both computations take less than 10s of CPU time. Simulations of the two systems are given in Figure 1 for mode u " 1 and starting point z " p´4, 3.8q, and mode u " 2 and starting point z " p0, 3q. The initial ball Bpz, δ 0 q is depicted in black, the final ball BpEr Xτ,z s, δ τ,δ0 q in red, and 200 random sampling trajectories in blue for t P r0, τ s.

Example 3. We now consider a nonlinear model of a pendulum on a cart borrowed from [START_REF] Borri | Symbolic models for nonlinear control systems affected by disturbances[END_REF][START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF]. The dynamics is given by:

dx 1 " x 2 dt `0.03x 1 dW 1 t dx 2 " p´g l sin x 1 ´k m x 2 `1 ml 2 uqdt `0.03x 2 dW 2 t
where the state px 1 , x 2 q represents the angular position and velocity of the point mass on the pendulum, u is the control input (torque applied to the cart). The parameters are g " 9.8, l " 0.5, m " 0.6, k " 2. The program gives for control input u " 0.5, τ " 0.1, ∆ t " τ {10 2 : q " 1, D " 2.56, L g " 0.03, λ " 1.446; and for z " p0, 0q:

E 2,z " 47.5, E 4,z " 1977.2 ¨10 5 , F 2,z " 9.8 ¨10 ´4, F 4,z " 9.73 ¨10 ´7.
The computation takes less than 10s of CPU time. Simulations of the system are given in Figure 2 for control input u " 0.5 and starting point z " p0, 0q. The initial ball Bpz, δ 0 q is depicted in black, the final ball BpEr Xτ,z s, δ τ,δ0 q in red, and 200 random sampling trajectories in blue for t P r0, τ s.

Sampled stochastic switched systems

Stochastic switched system as a finite collection of SDEs

We now consider a finite number of SDEs. Each SDE is referred to as a mode j, and the set of modes is referred to as U " t1, . . . , M u. We will denote by X j t,x0 the solution at time t of the system:

dxptq " f j pxptqq `gj pxptqqdW j t , xp0q " x 0 .

(

) 19
where x 0 is a random variable that is measurable in F 0 . Hypotheses (H1)-(H2)-(H3), as defined in Section 3, are naturally extended to every mode j of U . Accordingly, constants L g , λ, F associated to SDE (1) in Section 3, now become L gj , λ j , F j respectively, for each j P U . Likewise, for each j P U , the nonnegative real pδ t,δ0 q 2 becomes pδ j t,δ0 q 2 for each mode j; the approximate continuous-time solution of [START_REF] Maruyama | Continuous Markov processes and stochastic equations[END_REF] starting from z, is denoted by Xj t,z , and the approximate staircase solution by X j t,z .

Control patterns

The control laws that we now consider are "piecewise constant of duration τ " in the sense that, every τ seconds, they select a given mode (see ( [START_REF] Zamani | Symbolic models for stochastic switched systems: A discretization and a discretization-free approach[END_REF])). We call "(control) pattern of length k" a sequence of k modes (i.e., an element of U k ). Each pattern π of the form j 1 j 2 ¨¨¨j k corresponds to the selection of mode j 1 for time t P r0, τ q, then mode j 2 for t P rτ, 2τ q, and so on, until t " kτ . We assume that the solution of the system is continuous at sampling instants t " τ, 2τ, . . . (which means that there is no "reset" of the system at sampling instants).

Given a stochastic switched system, a pattern π of length k and an initial random variable z, one constructs the "approximate solution controlled by π" by composing together the approximations obtained by successive application of the modes of the pattern π. Formally, the "continuous" approximate solution Xπ t,z is defined at time t P r0, kτ s as follows:

• Xπ t,z " Xj t,z if π " j P U , k " 1 and t P r0, τ s, and

• Xπ pk´1qτ `t1 ,z " Xj t,z 1 with z 1 " Xπ 1 pk´1qτ,z if k ě 2, t 1 P r0, τ s, π " π 1 ˚j for some j P U and π 1 P U k´1 .

The "staircase" approximate solution X π t,z is defined analogously. Likewise, given an initial error radius δ 0 ą 0 and a pattern π of length k ě 1, one defines the error radius δ π t,δ0 using (7) as follows:

• δ π t,δ0 " δ j t,δ0 if π " j P U , k " 1 and t P r0, τ s, and

• δ π pk´1qτ `t1 ,δ0 " δ j t 1 ,δ 1 with δ 1 " δ π 1 pk´1qτ,δ0
, if k ě 2, t 1 P r0, τ s, π " π 1 ˚j for some j P U and π 1 P U k´1 .

Controlled probabilistic reachability

Given a ball R " BpO, ρq Ă R d with ρ ą 0, we define the problem of "controlled probabilistic reachability inside R" as follows:

Given a starting point z P R d and a positive real δ 0 , find a pattern π of length k such that, at t " kτ , the image of the points located initially in Bpz, δ 0 q be, on average, located as close as possible to O, and, in particular:

BpEr Xπ t,z s, δ π t,δ0 q Ď R for t " kτ.
The program mentioned in Section 3.5, has been extended in order to find, by exhaustive enumeration1 , such a pattern π. This guarantees that, for all x 0 P Bpz, δ 0 q:

P rpX π kτ,x0 P Rq ě 1 ´1 ρ 2 p b Er} Xπ kτ,z } 2 s `δπ kτ,δ0 q 2 ,
We now give two examples of application of this program.

Example 4. Consider the system of Example 1 with w " 0.01:

dx 1 " p´0.25x 1 `ux 2 `p´1q u 0.25qdt `0.01x 1 dW 1 t dx 2 " ppu ´3qx 1 ´0.25x 2 `p´1q u p3 ´uqqdt `0.01x 2 dW 2 t
where u " 1, 2.

For τ " 0.5, ∆ t " 10 ´4, one finds (for all modes u " 1, 2): q " 0, D " 1.36, L g " 0.01, λ " 0.25; for z " p´4, ´3.8q: E 2,z " 893.31, E 4,z " 2.14 ¨10 5 , F 2,z " 0.002, F 4,z " 4.9 ¨10 ´6; and for z " p0, 3q: E 2,z " 543.22, E 4,z " 7.94 ¨10 4 , F 2,z " 0.0442, F 4,z " 0.00178. Our program shows the probabilistic reachability inside R " Bpp0, 0q, ρq, for all point x 0 P Bpp5, 4q, δ 0 q and all point y 0 P Bpp´5, ´4q, δ 0 q, with δ 0 " 0.1 and ρ " 7. We have, for π 1 " p2 ¨2 ¨2 ¨2q

P rpX π1 4τ,x0 P Rq ě 1 ´1 ρ 2 p b Er} Xπ1 4τ,p5,4q } 2 s `δπ1 4τ,δ0 q 2 ,
P rpX π1 4τ,x0 P Rq ě 0.689, and for π 2 " p1 ¨1 ¨1 ¨1 ¨1q:

P rpX π2 5τ,y0 P Rq ě 1 ´1 ρ 2 p b Er} Xπ2 5τ,p´5,´4q } 2 s `δπ2 5τ,δ0 q 2 .
P rpX π2 5τ,y0 P Rq ě 0.713 Figure 3 depicts in black the initial balls (at t " 0) centered at p5, 4q and p´5, ´4q; and for each initial ball, the pattern that sends the ball inside R (at time t " kτ ); the intermediate balls (at t " τ, 2τ, . . . , pk ´1qτ ) are depicted in red, and 200 sampling trajectories drawn in black.

Example 5. (the slit problem)

The problem is adapted from ( [START_REF] Morzfeld | Implicit sampling for path integral control, Monte Carlo localization, and SLAM[END_REF]). The controlled dynamics is:

dX " udt `dW, X 0 " 1 with mode u P t´6, ´5, ´4, ´3, ´2, 1, 0, 1, 2, 3, 4, 5, 6u. We have (at t " 0.5) a slit at x P r´1, ´3s. The objective is thus to control the system so that xptq P R " r´1, ´3s (i.e., xptq P R " BpO 1 , ρq with O 1 " ´2 and ρ " 1q at t " 0.5. One has, for all modes: q " 0, D " 0, L g " 0, λ " 0. For δ 0 " 0.2, an initial point z " 1 and a sampling time τ " 0.5 with subsampling ∆ t " 2.5 ˆ10 ´3, one has for mode u " ´6: E 2,z " 144, E 4,z " 20736, F 2,z " 4, F 4,z " 16; and for mode u " 0: E 2,z " 0, E 4,z " 0, F 2,z " 4, F 4,z " 16.

Suppose that all the trajectories start at x 0 P Bpz, δ 0 q, with z " 1 and δ 0 " 0.2. When there is no control (u " 0), at time t " 0.5, the expected value of X t,x0 is in Bpz 1 , δ t,δ0 q with z 1 " 1 and δ t,δ0 " 0.29. From Markov's inequality, it follows that the trajectories pass by R " Bp´2, 1q at t " 0.5 with low probability: see Figure 4. On the other hand, with control u " ´6, at time t " τ " 0.5, the expected value of X t,x0 is now in Bpz 1 1 , δ 1 τ,δ0 q with z 1 1 " ´2 and δ 1 τ,δ0 " 0.29. This explains why the trajectories pass by R " Bp´2, 1q at t " 0.5 with high probability. We have, for x 0 P Bpz, δ 0 q with z " 1, δ 0 " 0.2, t " τ " 0.5, ρ " 1, O 1 " ´2:

P rpX u τ,x0 P R " BpO 1 , ρqq ě 1 ´1 ρ 2 p b Er} Xu τ,z ´O1 } 2 s `δu τ,δ0 q 2 . ( 20
)
Figure 4: Simulations of Example 5 without control (π " p0 ¨0q) and with control pattern p´6 ¨0q from the initial ball Bp1, 0.2q. The initial region is in green. The red regions are BpEr Xπ t,z s, δ π t,δ 0 q. with z " 1, δ 0 " 0.2 for the two patterns π " p0 ¨0q (top) and π " p´6 ¨0q (bottom).

For mode u " ´6, this gives: δ u τ,δ0 " 0.29 and Er} Xu τ,z ´O1 } 2 s " 0.048. It follows: P rpX u τ,x0 P Rq ě 0.871.

Optimal control by state space discretization

In this Subsection, we are interested in solving an optimal control problem in order to minimize the average distance of the state reached after a given number of steps to the origin. In order to solve such optimal control problems, it is classical to spatially discretize the set S into a finite number of cells of equal size. Consider a partition of a compact subset S of R d in cells. For the sake of simplicity, we suppose that each cell C is a hypercube of R d of side length 2ε{ ? d, and we consider the centre of the cell, denoted 7x, as the "representative" of the cell. Note that any point z located in the cell is such that }z ´7x} ď ε. In [START_REF] Bertsekas | Convergence of discretization procedures in dynamic programming[END_REF] (cf. [START_REF] Esmaeil | Adaptive and sequential gridding procedures for the abstraction and verification of stochastic processes[END_REF]), it is explained how a discretized procedure, based on classical techniques of dynamic programming, constructs a mapping u which associates to each representative 7x a mode up7xq P U , corresponding to the optimal policy at position 7x. More precisely, for any representative point 7x 0,ε , the discretized procedure P ROCp7x 0,ε q constructs a sequence of modes u 1,ε , . . . , u k,ε and a sequence of representatives 7x 1,ε , . . . , 7x k,ε such that, for all 1 ď i ď k:

• u i,ε " up7x i´1,ε q.
• 7x i,ε is the representative of the cell containing Er Xui,ε τ,7xi´1,ε s. For the sake of notation simplicity, we will now drop the index ε in symbols u i,ε and 7x i,ε . Let us define the function next ui by: next ui p7x i´1 q " 7x i , where 7x i is the representative of the cell containing Er Xui τ,7xi´1 s. We define the function next π as follows:

• if π " u for some u P U then next π p7xq " next u p7xq,

• if π is of the form u ¨π1 then next π p7xq " next π 1 pnext u p7xqq.

Given a point z 0 P R d and a positive integer k, we focus on the problem of finding the optimal pattern π ˚pz 0 q P U k , which, at time T " kτ , maps z 0 to a distribution located at a minimal distance of the origin, i.e., such that:

π ˚pz 0 q " argmin πPU k Er}X π T,z0 } 2 s.
Let 7x 0 be the representative of z 0 . The sequences tu i u 1ďiďk and t7x i u 1ďiďk given by P ROCp7x 0 q, are optimal in the following sense: (see, e.g., [START_REF] Bertsekas | Convergence of discretization procedures in dynamic programming[END_REF][START_REF] Esmaeil | Adaptive and sequential gridding procedures for the abstraction and verification of stochastic processes[END_REF]):

• the pattern u 1 ¨¨¨u k coincides with π ˚pz 0 q when ε is sufficiently small, and

• the discrete random variable 7x k converges in distribution to the random variable X u1¨¨¨u k T,z0

when ε tends to 0.

This justifies the employment of the pattern u 1 ¨¨¨u k generated by P ROCp7x 0 q as a (sub)optimal control pattern π for the original stochastic diffusion system2 .

Let us now, in our context, give a precise upper bound to Er}X u1¨¨¨u k T,z0

´7x k } 2 s, where 7x k is the representative point defined by next u1¨¨¨u k p7x 0 q. We first define λ i ě 0, for 1 ď i ď k:

λ i :" b Er} Xui τ,7xi´1 ´7x i } 2 s.
Let us also define ∆pu 1 ¨¨¨u k q recursively by:

• ∆pu 1 ¨¨¨u i q " δ ui τ,λi´1`∆pu1¨¨¨ui´1q for i " 2, . . . , k, and

• ∆pu 1 q " δ u1 τ,ε .

Without loss of understanding, we abbreviate ∆pu 1 ¨¨¨u i q by ∆ i . We have

Lemma 3. b Er}X u1¨¨¨u k kτ,z0 ´X u k τ,7x k´1 } 2 s ď ∆ k .
Proof. The proof is by induction on k. The case k " 1 is trivial. Suppose k ě 2.

We have:

b Er}X u1¨¨¨u k kτ,z0 ´X u k τ,7x k´1 } 2 s ď δ u k τ, b Er}X u 1 ¨¨¨u k´1 pk´1qτ,z 0 ´7x k´1 } 2 s
On the other hand, using the triangular inequality and the induction hypothesis, we have:

b Er}X

u1¨¨¨u k´1 pk´1qτ,z0 ´7x k´1 } 2 s ď b Er}X u1¨¨¨u k´1 pk´1qτ,z0 ´X u k´1 τ,7x k´2 } 2 s `bEr} Xu k´1 τ,7x k´2 ´7x k´1 } 2 s ď ∆ k´1 `λk´1 . Hence: b Er}X u1¨¨¨u k kτ,z0 ´X u k τ,7x k´1 } 2 s ď δ u k τ,∆ k´1 `λk´1 " ∆ k .
This achieves the proof of the induction step.

We have Theorem 3. Consider a point z 0 P R d of representative 7x 0 , and let 7x k " next u1¨¨¨u k p7x 0 q for sequence u 1 ¨¨¨u k P U k (generated by P ROCp7x 0 q). We have:

b

Er}X u1¨¨¨u k T,z0 ´7x k } 2 s ď ∆ k `λk .
Proof. We have

´7x k } 2 s ď ∆ k `λk .
Note that the term ∆ k of the inequality right-hand side takes into account the error introduced by the Euler scheme (with respect to the real solution of the SDE); in contrast, in the literature (e.g. [START_REF] Esmaeil | Adaptive and sequential gridding procedures for the abstraction and verification of stochastic processes[END_REF][START_REF] Pham | On some recent aspects of stochastic control and their applications[END_REF]), the Euler error is often disregarded, the analysis focusing just on the convergence of 7x k to the Euler approximation Xu1¨¨¨u k T,z0 . Example 6. We apply the above procedure P ROCp7x 0 q on Example 1 with w " 0.01. Consider the interest set R " Bpp0, 0q, ρq with ρ " 7, discretized with an accuracy ε " 0.57. We compute (sub)optimal patterns for the entire set, using different lengths of patterns, and simulate the induced controller for 200 initial conditions randomly selected in R. Simulations are given in Figure 5. The procedure took 11.8 seconds of computation for patterns of length 1, 47.4 seconds of computation for patterns of length 3, 152.0 seconds of computation for patterns of length 5. We now compare the results yielded by P ROCp7x 0 q with those of Example 4. The sub-optimal procedure returns patterns ensuring the probabilistic reachability inside R " Bpp0, 0q, ρq, for all point x 0 P Bpp5, 4q, δ 0 q and all point y 0 P Bpp´5, ´4q, δ 0 q, with δ 0 " 0.1 and ρ " 7. For x 0 P Bpp5, 4q, δ 0 q, the procedure returns pattern π 1 " p1 ¨2 ¨2 ¨2q ensuring a probability P rpX π1 4τ,x0 P Rq ě 0.801, and for x 0 P Bpp´5, ´4q, δ 0 q, the procedure returns pattern π 2 " p2 ¨2 ¨2 ¨1 ¨1q, with probability P rpX π2 5τ,y0 P Rq ě 0.802. The probabilities are thus higher than in Example 4. Simulations are shown in Figure 6, depicting in black the initial balls (at t " 0) centered at p5, 4q and p´5, ´4q; and for each initial ball, the pattern that sends the ball inside R (at time t " kτ ); the intermediate balls (at t " τ, 2τ, . . . , pk ´1qτ ) are depicted in red, and 200 sampling trajectories drawn in black. We consider the set of control inputs u P t´1.5, ´1, ´0.5, 0, 0.5, 1, 1.5u. In order to reproduce the results of [START_REF] Zamani | Symbolic control of stochastic systems via approximately bisimilar finite abstractions[END_REF], we slightly modify the procedure so that the interest set is R " Bpp0, 0q, ρq with ρ " 1.0, but the optimization objective is to minimize the average distance with the target state x obj " p0.5, 0q. The state space is discretized with an accuracy ε " 0.02. We compute (sub)optimal patterns for the entire set R, using patterns of length 1, and simulate the induced controller for 30 initial conditions randomly selected in R. Simulations are given in Figure 7. The procedure took 127 seconds of computation for patterns of length 1.

Final remarks and future work

We have established a mean-square error bound on an Euler-like method to simulate stochastic differential equations that satisfy a one-sided Lipschitz (OSL) condition. An advantage of our method is that the OSL assumption is weaker than the condition under which classical methods converge. From the mean square error bound, we have inferred a lower bound on the probability to reach a given set R at a given time. This allows us to determine a control pattern (in term of sequences of control modes) which makes the stochastic switched system reach a target with guaranteed minimum probability. Besides, instead of searching over all possible sequences of modes (which is a process exponential in the length of sequences), we have proposed an optimized method inferring a 'best' pattern, which may entail a significant gain in time and space. In the future, we plan to experiment our method with significantly larger examples than those given here.

Figure 1 :

 1 Figure 1: Simulations of Example 2 with mode u " 1 and initial ball Bpp´4, 3.8q, 0.5q, and mode u " 2 and initial ball Bpp0, 3q, 0.5q; τ " 1.

Figure 2 :

 2 Figure 2: Simulations of Example 3 with control input u " 0.5, initial ball Bpp0, 0q, 0.01q, τ " 0.1.

Figure 3 :

 3 Figure 3: Simulations of Example 4 from the initial balls Bpp5, 4q, 0.1q and Bpp´5, ´4q, 0.1q using patterns p2 ¨2 ¨2 ¨2q and p1 ¨1 ¨1 ¨1 ¨1q resp., with τ " 0.5.

Figure 5 :

 5 Figure 5: Simulations of Example 6 with the controller induced by P ROCp7x 0 q, for patterns of length 1 (top), length 3 (center), and length 5 (bottom). The blue circle is the set R " Bpp0, 0q, 7q, the red marker is the target state (the origin), the black lines are the controlled trajectories.

Figure 6 :Example 7 . 1 t dx 2 "

 6712 Figure 6: Simulations of Example 6 from the initial balls Bpp5, 4q, 0.1q and Bpp´5, ´4q, 0.1q using patterns p1 ¨2 ¨2 ¨2q and p2 ¨2 ¨2 ¨1 ¨1q resp., with τ " 0.5.

Figure 7 :

 7 Figure 7: Simulations of Example 7 projected on the first dimension, with the controller induced by P ROCp7x 0 q for patterns of length 1. The red line is the target state x obj " p0.5, 0q, the black lines are the controlled trajectories.

In Section 4.4 a strategy based on the search of an optimal pattern is suggested rather than an exhaustive search.

see, e.g.,[START_REF] Bertsekas | Convergence of discretization procedures in dynamic programming[END_REF]

This work is supported by the LASSO project financed by an ERC adv. grant; and the DiCyPS project funded by Innovation Fund Denmark.

License CC BY-NC-SA 4.0 : https://creativecommons.org/licenses/by-nc-sa/4.0/

where k " pk 1 , . . . , k 6 q and C k1,k2,k3 pzq " ´αpzq T αpzq ¯k1 ´2αpzq T βpzq ¯k2 ´βpzq T βpzq ¯k3 , and

where apzq " 2βpzq T γpzq, bpzq " 2αpzq T γpzq, Qpzq " γpzq T γpzq and D n1,...,ni pA 1 , . . . , A i q being the normalized Davis polynomials [START_REF] Davis | On the construction of a class of invariant polynomials in several matrices, extending the zonal polynomials[END_REF][START_REF] Hillier | Generating functions and short recursions, with applications to the moments of quadratic forms in noncentral normal vectors[END_REF] and

`k6´1 k 5 k 6 pk 4 ´k5 q, λ p3q κ " 2 k5`k6 k 5 k 6 . . Proof. First, note that we can write Xt,z as the sum of a mean value, a linear drift and a standard Wiener process Xt,z " αpzq `βpzqt `γpzqW t , for z P R d . Then, the squared norm of Xt,z is XT t,z Xt,z " pαpzq `βpzqt `γpzqW t q T pαpzq `βpzqt `γpzqW t q " αpzq T αpzq `2pαpzq T βpzq `βpzq T γpzqW t qt `βpzq T βpzqt 2 `2αpzq T γpzqW t `W T t γpzq T γpzqW t .

The Newton multinomial formula allows to express the 2q th power of the Xt,z norm

ˆq k ˙´αpzq T αpzq ¯k1 ´2αpzq T βpzqt ¯k2 ´βpzq T βpzqt 2 ¯k3 ´2βpzq T γpzqW t t ¯k4 ´2αpzq T γpzqW t ¯k5 ´W T t γpzq T γpzqW t ¯k6


Maud Charles
email: maud.charles@outlook.fr

Antonio Villalba

Gary Meyer

Suzanne Trancart

Coralie Lagy

Ismaël Bernard

Maryline Houssin

Francisella halioticida in mussels

Keywords: Francisella halioticida, mussels, granulocytomas, polymerase chain reaction, in situ hybridisation

This note describes the first detection of the bacteria Francisella halioticida in mussels (Mytilus spp.) from locations in Normandy and northern Brittany (France) experiencing high mussel mortalities, while it was not detected in Bay of St Brieuc (N Brittany), which was not affected by abnormal mussel mortality. The distribution of the bacteria in mussels seems to be restricted within inflammatory granulomas as observed in Yesso scallops Mizuhopecten yessoensis from Canada and Japan.

Francisella halioticida has been identified as being involved in mass (>80%) mortality of abalones Haliotis gigantea in Japan and high (up to 40%) mortality of Yesso scallops M. yessoensis in Canada as well as in lesions reducing marketability of Yesso scallops in Japan. The impact of this bacterium on the health of mussels needs to be investigated in future research, especially since the cause of high mussel mortalities that have been occurring in France for the past few years is still undetermined.

INTRODUCTION

The gram-negative bacteria Francisella halioticida was originally described in farmed abalone (Haliotis gigantea) as the cause of mass mortality that occurred in Japan [START_REF] Kamaishi | Mass mortality of giant abalone Haliotis gigantea caused by a Francisella sp. bacterium[END_REF], Brevik et al. 2011). Subsequently, this bacterium was identified in diseased Yesso scallops, Mizuhopecten (=Patinopecten) yessoensis, in conjunction with mortality events in Canada [START_REF] Meyer | Disease and mortality among Yesso scallops Patinopecten yessoensis putatively caused by infection with Francisella halioticida[END_REF] and Japan [START_REF] Kawahara | Francisella halioticida, Identified as the Most Probable Cause of Adductor Muscle Lesions in Yesso scallops Patinopecten yessoensis Cultured in Southern Hokkaido, Japan[END_REF]). In addition, [START_REF] Kawahara | Parallel studies confirm Francisella halioticida causes mortality in Yesso scallops Patinopecten yessoensis[END_REF] demonstrated that this bacterium can be highly virulent in Yesso scallops. Thus, F. halioticida has shown pathogenic potential in both gastropod and bivalve species of molluscs.

Since 2014, unexplained mass mortalities of mussels (Mytilus edulis, Mytilus galloprovincialis and their hybrids) have occurred in France from the Atlantic to the English Channel coasts [START_REF] Béchemin | Surmortalités de la moule bleue Mytilus edulis dans les Pertuis Charentais (mars 2014) -Rapport d'expertise[END_REF], 2015, Allain & Bernard 2016, Travers et al. 2016[START_REF] Bernard | Bilan de l'observatoire des mortalités de moules en Bretagne Nord pour la saison 2016-2017 et premiers éléments sur les organismes pathogènes présents -Rapport technique[END_REF], with mortality rates ranging from 40 % to 70 % [START_REF] Bernard | Bilan de l'observatoire des mortalités de moules en Bretagne Nord pour la saison 2016-2017 et premiers éléments sur les organismes pathogènes présents -Rapport technique[END_REF][START_REF] Bernard | Bilan de l'observatoire des mortalités de moules en Bretagne Nord pour la saison 2016-2017 et premiers éléments sur les organismes pathogènes présents -Rapport technique[END_REF][START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF]). Several hypotheses have been proposed to explain these mortalities, including the involvement of pathogenic bacteria belonging to the Vibrio splendidus clade [START_REF] Béchemin | Surmortalités de la moule bleue Mytilus edulis dans les Pertuis Charentais (mars 2014) -Rapport d'expertise[END_REF], 2015[START_REF] François | Bilan 2014 du réseau Repamo -Réseau national de surveillance de la santé des mollusques marins[END_REF], Soletchnik & Robert 2016, Travers et al. 2016), possibility of genomic abnormalities and disseminated neoplasia [START_REF] Benabdelmouna | The mass mortality of blue mussels (Mytilus spp.) from the Atlantic coast of France is associated with heavy genomic abnormalities as evidenced by flow cytometry[END_REF], Benabdelmouna et al. 2018,). Nevertheless, none of these proposals alone could explain such high mortality rates [START_REF] Charles | Étude des organismes pathogènes, des conditions physiologiques et pathologiques impliqués dans les mortalités anormales de moules (Mytilus sp.)[END_REF]). High mortality rates were recently observed among mussels in Northern Brittany and although V. splendidus strains and disseminated neoplasia were detected, the prevalence in mortality-affected versus non-affected areas did not correlate with mortality and no causative agent was identified [START_REF] Bernard | Bilan de l'observatoire des mortalités de moules en Bretagne Nord pour la saison 2016-2017 et premiers éléments sur les organismes pathogènes présents -Rapport technique[END_REF][START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF].

In a multi-parametric field study to identify potential causes of high mortality of mussels, [START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF] reported the occurrence of heavy haemocyte infiltra-tion and formation of granulomatous lesions in the connective tissue of various organs, however were not able to identify an aetiological agent. Similarities between the histological observations in mussels and the histopathology reported in Yesso scallops infected with F. halioticida [START_REF] Meyer | Disease and mortality among Yesso scallops Patinopecten yessoensis putatively caused by infection with Francisella halioticida[END_REF][START_REF] Kawahara | Francisella halioticida, Identified as the Most Probable Cause of Adductor Muscle Lesions in Yesso scallops Patinopecten yessoensis Cultured in Southern Hokkaido, Japan[END_REF]) is what prompted this investigation to test for the presence of F. halioticida in French mussel samples from previous studies.

To the best of our knowledge, F. halioticida has never been reported in mussels. The current study reveals the presence of F. halioticida, identified by polymerase chain reaction (PCR) and by in situ hybridisation (ISH), in cultured French mussels experiencing mortalities.

MATERIALS & METHODS

The 341 mussels analysed for detection of Francisella halioticida in the current study came from two overlapping multi-parametric studies in France: the first was conducted in 2017 with mussels taken from three sites in northern Brittany [START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF], two sites (Brest and Lannion) were affected by high mussel mortality while the third site (St Brieuc) did not experience mortalities. In this first study, mussel seed batches were periodically checked for mortality and sampled (live mussels randomly taken) for diagnosis. The second study was conducted in 2017 and 2018 by analysing mussels for obligatory reporting of increasing mortality (Council Directive 2006/88/EC) in Normandy and Hauts-de-France region [START_REF] Charles | Étude des organismes pathogènes, des conditions physiologiques et pathologiques impliqués dans les mortalités anormales de moules (Mytilus sp.)[END_REF]; in this second study, the sampling agents (from Administration and/or the Regional Shellfish Farming Committee) collected mostly moribund mussels and the information on mortality was obtained from Administration reports [START_REF] Canier | REPAMO 2017 -Bulletins de la surveillance -Janvier à Décembre[END_REF][START_REF] Combette | REPAMO 2018 -Bulletins de la surveillance -Janvier à Décembre[END_REF]).

In the two previous studies, mussels were analysed using histological, bacteriological and molecular diagnostic procedures; and although a number of pathogens and pathological conditions were detected none of them were believed to be the primary cause of high mussel mortality in northern France [START_REF] Charles | Étude des organismes pathogènes, des conditions physiologiques et pathologiques impliqués dans les mortalités anormales de moules (Mytilus sp.)[END_REF][START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF].

The 17 mussel groups analysed in the current study are listed in Table 1 together with some data derived from the two previous studies, namely the sampling place, the corresponding farming method, the cumulative mortality at sampling and the condition of mussels in the samples (either alive or moribund). It is noteworthy that the mussels collected during the first previous study, designed in advance to monitor mortality, were alive (valves remained closed when out of seawater), while the mussels collected during the second previous study, based on sampling immediately after knowing mortality events, were primarily moribund specimens (valves were gaping and non-responsive) or a mixture of live and moribund. Two molecular diagnostic procedures were used targeting F. halioticida DNA, polymerase chain reaction (PCR) and in situ hybridisation (ISH). Information on how the mussels had been processed in the previous two studies was reported by [START_REF] Charles | Étude des organismes pathogènes, des conditions physiologiques et pathologiques impliqués dans les mortalités anormales de moules (Mytilus sp.)[END_REF] and [START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF]. Briefly, for histological analysis, an approximately 5 mm thick transverse section of mussel tissue containing mantle lobes, visceral mass (gut, digestive gland) and gills was excised, fixed in Davidson's solution for 48h, kept in 70% ethanol and then dehydrated through an ascending ethanol series and embedded in paraffin wax. From each mussel a pool of minced tissues from foot, adductor muscle, gills, mantle, and digestive gland were preserved in a microcentrifuge tube at -20 °C.

A total of 332 mussels were used for PCR assays, 20 mussels per group except for groups 11 and 15 in which 16 mussels were used. DNA was extracted from preserved tissues of individual sample mussels, using a QIAamp DNA minikit ® (Qiagen, Courtaboeuf, France) following the manufacturer's protocol for blood or body fluids (except elution was performed using 60 μL Qiagen elution buffer AE). Following extraction, the DNA extracts of 5 mussels were pooled, except in groups 11 and 15 in which DNA pools derived from 4 mussels. Thus, 4 DNA pools were produced for each mussel group. PCR assays were conducted using the primer pair Fh-rpoB/F and Fh-rpoB/R (Brevik et al. 2011) designed to amplify the gene corresponding to bacterial DNA-directed RNA polymerase beta subunit (rpoB). For the conventional PCR, typical 25 μL reaction mixtures contained 12.5 μL of Premix Ex Taq ® 2 X Takara ® (Lonza, Verviers, Belgium), 0.5 μL of each primer (20 μM), 9.5 μL of purified water and 2 μL of DNA template (replaced with 2 μL of purified water in the negative control). PCR amplifications were carried out in a T100 TM Thermal Cycler (Bio-Rad, France) using the following thermal cycling conditions: 10 s at 95 °C; 30 cycles for 10 s at 95 °C, 30 s at 55 °C, 40 s at 72 °C and a final extension of 3 min at 72 °C. All PCR products were purified and sequenced at LABÉO Frank Duncombe (Caen Cedex, France), and then analysed using the National Center for Biotechnology Information (NCBI) Basic Local Alignment Search Tool (BLAST).

For ISH assays, only mussels from the first study could be retrospectively analysed; nine archived paraffin blocks of mussels showing different histological characteristics were selected: three healthy mussels (from groups 5, 6 and 7), two mussels with heavy haemocyte infiltration in the visceral mass (from groups 1 and 5), two mussels with granulomas (from group 5), one mussel with massive haemocyte diapedesis through stomach epithelium (from group 4) and one mussel with abundant degenerated digestive gland tubules filled with necrotic cells (from group 6).

Samples analysed by ISH corresponded to mussels different from those analysed by PCR. The ISH assays were performed as described by [START_REF] Meyer | Disease and mortality among Yesso scallops Patinopecten yessoensis putatively caused by infection with Francisella halioticida[END_REF] except that a 1 % Light Green counterstain was substituted for Bismark Brown to yield better

RESULTS

The PCR results as well as the ISH results for detection of Francisella halioticida DNA in French Mytilus spp., are summarized in Table 1. BLAST results of all the PCR amplified fragments were a 99.9 % match to the reference rpoB sequences for F. halioticida (Brevik et al. 2011) (GenBank accession numbers JF290381 and JF290374). Regarding the first study conducted in northern Brittany during 2017 (groups 1 to 8), F. halioticida DNA was detected by PCR in mussels collected in May at Lannion and Brest which had mortality rates of 20 % and 56.5 % respectively, after 3 months of culture, but it was not detected in samples from St Brieuc, an area that was not affected by abnormal mortality (only 4.5 %). Concerning the groups collected following reports of abnormal mortality in Normandy and Hautsde-France region in 2017 and 2018 (groups 9 to 17), PCR-positive results were obtained from every group except the one coming from Chausey's Island.

Although ISH assays were not conducted on all groups, positive reactions were observed within granuloma lesions from two specimens (Figure 1). These lesions were characterized by an intense accumulation of granulocytes surrounded by layers of flattened, epithelioid cells (Figure 1A & Figure 2A); and every tested mussel with granulomas showed positive reaction. Other types of pathology, such as heavy haemocyte infiltration and diapedesis, were negative via ISH; furthermore, mussels with no significant histopathology were also negative via ISH. Thus, infection with F. halioticida appears to be associated with the host's response and formation of granulomas (Figure 2), which is common in vertebrates infected with other Francisella species [START_REF] Colquhoun | Francisella infections in farmed and wild aquatic organisms[END_REF].

DISCUSSION

The occurrence of Francisella halioticida in French areas affected by high mussel mortality was revealed, whereas it was not detected in two areas, St Brieuc, which was not affected by abnormal mussel mortality, and Chausey's Island, where mussels are known to suffer from intense predation by ducks, gulls and spider crabs, which were most likely contributing to the mortalities (CRC Normandie-Mer du Nord 2017, pers. com. 2019).. The findings from this study support a hypothesis proposed by [START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF] concerning the possibility of missing an uncultivable bacterial pathogen with the classical medium used for bacterial culture. This could also be the case of the network for the monitoring of mollusc pathologies (REPAMO), coordinated by the French institute Ifremer, that monitors reports of mussel mortalities and conducts bacteriology using traditional culture media. Isolation and culture of Francisella halioticida requires a special culture medium (Modified Eugon Agar) and it can take several weeks for the bacteria to grow [START_REF] Kamaishi | Mass mortality of giant abalone Haliotis gigantea caused by a Francisella sp. bacterium[END_REF][START_REF] Kawahara | Francisella halioticida, Identified as the Most Probable Cause of Adductor Muscle Lesions in Yesso scallops Patinopecten yessoensis Cultured in Southern Hokkaido, Japan[END_REF]. The limited sample size of this study prevents us from giving conclusive results; nevertheless, given that F. halioticida is known to be highly pathogenic in abalone and Yesso scallops, its detection in mussels from Normandy and northern Brittany areas affected by high mussel mortality highlights the need for further research concerning the association between mussel mortalities and infection with F.

halioticida. Regarding the lack of geographical connection of the French mussel culture areas with the previous locations where F. halioticida was found, Canada and Japan, the importation of M. yessoensis in the beginning of the year 1987 in a hatchery of Brittany for several aquaculture trials [START_REF] Buestel | Programme de recherche pour la mise au point d'une méthode d'élevage des coquilles Saint-Jacques en Méditerranée -Rapport d'avancement des travaux[END_REF]) could be considered as a potential hypothetical connection. It would be interesting to retrospectively test for the presence of F. halioticida DNA in samples of past mortalities, given that since 2015, several studies highlighted the presence of numerous inflammatory granulomas in French mussels affected by mortalities in the Atlantic coasts without any cause being identified [START_REF] Robert | Réseau national d'observation de la moule bleue Mytilus edulis, Mytilobs, campagne[END_REF], Travers et al. 2016[START_REF] Pépin | Mortalités de moules bleues dans les secteurs mytilicoles charentais et vendéens: description et facteurs liés -MORBLEU[END_REF][START_REF] Combette | REPAMO 2018 -Bulletins de la surveillance -Janvier à Décembre[END_REF][START_REF] Charles | High mortality of mussels in northern Brittany -Evaluation of the involvement of pathogens, pathological conditions and pollutants[END_REF]. In fact, some of these authors have even suggested a link

Figure 1 .

 1 Figure 1. In situ hybridisation (ISH) assay showing positive reaction (dark blue deposits). (A) Section through the digestive gland of a mussel showing a granuloma (*) with tiny dark blue deposits within them. (B) Higher magnification of the inside of a granuloma showing dark blue deposits (arrows) indicative of positive ISH reaction within cells.

Figure 2 .

 2 Figure 2. Micrographs of a histological section consecutive to that of the Figure 1, stained with Harris' haematoxylin and eosin. (A) Large granuloma in the connective tissue of the digestive gland consisting of a mass of accumulated granulocytes (*) surrounded by layers of flattened, epithelioid cells (arrow). (B) Higher magnification of the inside of the granuloma showing numerous granulocytes with phagocytosed cells.

Table 1 :

 1 Summary of mussel (Mytilus spp.) samples collected from France and test results for Francisella halioticida detection by conventional Polymerase Chain Reaction (PCR) and in situ hybridisation (ISH) assays. Each pool in the PCR assays resulted from pooling DNA from 5 mussels except in mussel groups 11 and 15 in which DNA from 4 mussels was pooled. Individual mussels assayed by PCR were different from those assayed by ISH. Bouchot: intertidal mussel stakes; Mix: mixture of alive and moribund mussels; NA: missing data; NP: not performed.

	Group	Site	Region	Sampling date	Farming method	Cumulative mortality (%) at sampling time	Mussels condition	Number of pool of mussels analysed by PCR	PCR-results (positive pools/total pools)	ISH-results (positive mus-sels/total mussels analysed)
	1	Lannion	Brittany	24/02/2017	Longlines	5	(a) (b)	Alive	4	0/4	0/1
	2	Brest	Brittany	27/02/2017	Bouchot	12.5	(a) (b)	Alive	4	0/4	NP
	3	St Brieuc	Brittany	27/02/2017	Bouchot	2.6	(a) (b)	Alive	4	0/4	NP
	4	Brest	Brittany	24/04/2017	Bouchot	45	(a) (b)	Alive	4	0/4	0/1
	5	Lannion	Brittany	17/05/2017	Longlines	20	(a) (b)	Alive	4	1/4	2/4
	6	St Brieuc	Brittany	28/05/2017	Bouchot	4.5	(a) (b)	Alive	4	0/4	0/2
	7	Brest	Brittany	30/05/2017	Bouchot	56.5	(a) (b)	Alive	4	2/4	0/1
	8	Brest	Brittany	20/09/2017	Bouchot	66	(a) (b)	Alive	4	0/4	NP
	9	Oye-Plage	Hauts-de-France	20/03/2017	Bouchot	25	(c)	Moribund	4	3/4	NP
	10	Oye-Plage	Hauts-de-France	24/10/2017	Bouchot	10-50	(c)	Moribund	4	4/4	NP
	11	Tardinghen	Hauts-de-France	14/03/2018	Bouchot	10-50 (d)	Moribund	4	3/4	NP
	12	Donville-les-Bains	Normandy	19/06/2018	Bouchot	50	(d)	Mix	4	2/4	NP
	13	Agon-Coutainville	Normandy	16/07/2018	Bouchot	NA	(d)	Mix	4	2/4	NP
	14	Chausey	Normandy	25/09/2018	Bouchot	10-50	(d)	Mix	4	0/4	NP
	15	Wimereux	Hauts-de-France	15/10/2018	Wild	NA	(d)	Moribund	4	4/4	NP
	16	Tardinghen	Hauts-de-France	13/10/2018	Bouchot	NA	(d)	Alive	4	2/4	NP
	17	Dannes	Hauts-de-France	26/11/2018	Bouchot	25	(d)	Moribund	4	4/4	NP

(a)

[START_REF] Bernard | Bilan de l'observatoire des mortalités de moules en Bretagne Nord pour la saison 2016-2017 et premiers éléments sur les organismes pathogènes présents -Rapport technique[END_REF]

;

(b)

Charles et al. (2020) ; (c) Canier et al. (2018); (d) Combette (2018)

contrast.

between the presence of the granuloma lesions and mortality. In addition, since 2014, mussels were not the only molluscs in France to be affected by unexplained mortalities: abalone (Haliotis tuberculata), cockles (Cerastoderma edule), clams (Ruditapes philippinarum) and scallops (Pecten maximus, Chlamys varia) were also affected [START_REF] François | Bilan 2014 du réseau Repamo -Réseau national de surveillance de la santé des mollusques marins[END_REF], Lupo et al. 2016[START_REF] Amigo | Repamo 2016bulletins de la surveillance janvier à juin et décembre[END_REF].

Considering the contagious and virulent nature of the Francisella spp. bacteria, and particularly F. halioticida in scallops, the detection of this bacterial species in French mussels opens up a new path to follow for further research concerning the cause of high mussel mortalities that have been occurring in northern France for the past few years. A field approach combined with laboratory challenges assessing the Koch's classic criteria would contribute to clarify the role of F. halioticida in mussel high mortality in N France.

ACKNOWLEDGEMENTS

The authors thank the Regional Shellfish Committee of North Brittany for supporting sampling in northern Brittany and the local mussel farmers for making this possible.

The Regional Shellfish Committee of Normandy/ North Sea (CRC Normandie/Mer du Nord) and the Laboratory LABÉO financially supported this study. Maud Charles received cofounding from the Normandy region and from the laboratory LABÉO Frank Duncombe.

CONFLICT OF INTEREST

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

DATA AVAILABILITY STATEMENT

The authors confirm that one part of the data supporting the findings of this study is available within the article and the other part is available on request from the corresponding author.

LITERATURE CITED

Allain G, Bernard I ( 2016