Caterina Lupini
email: caterina.lupini@unibo.it

Mattia Cecchinato

Enrico Ricchizzi

Clive J Naylor

Elena Catelli

A Turkey rhinotracheitis outbreak caused by the environmental spread of a vaccine derived avian metapneumovirus

Keywords: Avian metapneumovirus, reversion to virulence, vaccine, TRT

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

Turkey rhinotracheitis outbreak caused by the environmental spread of a vaccine derived avian metapneumovirus

Introduction

Avian metapneumovirus (aMPV) causes an upper respiratory tract infection in turkeys and in some other avian species, including chickens. In laying or breeding birds the virus can cause drops in egg production and poor eggshell quality. These conditions are associated with serious economic losses in unprotected birds especially when secondary pathogens become involved [START_REF] Cook | Avian pneumovirus infections of turkeys and chickens[END_REF][START_REF] Gough | Avian Metapneumoviruses[END_REF]. aMPV is classified in the genus Metapneumovirus, family Paramyxoviridae (Pringle, 1998), and contains a non-segmented genome, negative stranded RNA of approximately 13,500 nucleotides, containing 8 genes (3'-N-P-M-F-M2-SH-G-L-5'), currently believed to code for nine proteins [START_REF] Easton | Animal pneumoviruses: molecular genetics and pathogenesis[END_REF]. aMPVs are grouped in four subtypes (A, B, C and D) based on nucleotide sequence differences [START_REF] Juhasz | Extensive sequence variation in the attachment (G) protein gene of avian pneumovirus: evidence for two distinct subgroups[END_REF][START_REF] Bayon-Auboyer | Nucleotide sequences of the F, L and G protein genes of two non-A/non-B avian pneumoviruses (APV) reveal a novel APV subgroup[END_REF]Seal et al., 2000). A and B subtypes have been prevalent in Europe since its first appearance [START_REF] Hafez | Presence of avian pneumovirus type A in continental Europe during the 1980s[END_REF][START_REF] Jones | Viral respiratory diseases (ILT, aMPV infections, IB): are they ever under control[END_REF].

The first aMPV virus isolation in Italy was made in 1987 [START_REF] Fabris | Rinotracheite infettiva del tacchino: osservazioni sul campo ed indagini di laboratorio[END_REF]. Since 2001, subtype B field viruses were frequently found in the most poultry populated areas in the North of Italy using sporadic and longitudinal, virological and molecular testing of turkeys and chickens (Catelli et al., 2004a;b;[START_REF] Catelli | Dati epidemiologici sulle infezioni da Pneumovirus Aviare in Italia. Atti del convegno Malattie virali respiratorie e problemi di produzione[END_REF][START_REF] Catelli | Field avian metapneumovirus evolution avoiding vaccine induced immunity[END_REF]Terregino, 2010). In general, subtype B vaccines were used for its control, although subtype A vaccine has also been used, but to a more limited degree from approximately 2002. At this time no subtype A virus had been detected in the field in Italy. Although Italian subtype B viruses have been widely studied (Catelli et al., 2010;[START_REF] Cecchinato | Avian metapneumovirus (AMPV) attachment protein involvement in probable virus evolution concurrent with mass live vaccine introduction[END_REF] less information has been available for Italian aMPV A subtypes. In the only published study, disease was observed at 3 weeks of age which was shown to be caused by a vaccine derived revertant virus (Catelli et al., We now report details of another Italian subtype A aMPV isolated during a turkey rhinotracheitis (TRT) outbreak on a large farm of seven-week-old male turkeys which had only been vaccinated with subtype B vaccine at seven days of age.

Materials and Methods

Farm. Samples were collected during an outbreak of respiratory disease from a farm of approximately 20,000 male turkeys in a lowly populated poultry area (Umbria Region), belonging to an Italian integrated poultry company. Clinical signs started when turkeys were 45 days of age and were characterized by sneezing, ocular discharge and nasal exudate. Mortality reached 10-11%. All the birds were vaccinated by eye drop administration at seven days of age with a live aMPV subtype B vaccine. Subtype A vaccine was not in use and had not been in use for at least six months in other turkey farms within a distance of approximately five km. Within that distance, all turkey farms belonged to the same poultry company and were subjected to the same vaccination scheme. Therefore subtype A vaccine was not in use in the locality. aMPV isolation. Virus isolation was performed in chicken embryo tracheal organ cultures (TOC) [START_REF] Cook | The use of chicken tracheal organ cultures for the isolation and assay of avian infectious bronchitis virus[END_REF]. The swabs pool was centrifuged at 1000g for 10 min. at +4°C. The supernatant was used to inoculate TOCs after filtration through 0.45µm membrane filters. Ciliostasis was taken as the initial indicator of aMPV detection. RNA was extracted from the TOCs (QIAamp viral mini kit, Quiagen) and subsequent RT-PCR was used to confirm virus presence and to determine the subtype.

Sampling

RT-PCR and sequencing. RNA was typed using a subtype specific RT-nested PCR based on AMPV attachment (G) gene sequence which was able to differentiate A and B subtypes [START_REF] Cavanagh | Longitudinal field studies of infectious bronchitis virus and avian pneumovirus in broilers using type-specific polymerase chain reactions[END_REF]. Wizard® SV Gel and PCR Clean-Up System (Promega) was used for amplicon purification. Sequencing was performed at BMR Genomics (Padua, Italy) with the 3730xl DNA Analyzer (Applied Biosystems, USA).

Comparison to available subtype A aAMPV G gene sequences and phylogenetic

analysis. The partial G gene sequence (268 bp) of the isolated virus was compared to available sequences. Two live subtype A aMPV vaccines were available in Italy, which had both been derived by empirical passage of the UK #8544 strain. The sequence of vaccine 1 had been previously determined [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF], while the partial G gene sequence of vaccine 2 was determined using the method described above. Sequence of 309/04, a virulent revertant of vaccine 1, was also included [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF]. The remainder of the sequences were sourced from GenBank and were: UK#8544 (accession number: DQ666911.1), CVL/14/1 (accession number: Sequence comparisons were made using the CLUSTAL W web interface [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF]. Phylogenetic analysis was carried out under distance criterion, with neighbour-joining as algorithm, using MEGA4 software. Bootstrap values were obtained with 1,000 replicates.

Full genome sequencing and characterization of the aMPV isolate. The entire genome of the isolated virus was amplified and sequenced (GenBank accession number JF424833) as previously described [START_REF] Naylor | Development of a reverse-genetics system for avian pneumovirus demonstrates that the small hydrophobic (SH) and attachment (G) genes are not essential for virus viability[END_REF]. This was compared using the CLUSTAL W web interface [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF] to the full sequence of vaccine 1 [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF], 309/04 strain [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF] and UK #8544 strain (accession number DQ666911).

Partial sequencing of vaccine 2.

Where vaccine 1 and UK#8544 sequences differed (genome positions in antigenome sense: 2941, 3553, 3825, 5055 PubMed , 5140, 5929, 6358, 10022 and 11624), the genome locations were sequenced in vaccine 2 using previously described methods [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF]. These were compared to the homologous sequences of vaccine 1 and UK#8544, using the CLUSTAL W web interface [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF] Virulence assessment of the aMPV isolate. The virulence of the isolated aMPV was determined by scoring clinical signs following eyedrop inoculation of 1-day-old poults. Thirty 1-day-old commercial turkeys were housed in three positive pressure isolators with 10 poults in each. Poults in isolator 1 were inoculated with the aMPV isolate at a dose of 3.5 log 10 median ciliostatic doses (CD 50 ) per bird. Poults in isolator 2 were inoculated with the vaccine 1 at a dose of 4.5 log 10 median tissue culture infectious doses (TCID 50 ) which represented a ten times overdose compared to the manufacturer's recommendations. Virus and vaccine 1 titres were calculated in TOCs and VERO cells respectively, using the method of Reed and Muench (Reed & Muench, 1938). The third group was inoculated with sterile water and kept as the negative controls. Clinical signs were scored daily up to 14 days post infection as previously described [START_REF] Naylor | Demonstration of a virulent subpopulation in a prototype live attenuated turkey rhinotracheitis vaccine[END_REF]. Briefly: 0. no clinical signs; 1. clear nasal exudate; 2. turbid nasal exudate; 3. swollen infra-orbital sinuses and/or frothy eyes. At the end of the trial all the birds were humanely killed.

The distribution of the clinical scores within the groups was first tested using the Kruskal-Wallis non parametric one-way ANOVA. The differences between groups were tested using the Mann-Whitney test. A p-value <0.05 was considered statistically significant.

Results

aMPV isolation. The ciliostatic agent isolated on TOCs was confirmed to be aMPV subtype A by RT nested PCR and was named IT/Ty/A/259-01/03. revealed that the strain IT/Ty/A/259-01/03 formed a separate cluster containing vaccine 1 and 309/04 (Figure 1), both derived from the UK#8544 strain. No other sequence, including that of vaccine 2, which was also derived from the UK#8544 strain, formed part of this cluster. At the detailed sequence level the vaccine 1, IT/Ty/A/259-01/03 and 309/04 strains possessed an exclusive cytosine residue at genome position 6358 which was a uracil in all the other strains with which it was compared (Table 1).

Full genome sequencing and characterization of the aMPV isolate. The complete sequencing of IT/Ty/A/259-01/03 showed that it possessed eight out of the nine nucleotide substitutions differentiating vaccine 1 and its progenitor, while the nucleotide at genome position 3553 occurring in the F gene matched the vaccine progenitor (Table 2). In addition there were four further mutations unique to IT/Ty/A/259-01/03 which occurred in the phosphoprotein (P) gene, matrix protein (M) gene, G gene and polymerase (L) gene. These are highlighted in Table 2. Only the mutations in the P, M and F genes coded for amino acid changes.

Partial sequencing of vaccine

Discussion

For the first time, a virulent, vaccine-derived aMPV has been detected in turkeys to which the vaccine had not been administered. Sequencing showed that IT/Ty/A/259-01/03 had derived from vaccine 1, and experimental infection of naïve turkeys showed it to be virulent. Furthermore the virulent nature of the virus and the further mutations seen in addition to those characteristic of the vaccine, suggest that this was not a contamination of unpassaged vaccine as perhaps might occur due to personnel and fomite movement from premises contaminated following vaccine application. Vaccine 1 had not been applied on the farm either in that period, or for at least six months previously. Hence, at the very least, either the vaccine derived virus had persisted in the immediate environment for some period prior to its detection, or had Field reversion to virulence of live vaccines is more likely to occur in avian than in other species because opportunities for horizontal transfer of vaccines are much increased due to the higher numbers of vaccinates. Convenience of mass live vaccine application increases the risk of vaccine virus circulation in birds. Also the precise nature of that application can have a bearing on persistence and disease. In a previous study, aMPV mass vaccine application by spray resulted in less than complete vaccine take, though unresponding birds did eventually receive virus, presumably due to bird to bird contact [START_REF] Ganapathy | Vaccination of commercial broiler chicks against avian metapneumovirus infection: a comparison of drinking-water, spray and oculo-oral delivery methods[END_REF]. It is generally believed, and has been anecdotally demonstrated (Wilding P., personal communication) that in this incomplete vaccination scenario, reversion to virulence can be encountered and it was noticeable that many birds in the Ganapathy study only started to show vaccine replication at three weeks of age. In another example, reversion to virulence has almost certainly occurred after live infectious laryngotracheitis vaccine use and in some cases, current field strains appear to be their direct descendants [START_REF] Oldoni | Pathogenicity and growth characteristics of selected infectious laryngotracheitis virus strains from the United States[END_REF][START_REF] Moreno | Epidemic of Infectious Laryngotracheitis in Italy: characterization of virus isolates by polymerase chain reaction -restriction fragment length polymorphism (PCR-RFLP) and sequence analysis[END_REF]. In the light of these findings, the testing of genetic stability of avian vaccines during the live vaccine registration process could be considered to be of the highest priority. Where acceptable stability has not been obtained by conventional empirical passage, rational attenuation by reverse genetics may be required.

In a previous study aMPV vaccine virus revertant 309/04 was isolated and characterized from a farm where vaccine had been applied [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF]. In the resultant discussion the view was taken that the revertant had probably arisen directly from the applied vaccine, but nonetheless it was accepted that it may have already been present in the environment. This was postulated because the five mutations, differentiating 309/04 from the vaccine seemed excessive to occur in just three weeks from vaccination. In the light of this information and the findings of the current study, it appears reasonable to conclude that both these revertant viruses deriving from vaccine 1 may have been present in the environment for some period prior to detection and isolation.

Circulation of aMPV vaccines in vaccinated or unvaccinated poultry flocks has been hypothesized, based on partial sequence analysis or epidemiological data [START_REF] Banet-Noach | Characterization of Israeli avian metapneumovirus strains in turkeys and chickens[END_REF][START_REF] Chacon | Detection by reverse transcriptase-polymerase chain reaction and molecular characterization of subtype B avian metapneumovirus isolated in Brazil[END_REF]Ongor et al., 2010). The complete sequencing of the IT/Ty/A/259-01/03 strain and comparison with vaccine and the vaccine progenitor showed the presence of eight out of the nine vaccine markers previously reported [START_REF] Catelli | Demonstration of loss of attenuation and extended field persistence of a live avian metapneumovirus vaccine[END_REF]. These mutations uniquely arose during the generation of vaccine 1 by passage of #8544, primarily in Vero cells. A simplistic calculation of permutations, where 13,373 genome bases could have one of 4 identities, shows that the probability of achieving the vaccine sequence by random mutation of a field strain would be 1 in 4 13,373 or 10 2,012. . In reality, mutation of most nucleotides would not result in viable virus hence this probability would actually be higher, though remain tiny. Furthermore, when another aMPV subtype A field strain was grown in Vero cells, none of the vaccine marker mutations were found (Naylor et al, 2004). This indicates that the marker mutations were not automatically generated simply by adaptation to Vero cell culture. By any reasonable assessment the isolated virus was without any doubt whatsoever, of vaccine origin. There can be no doubt that IT/Ty/A/259-01/03 and by similar argument 309/04, had derived from vaccine 1. In further support of this point, vaccine 2 which also derived via cell culture from the It can be deduced that the five mutations occurring during the field passage of vaccine 1 are responsible for the observed increase in virulence. None of the mutations occurred in recognized transcription controlling regions, hence only the coding changes in P, M and F genes are likely to have affected the phenotype. It is tempting to hypothesize that the mutation at genome position 3553, resulting in a predicted amino acid substitution in the F protein, is likely to have been of most significance because this had been mutated in the generation of the vaccine from UK#8544 and then back mutated during field passage. However consideration of the nature of the amino acid substituted makes this less likely because the V residue found in UK#8544 and in the revertant is of similar character to the A found in vaccine 1. Nonetheless the effects of this substitution may still remain significant. On the other hand, the A to E and S to Y predicted substitutions in the P and M proteins respectively are likely to be more influential because they represent charge and size differences which have greater potential to change protein properties. To determine which mutations increased the virulence, individual mutations would need to be made by reverse genetics as carried out for 309/04 in a previous study [START_REF] Brown | A single polymerase (L) mutation in avian metapneumovirus increased virulence and partially maintained virus viability at an elevated temperature[END_REF]. Whichever of these mutations eventually prove responsible for the virulence increase, it is interesting to make comparison with 309/04 because in that instance virulence was conclusively demonstrated to arise from a charged amino acid substitution in the polymerase gene. It therefore seems likely that reversion and attenuation might result from changes in any of a number, or possibly all, of the aMPV genes.

A further finding of this study is its confirmation of the limitations of cross protection between subtypes. Van de Zande et al. (2000) showed that the duration of heterologous protection declined while homologous protection was maintained. In this instance, the A subtype virus was found some six weeks after vaccination at a time when heterologous protection might have critically deteriorated.

The initial G gene sequencing, which showed the difference between vaccine, the progenitor and available sequences for numerous international subtype A aMPVs at genome position 6358, conveniently fell within the genomic region amplified by two commonly used diagnostic RT-PCR protocols [START_REF] Bayon-Auboyer | Comparison of F-, G-and N-based RT-PCR protocols with conventional virological procedures for the detection and typing of turkey rhinotracheitis virus[END_REF][START_REF] Cavanagh | Longitudinal field studies of infectious bronchitis virus and avian pneumovirus in broilers using type-specific polymerase chain reactions[END_REF]. Sequencing of this fragment during future routine diagnostic studies will give some indication of whether vaccine 1 or its revertants were present.

Based on the data currently available in Italy, and bearing in mind the relatively low number of subtype A viruses detected, it would currently be possible to hypothesize that most or all subtype A isolations have arisen from the spread of subtype A vaccine. Because the 6358 mutation occurs in this commonly used RT-PCR region, it would be relatively easy to sequence available amplicons to disprove or add weight to this notion. A similar single one base mutation approach has been conveniently used previously to distinguish between subtype B aMPV field viruses and vaccine and this was able to be coupled to a restriction endonuclease digestion for rapid differentiation [START_REF] Lupini | Discriminazione rapida fra ceppi vaccinali e di campo di Metapneumovirus aviare mediante tecnica RFLP[END_REF].

The study again highlights, and this time with perhaps greater urgency, the need for live aMPV vaccines that do not revert to virulence. It is generally, and probably correctly, believed that aMPV vaccine reversion to virulence is a rare event and also one which can be avoided by stringent vaccine application. The current study highlights the possibility that even where a farm follows such stringent procedures they may still be vulnerable from other farms working to a lower application standard. Ongor, H., Karahan, M., Kalin, R., Bulut, H. & Cetinkaya, B. (2010). Detection of avian metapneumovirus subtypes in turkeys using RT-PCR. The Veterinary Record, 166, 363-366. Pringle, C.R. (1998)

P C C A b A A E c 2197 M C C A b S S Y c 2941 T Transcriptionstart G A A none d none none 3553 F U C U V A V 3825 F C A A E K K 5055 M2 A G G K R R 5140 M2 U C C V N N 5929 SH-G intergenic A G G none none none 6358 G U C C F F F 7069 G U U C b Y Y Y 8766 L C C U b L L L 10022 L U G G V V V 11624 L U C C Y Y Y F o

  in that instance, the flocks had been vaccinated at one day of age with the same live subtype A vaccine.

.

  At seven weeks of age, when clinical signs were first observed, 15 birds were oro-pharingeally swabbed for virus isolation. Swabs were pooled into serum free Minimal Essential Medium (Invitrogen, Scotland) with antibiotics to suppress bacterial and Mycoplasma growth, kept refrigerated during transportation then stored at -80°C until processing.

  Comparison to available subtype A aMPV G gene sequences and phylogeneticanalysis. Phylogenetic analysis carried out on partial G gene sequences (268bp)

  2. Sequencing of UK#8544 derived vaccine 2 showed that it shared none of the mutations resulting from the preparation of vaccine 1 from the same progenitor. It was identical to UK#8544 in all nine positions.Virulence assessment of aMPV isolate. Strain IT/Ty/A/259-01/03 caused disease when poults were inoculated under experimental conditions. The daily mean clinical group is shown graphically in Figure2. Virus IT/Ty/A/259-01/03 caused severe clinical signs in the inoculated birds from day four. Birds in the other two groups showed only very minor signs. The mean cumulative score was 15.4 in the IT/Ty/A/259-01/03 infected group and 1.54 and 0.64 respectively in the vaccine and uninoculated groups. A very highly significant statistical difference in clinical scores was found between the IT/Ty/A/259-01/03 inoculated group and the vaccine group (U= 18531; p<0.0001). There was also a significant difference between the IT/Ty/A/259-01/03 inoculated group and the uninoculated control group (U=18937; p<0.0001). The difference between the vaccine and uninoculated control group was not significant (p=0.125).

  as a contaminant in an already back-mutated form from sources further away, by feed trucks, other vehicles or personnel.

  8544 progenitor, though not in Vero cell culture shared none of the nine mutations which arose during the generation of the otherwise similar vaccine 1.

Figure 2 .Figure 1 .

 21 Figure 2. Mean daily clinical scores following inoculation of 1-day-old poults with aMPV isolate IT/Ty/A/259-01/03, vaccine 1 or sterile water at 1-day-old.

Figure 2 .

 2 Figure 2. Mean daily clinical scores following inoculation of 1-day-old poults of isolate IT/Ty/A/259-01/03,

Table 1 .

 1 . Virus taxonomy--San Diego 1998. Archives of Virology, 143, 1449-1459. Reed, L.J. & Muench, H. (1938). A simple method of estimating fifty percent end points. American Journal of Hygiene, 27, 493-497. Seal, B.S., Sellers, H.S. & Meinersmann, R.J. (2000). Fusion protein predicted amino acid sequence of the first US avian pneumovirus isolate and lack of heterogeneity among other US isolates. Virus Research, 66, 139-147. Nucleotide sequences of partial G gene of subtype A aMPVs. Vaccine 1 sequence is shown in full and other residues are only specified where different from this sequence. Vaccine 1, IT/Ty/A/259-01/03 and 309/04 strains possess an exclusive cytosine residue at genome position 6358 which is a uracil in all other comparator strains.

	Terregino, C. (2010). Le principali patologia respiratorie (e non) del 2009 -I°
	quadrimestre 2010. XLIX Convegno annuale della Società Italiana di
	Patologia Aviare, Forlì, Italy. Available at:
	www.patologiaviare.org/hipgimmaginegrande.php?HIGNIdFile=972
	Van de Zande, S., Nauwynck, H., Naylor, C. & Pensaert, M. (2000). Duration of
	cross-protection between subtypes A and B avian pneumovirus in turkeys. The
	Veterinary Record, 147, 132-134.
	Seal, B.S., Sellers, H.S. & Meinersmann, R.J. (2000). Fusion protein predicted amino
	acid sequence of the first US avian pneumovirus isolate and lack of
	heterogeneity among other US isolates. Virus Research, 66, 139-147.

Table 2 .

 2 Comparison of the nucleotide and amino acid sequences of aMPVs subtype A IT/Ty/A/259-01/03, vaccine 1 subtype A and vaccine progenitor (UK#8544).

	Genome Positions a	Region of genome	UK#8544	Nucleotide Vaccine 1 IT IT/Ty/A/ 2 259-01/03	Coded amino acid UK#8544 Vaccine 1	IT/Ty/A/ 259-01/03
	1932

a Antigenome sense. b Nucleotide mutations additional to vaccine markers, unique to strain IT/Ty/A/259-01/03. c Amino acid mutations unique to strain IT/Ty/A/259-01/03. d Non coding regions.

  Figure1. Phylogenetic tree based on the alignment of nucleotide sequences (268 bp) of partial G genes of subtype A aMPVs constructed using Neighbor Joining method and bootstrap analysis (1,000 replicates) by MEGA 4 software. Strain IT/Ty/A/259-01/03 forms a separate cluster containing Vaccine 1 and the 309/04 strain. The bootstrap value is indicated at the left of the supported node.

	r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y


Antoine Lejay
email: antoine.lejay@iecn.u-nancy.fr

Géraldine Pichot
email: geraldine.pichot@inria.fr

Simulating diffusion processes in discontinuous media: a numerical scheme with constant time steps

Keywords: divergence form operators, stochastic differential equation, skew Brownian motion, Monte Carlo simulation, Euler scheme, geophysics, diffusive media with interfaces

In this article, we propose new Monte Carlo techniques for moving a diffusive particle in a discontinuous media. In this framework, we characterize the stochastic process that governs the positions of the particle. The key tool is the reduction of the process to a Skew Brownian Motion (SBM). In a zone where the coefficients are locally constant on each side of the discontinuity, the new position of the particle after a constant time step is sampled from the exact distribution of the SBM process at the considered time. To do so, we propose two different but equivalent algorithms: a two-steps simulation with a stop at the discontinuity and a one-step direct simulation of the SBM dynamic. Some benchmark tests illustrate their effectiveness.

Introduction

In many domains of physics and biology, one faces parabolic partial differential equations of type

∂C(t, x) ∂t = 1 2 ∇(D(x)∇C(t, x)) (1)
where C represents for example the concentration or the density of a population, and D is a diffusivity tensor. This equation arises naturally as it relates the variation of C(t, x) in time to the divergence of the flux (case of an incompressible fluid) in a media where the Fick law occurs (the flux J is related to the gradient of the concentration through J = 1 2 D∇C). This is for example the case of a fluid in a porous medium and ( 1) is then the Darcy law.

Equations of type (1) may be naturally solved by a Monte Carlo method by relating the value of C(t, x) with the density of particles moving randomly in the medium with a dynamic prescribed by the second-order differential operator ∇(D∇ • ). If D is smooth, then the particles move as the solutions to some Stochastic Differential Equation (SDE) X t = x + t 0 B(X s )dW s + t 0 ∇D(X s )ds, where the matrix B is such that BB T = D and W is a Brownian motion.

There are many situations where D is discontinuous, and the discontinuity may have several order of magnitude. Even if D is not smooth, a stochastic process X is associated to 1 2 ∇(D∇ • ) [START_REF] Stroock | Diffusion semigroups corresponding to uniformly elliptic divergence form operators[END_REF]. The simulation techniques we present here rely on the properties of X in a one-dimensional random media. As explained in Remark 3, they may be applied as well as the solution to the Stochastic Differential Equation (SDE) X t = x+ t 0 D(X s )dW s for a Brownian motion W , when D is discontinuous. The process is associated to the non-divergence form operator 1 2 D . Monte Carlo methods are sometimes the only available methods and have the advantage to be simple to implement in general and may be combined with other effects such as physical or chemical interactions with the medium, ... Discontinuities of D happen in many scientific domains, where the question of the simulation of X has been treated: geophysics [1-3, 5, 11, 17, 27, 31-33, 41, 43, 54-56, 59, 62], ecology [START_REF] Cantrell | Diffusion Models for Population Dynamics Incorporating Individual Behavior at Boundaries: Applications to Refuge Design[END_REF][START_REF] Ovaskainen | Biased movement at a boundary and conditional occupancy times for diffusion processes[END_REF][START_REF] Salamon | A review and numerical assessment of the random walk particle tracking method[END_REF], brain imaging [START_REF] Fieremans | Monte Carlo study of a two-compartment exchange model of diffusion[END_REF], astrophysics [START_REF] Marcowith | Postshock turbulence and diffusive shock acceleration in young supernova remnants[END_REF][START_REF] Zhang | Calculation of diffusive shock acceleration of charged particles by skew Brownian motion[END_REF], meteorology [START_REF] Thomson | Treatment of Interfaces in Random Walk Dispersion Models[END_REF], molecular dynamic [START_REF] Bossy | Probabilistic interpretation and random walk on spheres algorithms for the Poisson-Boltzmann equation in Molecular Dynamics[END_REF][START_REF] Mascagni | Monte Carlo methods for calculating some physical properties of large molecules[END_REF], ... As backed by numerical and experimental results [START_REF] Berkowitz | Laboratory experiments on dispersive transport across interfaces: The role of flow direction[END_REF], specific schemes have to be designed to take the discontinuities into account.

Although the multi-dimensional case remains open, the one-dimensional case is tractable and already arouses a lot of interest. Several approaches respect to ours, this method may be applied to the case of a variable diffusion coefficient with a drift. But it is more complex to set up.

The novelty of this article is to explain how to generate efficiently a random variate with density q(t, x, • ). Note that our algorithms are exact and do not introduce discretization errors. Generating such a random variate allows one to draw the position X (k+1)δt of the particle when its position X kδt at time kδt is known. This way, an approximation of the path of the particle may be constructed and we call this an Euler scheme, although it is different in essence from the Euler scheme for solutions to SDE.

Te paper is organized as follows. Sect. 2 presents the Euler scheme and introduces the SBM process. A recall of the dynamic of the SBM is done in Sect. 3. Then we present the main contributions of this paper which are two different but equivalent algorithms to simulate the SBM dynamic. The first one, presented in Sect. 4, relies on the computation of the first passage at the interface which is modelled by a discontinuity. The reason is that once at the interface, the simulation of the position of the particle at a given time is easily deduced from the one of a normal random variable. The second algorithm, presented in Sect. 5, relies on a direct simulation of the random variable by using the special form of the density. Sect. 6 illustrates the effectiveness of the proposed algorithms on some benchmark tests.

Notations and definitions A random variate is a realization of a random variable. In this article, every time an algorithm relies on the generation of several random variates, they are assumed to be independent. In practice, we use only a pseudo-random number generator.

We denote by N (µ, t) the Gaussian distribution with mean µ and variance t. The Gaussian density is

g(t, x) = 1 √ 2πt exp -x 2 2t ,
and we set g(x) = g(1, x).

The cumulative distribution function of the normal distribution N (0, 1) is

Φ(x) = P[G ≤ x] = x -∞ 1 √ 2π exp -y 2 2 dy = x -∞ g(y)dy,
and with set

Φ(x) = 1 -Φ(x) = P[G > x]
with G ∼ N (0, 1).

The uniform distribution on an interval [a, b) is denoted by U([a, b)). The Bernoulli distribution of parameter p is denoted by Ber(p).

A Brownian motion B with B 0 = x is simulated very easily by iteratively computing B t+δt = B t + ξ where ξ is a realization of N (0, δt) and is independent from B t .

The Skew Brownian motion and the Euler scheme

We consider the divergence form operator

L = 1 2 ∇(D∇ • )
in dimension one (note the one-half), where

D(x) = D + if x ≥ 0 and D(x) = D -if x < 0.
Following standard results in PDE (See [START_REF] Ladyženskaja | Solvability of diffraction problems in the classical sense[END_REF] for example), for a square integrable function f (resp. g), solutions of the elliptic (resp. parabolic) PDE

     Lu = f on (-1 , 2 ), u(-1 ) = α, u( 2 ) = β (resp. ∂ t v(t, x) = Lv(t, x), t > 0, x ∈ R, v(0, x) = g(x) ) is smooth on (-1 , 2 ) \ {0} (resp. R + × (R \ {0}))
, and the conditions

u(0+) = u(0-) and D + ∇u(0+) = D -∇u(0-) (2)
(resp. u(t, 0+) = u(t, 0-) and D + ∇u(t, 0+) = D -∇u(t, 0-), ∀t > 0).

The condition on the gradients in (2) is called the flux condition.

Remark 1. More generally, we may consider a problem with diffusivity D + on the right side of 0, D -on the left side of 0 and a solution which is continuous at 0 with a flux condition (1 + q)∇u(0+) = (1 -q)∇u(0-), q ∈ (-1, 1). For this, select four values (D + , D -, D + , D -) such that

D -= 1, q = D + -D - D + + D -, D + = D + D + and D -= D - D -
and use the operator L = 1 2D ∇(D∇ • ). See [START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF] for details. Following [START_REF] Marseguerra | Normal and anomalous transport across an interface: Monte Carlo and analytical approach[END_REF], D may be interpreted as a diffusivity factor (in meters) while 1/D may be interpreted as a velocity called the parametric velocity (in m/s). Thus D = D/2D is the diffusivity coefficient (in m 2 /s), while DD (in s) has the dimension of a time and is used as the mean time of the time jumps in the Continuous Time Random Walk technique.

For the sake of simplicity, we consider only divergence form operators in this article but everything could be carried to general framework as shown in Remarks 3 and 4 below.

Our aim is to design a numerical scheme for the stochastic process X associated to L with a constant time step. This means that for a given time step δt and a given position X t of the particle, how to simulate the next position X t+δt ?

Let us recall that the process X is in some sense "memoryless", which is rigorously expressed by the Markov property. This means that the future positions of (X s ) s≥t when X t is known depend only on X t and not on (X r ) r<t . Besides, these processes are homogeneous in time, which means that the distribution of (X t+θ ) θ≥0 given X t = x is the same as the distribution of (X θ ) θ≥0 given X 0 = x. These two properties are pretty convenient for writing numerical schemes with constant time steps: one has only to prescribe, when X 0 = x, how to find exactly or approximately the position y of the particle at time X δt , and then to start again by replacing x by y.

We will abusively called a constant time step scheme an Euler scheme, by similarity with the most common and widely used Euler scheme for solving SDE, although our process X belongs to a more general class of stochastic processes.

For convenience, we have placed the discontinuity at x d = 0. In the case of diffusivity coefficient D discontinuous at some point x d = 0, it is sufficient to shift the old and new positions respectively by -x d and x d for reducing this case to x d = 0. Remark 2. According to the respective values of D + and D -, it may be convenient to use different time steps when the particle is on the right or on the left of the discontinuity. As we deal here with the dynamic of the particle close to the discontinuity, we do not address this problem here.

Reduction to the Skew Brownian motion

It has been proved that the process X is solution to a stochastic differential equation (SDE) with local time (See [START_REF] Étoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF] for example):

X t = x + t 0 D(X s )dW s + γL 0 t (X), γ = D + -D - D + + D -.
Here, L 0 t (X) is the symmetric local time of X. This is a particular process which describes the amount of time the process spends to 0. It may be formally defined as

L 0 t (X) = lim →0 1 2 t 0 1 {Xs∈[-, ]} d X s ,
where X is the quadratic variation of X, which is X t = t 0 D(X s )ds. This is a continuous, non-decreasing process. However, it increases only on the closure of the set {t ≥ 0|X t = 0}, which is a set of zero Lebesgue's measure. [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]Theorem 3.7.1] and the properties of the local time yields that Y t = Φ(X t ) is solution to

With Ψ(x) = x/ √ D + if x ≥ 0 and Ψ(x) = x/ √ D -if x < 0, the Itô- Tanaka formula
Y t = Ψ(x) + W t + θL 0 t (Y ), θ = √ D + - √ D - √ D + + √ D - . (3)
The process Y is called the skew Brownian motion (SBM) of parameter θ.

The survey article [START_REF]On the constructions of the Skew Brownian motion[END_REF] contains an account about the possible constructions and the main properties of this process which was introduced in the 70's.

The change of variable Ψ has been used in [START_REF]Simulating a diffusion on a graph. Application to reservoir engineering[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF] and in [START_REF] Ramirez | Skew Brownian Motion and Branching Processes Applied to Diffusion-Advection in Heterogenous Media and Fluid Flow[END_REF][START_REF] Ramirez | A generalized Taylor-Aris formula and skew diffusion[END_REF][START_REF] Ramirez | A Note on the Theoretical Foundations of Particle Tracking Methods in Heterogeneous Porous Media[END_REF].

The density transition function q(t, x, y) of the process X is linked to the density transition function p θ of the SBM of parameter θ given by (3) by

q(t, x, y) = 1 D(y) p θ t, x D(x) , y

D(y) .

A simple, analytic expression of p θ is given by ( 5) below. The meaning of the function q(t, x, y) is that when being at x at time 0, the probability that the particle moving according to L with a discontinuity at 0 belongs to the infinitesimal interval of size dy around y is q(t, x, y)dy.

This change of variable is pretty simple and may be applied only locally. Algorithm 1 below means that when the particle X is close enough to the boundary (See Sect. 2.2), we first transform it into a Skew Brownian motion by Y t = X t / D(X t ). Second, we apply one of the schemes given below to simulate Y t+δt . Third, we find the position of the particle at time t + δt by

X t+δt = D(Y t+δt )Y t+δt .
Remark 3. Similar computations may be applied to the solution of the SDE X t = x + t 0 D(X s )dB s whose infinitesimal generator is 1 2 D . More generally, as noted in Remark 1, we may consider the process X generated by

L = 1 2D ∇(D∇ • ). In this case, define Φ(x) = β + x for x ≥ 0 and Φ(x) = β -x with β ± = D ± / D ± . Hence Φ(X) is a SBM with θ = D + D + -D -D - D + D + + D -D -
as seen in [START_REF] Étoré | On random walk simulation of one-dimensional diffusion processes with discontinuous coefficients[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF]. In particular, the SBM of parameter θ may be associated to D

+ = D + = (1 + θ)/(1 -θ) and D -= D -= 1 [38].
Remark 4. The density transition function q(t, x, y) of the process X generator by L = 1 2D ∇(D∇ • ) satisfies both the Kolmogorov backward equation

∂q(t, x, y) ∂t + 1 2D(x) ∇ x (D(x)∇ x q(t, x, y)) = 0
as a function of x and the Kolmogorov forward equation

∂q(t, x, y) ∂t = 1 2 ∇ y D(y)∇ y 1 D(y) q(t, x, y)
as a function of y. This expresses that at 0 q(t, 0+, y) = q(t, 0-, y) and

D + ∇ x q(t, 0+, y) = D + ∇ x q(t, 0-, y)
for any y, while

1 D + q(t, x, 0+) = 1 D -q(t, x, 0-) and D
+ D + ∇ y q(t, x, 0+) = D - D -∇ y q(t, x, 0-) (4
) for any x. As noted in the end of Remark 3, the SBM of parameter θ belongs to this class of processes. Thus, the first condition in ( 4) becomes (1θ)p θ (t, x, 0+) = (1 + θ)p θ (t, x, 0-) and explains the discontinuities observed in Figure 2.

The interface layer

The SBM behaves like a Brownian motion except when it reaches 0. Let τ be the first time at which the Brownian motion starting from x < 0 reaches 0. Then P x [τ < t] = 2Φ(x/ √ t) when x < 0. Fix a probability α ∈ (0, 1) and find c α > 0 such that 2Φ(-c α ) = α. Then with probability α, the SBM does not reach the interface at position 0 for x < -c α √ t. For α = 5%, c α = 1.95. For α = 1%, c α = 2.57. For α = 0.5%, c α = 2.80.

As summarized in Figure 1, the idea is then to use the usual Euler scheme with step δt for the Brownian motion when the particle does not belong to

[- √ D -c α √ δt, √ D + c α √ δt]
, and a biased Euler scheme with step δt when the particle belongs to this interval [-

√ D -c α √ δt, √ D + c α √ δt],
which we call the interface layer. By a biased Euler scheme, we mean a scheme that replicates the dynamic of the SBM.

The Euler scheme

A scheme for X will be deduced from a scheme on Y by simply applying the function Ψ -1 to the next position of the particle. Away from the interface, the usual Euler scheme (or any other efficient scheme) may be used whereas

c α √ D -δt c α √ D + δt Interface layer ≡ SBM step normal step of s. dev. √ D -δt normal step of s. dev. √ D + δt Figure 1: The interface layer [- √ D -c α √ δt, √ D + c α √ δt];
Outside the interface layer: normal step; Inside the interface layer: SBM step (notice the non Gaussian character of the corresponding density).

inside the interface, a biased scheme is proposed: See Algorithm 1. Also, there is no problem in dealing with multiple interfaces.

We now concentrate on the simulation of the biased Euler scheme that gives the next position of a particle starting within the interface layer.

The dynamic of the Skew Brownian motion

We have to focus on an approximation of the Skew Brownian motion at a given time t. This process has the Markov properties and is homogeneous in time, which means that the distribution of (Y s+t ) t≥0 given (Y r ) 0≤r≤s has the same distribution as (Y t ) t≥0 given Y 0 .

Let P x be the distribution of Y with Y 0 = x. With an abuse of notation, we use the same convention for the Brownian motion W .

The density transition function of the SBM for a discontinuity at x d = 0 has a simple expression:

p θ (t, x, y) = g(t, y -x) + sgn(y)θg(t, |y| + |x|) (5)
with sgn(x) = 1 if x ≥ 0 and sgn(x) = -1 if x < 0. This expression may be derived from the reflection principle [START_REF]On the constructions of the Skew Brownian motion[END_REF][START_REF] Walsh | A diffusion with discontinuous local time, Temps locaux, Astérisques[END_REF]. Figure 2 shows typical plots of the density. More generally, for a discontinuity at x d , the density transition function writes p θ,x d (t, x, y) = p θ (t, x -x d , y -x d ).

From (5), we deduce immediately that the SBM enjoys the same scaling property as the Brownian motion.

Data: The position X t at time t of the particle. Result: The position X t+δt at time t + δt of the particle.

if X t < - √ D -c α √ δt then /*
H = Y δt with Y 0 = x, where Y is the SBM of parameter ( √ D + - √ D -)/( √ D + + √ D -); return X t+δt = Ψ -1 (H); end
Algorithm 1: A step of the Euler scheme. Property 1 (Scaling property). For any c > 0 and any t > 0,

Y t = √ cY t/c with Y 0 = √ cx is equal in distribution to Y t with Y 0 = x.
The "symmetry" property is also useful and follows from the equality p θ (t, x, y) = p -θ (t, -x, y). Note that θ is transformed into -θ.

Besides P 0 [Y t > 0] = (1 + θ)/2.
We recover the experimental result of [START_REF] Semra | Modélisation tridimensionnel le du transport d'un traceur en milieux poreux saturé: évaluation des théories stochastiques[END_REF]: The average number of particles at time t that are on the right of the discontinuity is equal to

√ D + /( √ D + + √ D -) time the number of particles.
The next two sections present the two algorithms we propose for the simulation of the dynamic of the SBM. The first one is a two-steps algorithm that contains a stop at the interface of discontinuity, the second one yields a direct simulation of the SBM process according to its known density transition function ( 5). An alternative but more complex algorithm, that relies on some fine properties of Brownian paths, may be found in [START_REF] Lejay | Simulation of a stochastic process in a discontinuous layered medium[END_REF].

A two-steps algorithm for the simulation of the SBM evolution

This algorithm requires to simulate the first time τ the process reaches zero, if it happens before δt. This aspect is developed in B. In particular, we present some existing methods that allow one to simulate the exact time it takes for the particle to reach the interface [START_REF] Baldi | Exact asymptotics for the probability of exit from a domain and applications to simulation[END_REF][START_REF] Buchman | Simulation of stopped diffusion[END_REF][START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF][START_REF] Lerche | Approximate exit probabilities for a Brownian bridge on a short time interval, and applications[END_REF]. We also discuss the choice of a linear interpolation.

Starting at a position say x < 0, once the time τ = inf{t > 0|W t = 0} is known for a Brownian motion W with W 0 = x, one faces two options: either τ is greater than δt, which means the particle did not cross the interface of discontinuity during the time δt, either it is lower, then the particle is moved at the interface and one has to simulate its next position, during the remaining time δt -τ .

In this last case, the starting point x is then 0, and Y t = |G|, where G ∼ N (0, t), that is a normal of mean 0 and variance t, and ∈ {-1, 1} is independent from G and P[ = 1] = (1 + θ)/2. This follows immediately from the symmetry property of the Gaussian distribution and the expression [START_REF] Bechtold | Efficient random walk particle tracking algorithm for advective-dispersive transport in media with discontinuous dispersion coefficients and water contents[END_REF].

The SBM evolution step of Algorithm 1 is given by Algorithm 2.

Data: An initial position x < 0 and a time δt > 0.

Result: A realization of Y δt where Y is the SBM of parameter θ and Y 0 = x. Generate a realization (s, y) of (δt ∧ τ, W δt∧τ ) where τ = inf{t > 0|W t = 0} for a Brownian motion W with

W 0 = x; if s = δt then /*

A direct algorithm for the simulation of the SBM evolution

We give an algorithm that does not rely on computing the first passage time to the interface, and which is then an alternative to Algorithm 2. Now, let us consider the case where the particle is at a given point x < 0. The case x > 0 may be considered by using the symmetry property (Property 2). In this case, one should change x into -x, y into -y and θ into -θ.

It follows from the scaling property (Property 1) that

Y t given {Y t > 0} with Y 0 = x is distributed as √ tY 1 given {Y 1 > 0} with Y 0 = x/ √ t. We set z = x/ √ t.
Since {Y 1 = 0} is a null event under P z , it does not matter to consider the event {Y 1 > 0} instead of {Y 1 ≥ 0}. With g(y) = g(1, y) the density of the normal random variable N (0, 1), the probability that

Y 1 > 0 is κ(z) = P z [Y 1 > 0] = +∞ 0 p θ (1, z, y)dy = (1 + θ) +∞ -z g(y)dy = (1 + θ)Φ(-z).
With the scaling property, P x [Y t > 0] = κ(z) = κ(x/ √ t). If x > 0, then a similar computation or a symmetry argument (Property 2) shows that

P x [Y t > 0] = 1 -(1 -θ)Φ(z). Let us introduce the random variable Z + (z) distributed as Y 1 given Y 1 > 0 and Y 0 = z, as well as Z -(z) distributed as Y 1 given Y 1 < 0 and Y 0 = z.
Below, we study the main properties of these random variables (simulation, moments, approximations). The simulation of Y δt is then given by Algorithm 3, which is an alternative to Algorithm 2. The simulation of Z + (z) and Z -(z) are given in Algorithms 4, 5 and 6 below.

Data: An initial position x < 0 and a time δt > 0. Result: A realization of Y δt where Y is the SBM of parameter θ and

Y 0 = x. Generate a random variate U ∼ U([0, 1)); if U ≤ (1 + θ)Φ(-x/ √ δt) then Generate a random variate H = √ δtZ + (x/ √ δt); else Generate a random variate H = √ δtZ -(x/ √ δt
); end return H; Algorithm 3: The direct algorithm simulating the dynamic of the SBM.

Simulation of Z +

Assume without loss of generalities that x < 0. With (5), p θ (t, x, y) = (1 + θ)g(y -x) when y > 0 and κ(x) = (1 + θ)Φ(-x).

For a bounded, measurable function f with a support included in R + ,

E z [f (Y 1 )|Y 1 ≥ 0] = E z [f (Y 1 ); Y 1 ≥ 0] P z [Y 1 ≥ 0] = 1 Φ(-z) +∞ 0 f (y)g(y -z)dy = 1 Φ(-z) +∞ -z f (y + z)g(y)dy = E[f (z + G)|G > -z] with G ∼ N (0, 1). This means that Z + (z) is distributed like G + z given {G ≥ -z}. The distribution of G ∼ N (0, 1) given {G ≥ m} is called the Gaussian tail distribution GT (m).
In A, we present some ways to simulate a random variate with the Gaussian tail distribution GT (z). Hence, the simulation of Z + (z) is pretty easy, as shown in Algorithm 4.

Data: A initial position z < 0. Result: A random variate of Z + (z). Generate a random variate G ∼ GT (-z); return z + G;

Algorithm 4: Simulation of Z + (z) when z < 0.

Simulation of Z -

The probability that

Y 1 < 0 with Y 0 = z < 0 is 1 -κ(z) = 1 -(1 + θ)Φ(-z).
The distribution function of Z -(z) is

P -(y) = P z [Z -(z) ≤ y] = P z [Y 1 ≤ y|Y 1 ≤ 0] = 1 1 -(1 + θ)Φ(-z) (Φ(y -z) -θΦ(y + z)) , y ≤ 0,
while its density is

p - θ (y) = 1 1 -κ(z) (g(y -z) -θg(y + z)) , y ≤ 0. Yet g(y -z) is the density of z + G with G ∼ N (0, 1).
It is then possible to use a Newton algorithm in order to compute (P -) -1 (U ) with U ∼ U([0, 1)), which gives the desired random variable. However, this method may be time consuming.

An alternative way of proceeding is the following and depends on the sign of θ.

If θ > 0, then one may use a rejection method. Note that

p - θ (y) ≤ 1 1 -(1 + θ)Φ(-z) g(y -z). Indeed, g(y -z) is the density of G + z with G ∼ N (0, 1) and P[G ≤ -z] = Φ(-z) = 1-Φ(-z). Hence h(y) = g(y-z)/Φ(-z) is the density of H = z-G with G ∼ GT (z). Thus, p - θ (y) ≤ β(z)h(y) with β(z) = 1 -Φ(-z) 1 -(1 + θ)Φ(-z) . ( 6
)
A variant of the rejection techniques is the following, when one has a bound of type (6) on densities (See [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF]II.3,p. 40]): Simulate a random variate H with density h as well as a uniform random variate U until p - θ (H)/h(H) ≤ β(z)U . Then H is a realization of a random variable with density p - θ . Note that in our case,

p - θ (y) h(y)β(z) = 1 -θ g(y + z) g(y -z) = 1 -θ exp(-2yz).
Since 1 -U and U have the same distribution when U ∼ U([0, 1)), the stop condition in the rejection method takes a simple form. Algorithm 5 describes the simulation of Z -(z) when θ > 0, z < 0.

Data: A point z < 0 and a parameter θ > 0.

Result: A random variate of Z -(z). repeat Generate a random variate H = z -G with G ∼ GT (z); Generate a random variate U ∼ U([0, 1)); until U ≥ θ exp(-2Hz); Algorithm 5: Simulation of Z -(z) when θ > 0, z < 0.
If θ < 0, then one may use a composition method. Let us note that

P -(y) = Φ(-z) 1 -κ(z) Φ(y -z) Φ(-z) -θ Φ(z) 1 -κ(z) Φ(y + z) Φ(z) .
The function Φ(y -z)/Φ(-z) is the distribution function of G + z given {G ≤ -z} while Φ(y + z)/Φ(z) is the distribution function of G -z given {G ≤ z}.

Using the symmetry of the Gaussian distribution, Φ(z) = Φ(-z) for any z ∈ R and 1 -κ(z) = Φ(-z) -θΦ(-z) since Φ(-z) + Φ(-z) = 1 for any z ∈ R.

For

α(z) = Φ(-z) 1 -κ(z) = 1 -Φ(-z) 1 -(1 + θ)Φ(-z) , (7)
P -(y) may be written

P -(y) = α(z) Φ(y -z) Φ(-z) + (1 -α(z)) Φ(y + z) Φ(z) .
One may then simulate Z -(z) using the composition method [18, II.4, p. 66] given by Algorithm 6.

As previously, we are reduced to simulate Gaussian tail distributions introduced in Section 5.1 and A.

Data: A point z < 0 and a parameter θ < 0. Result: A random variate of Z -(z). Generate a uniform random variate U ∼ U([0, 1)); Compute α(z) given by ( 7

); if U < α(z) then Generate y = z -G with G ∼ GT (z); else Generate y = -z -G with G ∼ GT (-z); end return y;
Algorithm 6: Simulation of Z -(z) when θ < 0.

5.3

The moments of the random variables Z + and Z -

The moments of the random variables Z + and Z -will be used to get some uniform or Gaussian approximations of them. After standard computations, for z < 0,

E[Z + (z)] = g(z) Φ(-z) + z, E[Z -(z)] = (θ -1)g(z) 1 -(1 + θ)Φ(-z) + z, and
E[(Z + (z)) 2 ] = z g(z) Φ(-z) + z 2 + 1, E[(Z -(z)) 2 ] = -zg(z)(1 + θ) 1 -(1 + θ)Φ(-z) + 1 + z 2 .

Uniform distribution of Z + and Z -

A first approximation consists in replacing the random variables Z + (resp. Z -) by a uniform random variable

L + U (resp. -L -U ) with U ∼ U([0, 1)).
This means that the density p + θ (resp. p - θ ) is replaced by the density

(L + ) -1 1 [0,L + ) (resp. (L -) -1 1 [0,L -) ).
A natural way to choose the parameter L ± is by matching the means of Z ± and ±L ± U . Since the expectation of U is 1/2, we take

L ± = ±2E[Z ± (z)].

Half-normal approximation of Z + and Z -

The density p ± θ has a support on R ± . In order to get faster computations, one may want to replace Z + (resp. Z -) by a simpler random variable. An obvious choice consists in using a half-normal distribution with parameter σ + . By this, we mean |G| (resp. -|G|), where G ∼ N (0, σ + ) (resp. G ∼ N (0, σ -)).

The density of |G| is 2g(σ + , • ) on R + . The Kullback-Leibler divergence is a way to characterize the distance between two distributions. Hence, we choose the parameter σ + in order to reduce the Kullback-Leibler divergence between Z + and |G| with G ∼ N (0, σ + ). This divergence is defined by

K(p + θ , 2g(σ + , • )) = +∞ 0 p + θ (x) log p + θ (x) 2g(σ + , x) dx.
Optimizing over σ + , one gets that

σ + = +∞ 0 x 2 p + θ (x)dx = E[(Z + (z)) 2 ].

Benchmark tests

To show the effectiveness of our methods, we consider some test cases presented in [START_REF] Labolle | Random-Walk Simulation of Transport in Heterogeneous Porous Media: Local Mass-Conservation Problem and Implementation Methods[END_REF] with the discontinuity at x d = 0. We consider a 1D domain,

[-L, L] with L = 49. We have D = D + on [0, L] and D = D -on [-L, 0]
with reflecting boundary conditions on x min = -L and x max = L. The time step of the algorithm is δt = 0.005.

Test 1: density test

On Figure 3, we plot an approximation of the density of the particles at time T = 6 when the starting point x 0 is -5 for D + = 0.5 and D -= 10, δt = 5 × 10 -3 , and c α = 3. We use 100,000 particles and the simulation is performed using Algorithm 1 coupled with Algorithm 3 (direct SBM algorithm). The interface layer is, in that case, [-

√ D -c α √ δt, √ D + c α √ δt] = [-0.6708, 0 .15]
. At time T = 6, the probability that a particle has reached x min or x max is very small. Hence, the density may be compared to the density q(t, x, y) on the whole space and we see a perfect match. We obtain the same good result by using the two-steps Algorithm 2 instead of Algorithm 3 (not drawn).

Remark 5. In order to compare these expression with the analytic expressions found for example in [START_REF] Labolle | Random-Walk Simulation of Transport in Heterogeneous Porous Media: Local Mass-Conservation Problem and Implementation Methods[END_REF][START_REF] Marseguerra | Normal and anomalous transport across an interface: Monte Carlo and analytical approach[END_REF], one has to take into account the one-half factor so in these articles, correspond to D 2 = D + /2 and D 1 = D -/2.

Test 2: long time behavior

We also study the long time behavior of the scheme, as in [START_REF] Labolle | Random-Walk Simulation of Transport in Heterogeneous Porous Media: Local Mass-Conservation Problem and Implementation Methods[END_REF]. It is known that in the steady state regime, the particles are uniformly spread over the domain, whatever the values of the diffusivity coefficient D. If N + (t) is the number of particles on [0, L] and N -(t) is the number of particles on [-L, 0] at time t, then N -(t)/N + (t) converges to 1. Figure 4 presents the evolution of R(t) = N -(t)/N + (t) for D -= 5 and D + = 0.25 for 10,000 particles starting at the position x min = -L with Algorithm 1 coupled with Algorithm 3 or 2. The curve "No correction" corresponds to the result obtained with the use of a classical Euler scheme on the entire domain, even in the interface layer. It emphasizes the need of a specific algorithm inside the interface layer, like one of the two proposed here, to take the discontinuity into account.

On Figure 4, the horizontal red dashed lines represent the 99%-confidence interval on R(t) based on the quantiles of 4Z and the above approximation. Indeed, by the Central Limit Theorem, assuming that the steady state is reached, then P[X t ∈ [0, L]] = 1/2 and then N + (t) may be approximated by N + (t) = N (1/2 + Z(t)) where Z(t) ∼ N (0, 1/4N ), Hence,

R(t) = 1 2 + Z 1 2 -Z ≈ 1 + 4Z.
Hence, for t large enough, R(t) fluctuates around 1 but these fluctuations stay within the error range of the Monte-Carlo method.

We also check the convergence of the ratio R(t) towards 1 for different ratio D -/D + , D -/D + ∈ {2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20}. Figure 5 shows the ratio of N -/N + averaged over the time 20 000-30 000, for which the steady state is reached. Without any appropriate method to take the discontinuity into account, the results are not satisfactory as attested by the curve named "No correction". On the contrary, all simulations performed using one

Summary and future work

We have proposed two effective algorithms for the simulation of the SBM dynamic. These algorithms are exact ones and are essential tools for the Lagrangian simulation of diffusion in discontinuous media. Extensive test cases as well as comparison with other algorithms (like the ones proposed in geosciences [START_REF] Hoteit | Three-dimensional modeling of mass transfer in porous media using the mixed hybrid finite elements and the random-walk methods[END_REF][START_REF] Uffink | A random walk method for the simulation of macrodispersion in a stratified aquifer, Relation of Groundwater Quantity and Quality[END_REF]) will be provided in a forthcoming paper.

Another interesting future work will be to introduce immobile times in order to consider more complex particle fates such as absorption, stagnation and desorption [START_REF] Benson | A simple and random walk solution of multi-rate mobile/immobile mass transport equations[END_REF][START_REF] Haggerty | Multiple-Rate Mass Transfer for Modeling Diffusion and Surface Reactions in Media with Pore-Scale Heterogeneity[END_REF][START_REF] Salamon | Modeling mass transfer processes using random walk particle tracking[END_REF].

Acknowledgments

We sincerely wish to thank the Reviewers for their fruitful comments.

A Simulation of the Gaussian tail distribution

The Gaussian tail distribution GT (z) is the distribution of the normal random variable G given to {G ≥ z}. If z < 0, then a simple rejection technique may be used: one draws a normal random variate until it is greater that z.

B Simulation of the first time the particle reaches the interface

to perform such simulations. But some of them (See [START_REF] Gobet | Stopped diffusion processes: boundary corrections and overshoot[END_REF] for example) show their relevance in term of weak error, that is the error induced in functionals of type E[φ(τ, X τ )] where the first hitting time τ of 0 is replaced by its approximation. Other methods exist where the time step is reduced (See for example [START_REF] Jansons | Exponential timestepping with boundary test for stochastic differential equations[END_REF], which also use the Brownian bridge test in a way similar to the one presented in B.2 but with exponential time steps instead of constant time steps).

What is said here may be applied to more general diffusion processes. Using the scale, we work on the level of the Brownian motion which we may assume to start from x. Again, the results of Sect. 2.2 may be applied there. If x is far enough from 0, then the probability the particle crosses 0 before δt may be neglected.

B.1 Linear approximation

We perform an Euler scheme for the Brownian motion with diffusivity D + or with diffusivity D -. If the interface is crossed, the new position does not have the same sign as the starting position. This new position is kept to perform a linear approximation of τ and the new position is set to be 0 as described by Algorithm 7.

Data: An initial position x = W 0 . Result: A couple (τ, 0) or (δt, W δt ). Generate a random variate y ∼ N (0, δt); if sgn(x) = sgn(y) then return (δt, y); else Set τ ←δt|x|/(|x| + |y|); return (τ, 0); end Algorithm 7: Linear approximation of τ

The choice of τ in case of a crossing is justified by Remark 6 below in B.2.

The drawback with this algorithm is to over-estimate the first hitting time of 0: if X k and X k+1 have the same sign, it is possible that the part of the path between times t k and t k+1 crosses the boundary 0. We refine this method in the next section.

B.2 Using Brownian bridge's properties

A way to modify the naive previous scheme consists in using the Brownian bridge distribution (see Algorithm 8). Indeed, the Brownian path (W t ) t∈[0,δt] when X 0 = x and X δt = y is called a Brownian bridge. As δt goes to zero, the deviation of a Brownian bridge path from the segment linking x to y becomes smaller and smaller.

The idea is then to use its property to simulate the first time the Brownian bridge crosses 0, if it does. For this method, see [START_REF] Baldi | Exact asymptotics for the probability of exit from a domain and applications to simulation[END_REF][START_REF] Buchman | Simulation of stopped diffusion[END_REF][START_REF] Gobet | Weak approximation of killed diffusion using Euler schemes[END_REF][START_REF] Lerche | Approximate exit probabilities for a Brownian bridge on a short time interval, and applications[END_REF], ... There are two cases to consider: xy > 0 and xy < 0. Without loss of generalities, we assume that x < 0. The computations are borrowed from [START_REF] Beskos | Exact simulation of diffusion[END_REF].

The case x < 0, y < 0. Here, the process may have not crossed the boundary.

The probability that the Brownian bridge goes above a level β > y > x is

P x sup t∈[0,δt] W t ≥ β W δt = y = exp -2 δt (β -x)(β -y) .
If the Brownian bridge crosses 0, then the first time τ = δtξ/(1 + ξ) where ξ follows an inverse Gaussian distribution IG(x/y, x 2 /δt). The density r(µ, λ,

• ) of IG(µ, λ) is r(µ, λ, x) = λ 2πx 3 exp -λ(x -µ) 2 2µ 2 x .
This distribution is easily simulated using the algorithm [START_REF] Michael | Generating random variates using transformations with multiple roots[END_REF] (See also [18, p. 148]).

The case x < 0, y > 0. In this case, the process crosses the boundary. Its first hitting time τ is equal to δtξ/(1 + ξ) where ξ follows an inverse Gaussian distribution IG(-x/y, x 2 /δt).

Remark 6. The mean of the IG(µ, λ) is µ and its variance is µ 3 /λ. Thus, for ξ ∼ IG(-x/y, x 2 /δt) with x < 0 and y > 0,

ξ = |x| |y| + √ δtα
where α is a random variable with zero mean and variance |x|/|y|

5 Figure 2 :

 52 Figure 2: The density p θ (t, x, • ) of the Skew Brownian motion for θ = 0.5 at time t = 1 for three values of the initial point x.

Property 2 (

 2 "Symmetry" property). For a SBM process Y of parameter θ with Y 0 = x, the process Y = -Y is a SBM of parameter -θ with Y 0 = -x. Under P x , |Y | has the same distribution as a reflected Brownian motion |W |. This is an immediate consequence of the Itô-Tanaka formula applied to | • |.

Figure 3 :

 3 Figure 3: Normalized histogram of particles at time T = 6 obtained with Algorithm 1 coupled with Algorithm 3 (direct SBM algorithm), with 100,000 particles when x 0 = -5 for D -= 10, D + = 0.5, δt = 5 × 10 -3 , and c α = 3 in the definition of the interface layer. The black line is the density q(t, x, y).

Figure 4 :

 4 Figure 4: Convergence of the ratio R(t) = N -(t)/N + (t) vs time with the two proposed algorithms for D -/D + = 20 and 10,000 particles starting at the position x min = -L. The black curve "No correction" corresponds to the result obtained with the use of a classical Euler scheme on the entire domain, even in the interface layer.

Figure 5 :

 5 Figure 5: Convergence of the averaged ratio of N -/N + over time 20000-30000 with the two proposed algorithms for D -/D + ∈ {2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20} and 10,000 particles starting at the position x min = -L. The black curve "No correction" corresponds to the result obtained with the use of a classical Euler scheme on the entire domain, even in the interface layer.

  The two-steps algorithm simulating the dynamic of the SBM.

	The interface is not reached before t + δt	*/
	return y;
	else
	/* The interface is reached at time t + τ with τ < δt */
	Generate a Bernoulli random variate ∈ {-1, 1} with
	P[ = 1] = (1 + θ)/2;
	Generate a random variate G ∼ N (0, 1); return √ δt -s|G|;
	end
	Algorithm 2:

  3 . Hence, with α = α|y| has a variance |x|/|y|.Hence, for small values of δt, the random variable τ may be approximated by δt|x|/(|x| + |y|).Data: An initial point x and a time δt > 0. Result: A random variate (δt ∧ τ, W δt∧τ ). Generate a random variate G ∼ N (0, 1);Set y ← x + √ δtG; if sgn(x) = sgn(y) thenGenerate a random variate ξ ∼ IG(|x|/|y|, x Simulation of (δt ∧ τ, W δt∧τ ) using Brownian bridge's properties.

			return (τ, 0);
			else
			return (δt, y);
			end
		end
	Algorithm 8:
	τ δt	=	|x| |y| + 1 + |x| |y| + √	δtα √ δtα	=	|x| + |x| + |y| + √ δtα √ δtα	≈	|x| |x| + |y|	+ √	δtα	|y| |x|(|x| + |y|)	+• • • ,

2 /δt); Set τ ← δtξ/(1 + ξ); return (τ, 0); else Generate a random variate U ∼ U([0, 1)); if U < exp(-2xy/δt) then Generate a random variate ξ ∼ IG(|x|/|y|, x 2 /δt); Set τ ← δtξ/(1 + ξ);

We consider in this appendix how to simulate the process X until the (random) time τ at which it reaches the discontinuity at 0 for the first time. Until this time τ , then for a Brownian motionX t = √ D + W t (resp. X t = √ D -W t ) if X 0 > 0 (resp. X 0 < 0).The problem then reduces to the simulation of the first time the Brownian motion reaches 0.Given δt and X k , we then present several ways to get either X k+1 with X k+1 X k > 0 (the two successive positions are on the same side of the discontinuity) or the first time τ the process hits 0 when starting from X k (or equivalently, the Brownian motion W hits 0 with starting point X k / D(X k )) with τ < δt.Three methods are presented here. Other algorithms have been proposed

 [START_REF] Bechtold | Efficient random walk particle tracking algorithm for advective-dispersive transport in media with discontinuous dispersion coefficients and water contents[END_REF]This work has been supported by the ANR SIMUDMRI (ANR-10-COSI-SIMUDMRI).

B.3 Controlling exit times from intervals

If the particle is initially at position y, then let us fix x < y < z, and let τ = inf{t > 0|W t ∈ [x, z]} be the first time the Brownian motion W reaches either x or z.

The couple (τ ∧δt, W τ ∧δt ) may be simulated exactly for any δt > 0 [START_REF]Simulating a diffusion on a graph. Application to reservoir engineering[END_REF][START_REF] Lejay | A scheme for simulating one-dimensional diffusion processes with discontinuous coefficients[END_REF]. A C library is available for this purpose [START_REF] Lejay | exitbm: a library for simulating exit times and positions for the Brownian motion from simple domains[END_REF].

If the interface is at z = 0 and y < 0, then this algorithm may return (τ, X τ ) with τ < δt and X τ = x. It is then possible to iterate by replacing δt by δt -τ until the particle reaches the interface or the sum of the time is equal to δt, as seen in Algorithm 9.

Data: W 0 = y. Result: The hitting time τ of 0 if smaller than δt or the next position of the Euler scheme. Set s←0; repeat Select x < y < z with z = 0 if y < 0 and x = 0 if y > 0; Generate a random variate (τ , y ) of (τ

} and W 0 = y for a Brownian motion W ; Set s←s + τ ; Set y←y ; until y = 0 or s = δt; Algorithm 9: Simulation of (τ ∧ δt, W τ ∧δt ) using the exit time from intervals.

Algorithm 9 is more costly than Algorithm 8. Yet this algorithm is suitable in the case of multiple interfaces, since this algorithm ensures that before δt ∧ τ , the particle does not leave the interval [x, z].

For a multi-dimensional generalizations of these simulations, using squares, rectangles or spheres, see [START_REF] Campillo | A Monte Carlo method without grid for a fractured porous domain model[END_REF][START_REF] Deaconu | Hitting time for Bessel processes -walk on moving spheres algorithm[END_REF][START_REF] Deaconu | A random walk on rectangles algorithms[END_REF][START_REF] Milstein | Simulation of a space-time bounded diffusion[END_REF].