Background

Small-on-top theory: In a binary system, small size colloids prefer to accumulate on the external skin layer in the drying process

•▪

  First stage < 6 months • Follow up < 36 months https://www.adroitmarketresearch.com/industry-reports/baby-infant-formula-Particle shape depends on WPI% R with a critical percentage in the range 20% to 50% ▪ WPI colloids have a significant impact on interface rheological propertiesTo understand the skin layer formation in the WPI/NPC droplet drying, and consequently the evolution of surface mechanical properties

•

  Drying kinetics curves correspond to the morphology evolution stages • WPI% R =50% curve close to pure WPI How protein physico-chemical properties affect skin structure and mechanical properties? duration of surface buckling or border delamination • WPI% R ≤20%: deformable skin layer → buckling • WPI% R ≥50%: rigid skin layer → delamination WPI and NPC segregated during sol-gel transition → WPI accumulates on the outer layer ???? showed good agreement with our experiment: in the skin layer, WPI colloid are localized on the top of NPC micelles Main condition for particles segregation Small-on-top theory ------------ZJD model Zhou et al. -2017

➢

  In droplet of mixes, the WPI% R leads the morphology showing two shape categories: surface buckling (WPI% R ≤20%) or rigid convex shell (WPI% R ≥50%)➢ The drying kinetics and the drying time of mixes samples are governed by whey proteins if the WPI% R ≥ WPI% Rc ( 20%≤ WPI% Rc ≤50%) ➢ In the mixed sample, WPI molecules tend to accumulate on the particle outer layer (SEM image, confirmed by ZJD modelling) ➢ Whey proteins have a significant influence on the interfacial rheology properties if WPI% R

	Summary
	WPI%	20%	50%	80%
	R
	Ψ	0.01	0.03	0.05
	WPI
	29/06/2021		Material &	p. 12
			Background	Results	Summary
				methods

≥WPI% Rc 29/06/2021 Background Material & methods Results Summary p. 13 Merci 谢谢 Thanks for your attention 29/06/2021


Lea Salelles

e. early stages of gastric digestion)

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

p. 1

Juin 2021/ Journées scientifiques ED EGAAL/Salelles Léa Pepsin, secreted by the stomach's parietal cells, is the first protease encountered during digestion. The stomach slowly acidifies the food intake and it is generally said that pepsin is active at acidic pH only (pH ≤ 4), with an optimum around pH 2 1 . However, the impact of pH on gastric proteolysis has not been studied on complex mixes of alimentary proteins, only on few neat globular proteins. The goal of this study was to compare the kinetics of peptic hydrolysis of dairy caseins and egg white proteins as a function of both pH and in-vitro gastric digestion time. • The hydrolysis kinetics could be fairly modelled by power law, except for caseins pH 4 and 5 (Fig. 1).

Poteolysis kinetics by pepsin during gastric in vitro digestion at various pH.

• Final degrees of hydrolysis (at 2h) and initial rates of reaction both follow a similar trend as a function of pH (Fig. 23).

• Pepsin activity on egg white does not show a bell-shaped curve as a function of pH (Fig. 2), contrarily to what is commonly stated from the results observed with haemoglobin.

• Pepsin shows an optimum pH at 3-4 and remain active up to pH 6-7 on caseins (Fig. 2-3), that is the initial pH of cow milk ~6,7 3 and of gastric digestion.

Fig. 1 :

 1 Fig.1 : Degree of hydrolysis of egg-white and caseins by pepsin as a function of time, measured by the pH-STAT method. Data represent means ± standard deviation on 3 replicates. Full lines represent the model fittings of experimental values.

Fig. 2 :

 2 Fig. 2 : Degree of hydrolysis at 2 h of digestion

  et al. (2006); Cornish-Bowden Biochem. J. 113 (1969) 2. Brodkorb et al. Nat. Protoc. 14 (2019) ; Minekus et al.(2014) 3. Dalgleish & Corredig, 2012

  The kinetics of gastric proteolysis is strongly influenced by the dietary proteins, with a remarkably high pepsin activity on dairy caseins for pH ≥ 4 (i.e. early stages of gastric digestion)

		14	Power law model 𝐃𝐇 = 𝒂 * 𝑻𝒊𝒎𝒆 𝒃	Egg-white		Power law model 𝐃𝐇 = 𝒂 * 𝑻𝒊𝒎𝒆 𝒃	Caseins
	Hydrolysis (%)	12 8 10	pH a 2 0.33 3 0.22 4 0.15 5 0.08	b 0.41	R2 >0.99 >0.99 >0.99 0.96			pH a 2 0.88 3 1.09 4 1.13 5 1.00 6 0.46	b 0.23	R2 >0.99 >0.99 0.91 0.74 0.97
	of	6
	Degree	4
		2
		0
		0		2000	4000 Time (s)	6000	8000
			pH 2		pH 3	pH 4	pH 5	pH 6	pH 2	pH 3	pH 4	pH 5	pH 6

Léa Salelles, Steven Le Feunteun, Juliane Floury lea.salelles@inrae.fr ,steven.le-feunteun@inrae.fr STLO, INRA, AGROCAMPUS OUEST, Rennes, France p. 2 Juin 2021/ Journées scientifiques ED EGAAL/Salelles Léa


Elodie Louisy

Gaelle Fontaine

Valérie Gaucher

Fanny Bonnet

Grégory Stoclet


Henri Benisty
email: henri.benisty@institutoptique.fr

Vincent Brac

Abderahim Ramdane

Anatole Lupu

Parity-time Symmetric gratings in 1550 nm Distributed-Feedback lasers diodes: insight on device design rules

We build up on former results from our work on parity-time symmetric gratings implemented in 1550 nm distributed feedback laser diodes to address the design issues raised by the first observed trends. These laser diodes are of "complex-coupled" nature, with modulations of both real and imaginary part of effective index, with relative phase π related to the parity-time symmetry. The unidirectionality of the photonic behavior in reflection mode is dependent on the level of extra losses incurred by the metallic grating used to implement fixed loss modulation onto a nearly uniform gain, either as a first order or as a third order grating. The observed behavior suggests that facets play a large role in setting the desired "unidirectional" lasing operation points, with preferential emission on one side. We explore this issue, of generic interest for the coupling of parity-time symmetric structures to open space.

Introduction

Parity-time symmetry (PTS) has been a remarkably successful concept in optics [START_REF] Bender | Real spectra in non-Hermitian Hamiltonians having PT symmetry[END_REF][START_REF] Rüter | Observation of parity-time symmetry in optics[END_REF][START_REF] El-Ganainy | Theory of coupled optical PTsymmetric structures[END_REF][START_REF] Makris | Beam dynamics in PT symmetric optical lattices[END_REF][START_REF] Musslimani | Optical solitons in PT periodic potentials[END_REF][START_REF] Guo | Observation of PT-Symmetry Breaking in Complex Optical Potentials[END_REF][START_REF] Longhi | PT-symmetric laser absorber[END_REF][START_REF] Chong | Coherent perfect absorbers: time-reversed lasers[END_REF][START_REF] Feng | Single-mode laser byparity-time symmetry breaking[END_REF][START_REF] Hodaei | Parity-time-symmetric microring lasers[END_REF][START_REF] Peng | Loss-induced suppression and revival of lasing[END_REF][START_REF] Feng | Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies[END_REF][START_REF] Ma Miri | Large area single-mode parity-timesymmetric laser amplifiers[END_REF][START_REF] Hodaei | PTsymmetric large area single mode DFB lasers[END_REF][START_REF] Gu | Parity-time-symmetric circular Bragg lasers: a proposal and analysis[END_REF]. Blended with periodicity, it naturally introduces an asymmetry or "unidirectionality" in the way dispersion relations are considered [START_REF] Kulishov | Nonreciprocal waveguide Bragg gratings[END_REF][START_REF] Lin | Unidirectional invisibility induced by P T-symmetric periodic structures[END_REF][START_REF] Jalas | What is-and what is not-an optical isolator[END_REF]. The physics and optoelectronics communities have learnt from several pioneers in the last few decades how interesting it was to push periodicity and symmetries to new limits, (notably through metamaterials, active nanophotonics, and chirality, see e.g. Refs. [START_REF] Droulias | Chiral Metamaterials with PT Symmetry and Beyond[END_REF][START_REF] D'amico | Nanoscale quantum optics[END_REF][START_REF] Zhao | Chiral metamaterials: retrieval of the effective parameters with and without substrate[END_REF]). This is basically because a modulation such as the modulation of the complex refractive effective index in a distributed-feedback (DFB) laser diode takes in an ideal limit the simple form n eff (z) = n eff,0 + δn exp(iKz), with real part δn cos(Kz) and imaginary part δn sin(Kz) in quadrature, causing a scattering channel only for +K and none for -K, unlike classical real-index grating for instance [START_REF] Yariv | Periodic structures for integrated optics[END_REF][START_REF] Santis | High-coherence semiconductor lasers based on integral high-Q resonators in hybrid Si/III-V platforms[END_REF].. It must be said that the idea of complex-coupling in DFB gratings was pervasive in the late stage of DFB developments in the 90s [START_REF] Kim | Analysis of single-mode yields above threshold for complex-coupled distributed feedback lasers with asymmetric facet reflectivities[END_REF][START_REF] Champagne | Global and local effects in gain-coupled multiple-quantum-well DFB lasers[END_REF][START_REF] Baets | On the distinctive features of gain coupled DFB lasers and DFB lasers with second-order grating[END_REF][START_REF] Luo | Purely gain-coupled distributed feedback semiconductor lasers[END_REF][START_REF] Kapon | The effect of complex coupling coefficients on distributed feedback lasers[END_REF][START_REF] Morthier | Improved performance of AR-coated DFB lasers by the introduction of gain coupling[END_REF][START_REF] Cardimona | Dephased index and gain coupling in distributed feedback lasers[END_REF][START_REF] Boucher | Quasi-periodic complex-coupled distributed-feedback structures with an exponential-like gradient of coupling[END_REF]. Gain modulation at those times was envisioned as stemming from the periodic etching of multi-quantum-well (MQW) stacks [START_REF] Champagne | Global and local effects in gain-coupled multiple-quantum-well DFB lasers[END_REF][START_REF] Luo | Purely gain-coupled distributed feedback semiconductor lasers[END_REF], but the material issues and device performance penalties incurred by such a choice made the outcome not really profitable, and more fundamentally, such strategies set the choice of the relative phase of real and imaginary parts of index to 0 or π, leaving aside the PTS window of opportunity, already envisioned in those years but without the current terminology [START_REF] Boucher | Quasi-periodic complex-coupled distributed-feedback structures with an exponential-like gradient of coupling[END_REF][START_REF] Boucher | Non-reciprocal effects of complex-coupled distributed-feedback structures resulting from the phase difference between the coupling constants[END_REF].

We recently attempted to implement such PTS DFB lasers diodes [START_REF] Brac De La Perrière | Electrically injected parity-time symmetric distributed feedback laser diodes (DFB) for telecom applications[END_REF]. Our aim is to use this fundamental PTS concept and the associated "unidirectional" reflection to address an important issue of real-life telecommunication laser diode systems in fiber networks: the sensitivity of laser coherence to optical feedback, cast as an "ORL" or "optical return loss" [START_REF] Azouigui | Optical feedback tolerance of quantum dot and quantum dash based semiconductor lasers operating at 1.55 µm[END_REF][START_REF] Ke | Parity-time symmetric complex-coupled distributed feedback laser with excellent immunity to external optical feedback[END_REF]. If, for some reason along a fiber link, sizable reflections arise, this means return with rather high signal and small "return losses". Such a delayed feedback is well-known to compromise the laser spectral coherence. As coherent detection is having, since a decade or so, an edge in terms of pushing the bandwidth capability of optical networks to everincreasing heights (to confront world demand for bandwidth), the attainment of optimal coherence is a prerequisite. The issue is mostly tackled through expensive optical isolators. These devices fundamentally break the transmission matrix reciprocity thanks to magnetooptics. But their cost is a burden on the broader deployment of high-end DFB lasers beyond the core facilities and infrastructures, closer to homes or in the wide and demanding domain of datacom inside data centers (where copper eventually meets its limits at few meters link distances).

Hence, isolator-free solutions are a hallmark of recent device research. Using a fundamental concept to meet this class of goals is our general aim. Our previous work described our first results in that domain. We are going to summarize them briefly in the next section, Sec.2, as the audience is still mostly unaware of this branch of applied research, so as to connect it with the more popular PTS "lab-devices" that have been the center of the focus in the last few years. The issues at stake are the grating order, the loss level, and the impact on spectral purity and on "unidirectionality", as well as the actual impact, in our first tests, on "ORL" level, the key performance, whose values are the subject of regularly updated norms such as IEEE 802.3. In the following, Sec.3, we describe our efforts to understand the combined role of grating characteristics and lasers facets (high or low reflection) in the lasing behavior: since a high-efficiency PTS laser has to operate at moderate level of distributed reflection [START_REF] David | Gain-coupled DFB lasers versus index-coupled and phase shifted DFB lasers: a comparison based on spatial hole burning corrected yield[END_REF], the way the facets influence the feedback appears as more crucial to our scope than usual. The next section, Sec.4, is specifically devoted to the trends in the role of facets for unidirectionality. We combine the simulation results with our experience from previous work to discuss what an optimal facet design could be in section 5. We then conclude the paper.

Experimental findings summary

The InP platform is the reference for 1550 nm telecom laser diode. Few PTS have been investigated on this platform, as it is not the simplest to tailor gain and losses at the wavelength scale. One aspect of good laser diodes is that the active structure is buried, and isolated laterally by epitaxial regrowth of InP of the proper type to channel current from the top electrode to the few-µm-wide MQW active region. The grating generally sits on top of the active layer. A preferred alternative design to access the optical mode's longitudinal landscape periodic modulation is a relatively tall ridge structure, with electrode on top. As show in Fig. 1(a), in such a geometry, the MQW is buried just beneath the bottom of the ridge, and the ridge width modulation does translate into a sizable effective index modulation of the laser optical mode, but only of its real part (neglecting radiation losses for third order gratings). Then, as the optical mode still overlaps the bottom region the "cliffs" of the modulated ridge, metal stripes of the same periodicity (or of different periodicity if a different order is chosen) can be implemented on the adjacent area (as done in Ref. [START_REF] Santis | High-coherence semiconductor lasers based on integral high-Q resonators in hybrid Si/III-V platforms[END_REF] for commercial DFB lasers targeting a market of specialty wavelengths), with controlled relative phase.

The essential test of the benefits of PTS systems can then be investigated by tuning the relative phase of the loss and gain parts, Fig. 1(b). This amounts to choose the phase ϕ in a general formula of the modulation reading δn(z) = δn' cos(Kz)+ iδn"cos(Kz+ϕ). Of course this phase correspond to a physical shift of the fraction (ϕ/2π) of the grating period Λ = 2π/K (≈ m λ laser /2n eff at order m = 1 or 3 here). The general case is a case of "complex-coupling". In the usual notation, the coupling constant κ is essentially the complex amplitude of δn(z), but here its +K and -K Fourier components should be distinguished: within a factor of 2, they are δn'+ iδn"e iϕ and δn'+ iδn"e -iϕ , so that the choice ϕ = π/2 and δn' = δn" >0 cancels the +K Fourier component of κ and doubles its -K one. From a practical perspective, we have three continuous parameters at hand, the two index modulations and a single phase since we reference the phase to the real part of the grating index modulation (in other words, the four real numbers that control the two Fourier components at K and -K are not independent given the phase reference to the real part of the grating index modulation). Only when it comes to the addition of facets at gratings end points does this reference convention entail some consequences. But for the grating alone it is neutral. A large part of the experimental work was to find the proper technological parameters to get a good laser from this structure, with a reasonably accurate knowledge of the coupling constants (±15%), which was checked by several means. A typical difficulty is that the metal and the optical field overlap most at the critical region of the cliff's foot, where a 30 nm change in metal position, easily found due to the delicate process of defining stripes close to cliffs, has a substantial effect on the inferred (calculated) κ values.

Our study in Ref. [START_REF] Brac De La Perrière | Electrically injected parity-time symmetric distributed feedback laser diodes (DFB) for telecom applications[END_REF] compared index-coupled lasers (no metallic grating), gain-coupled lasers (no ridge modulation, only the weaker modulation of real index by the metal survives, but gold at 1550 nm has a small real part of the index), and complex-coupled lasers such as those of Fig. 1. However we did not pursue the goal of pure PTS because it would entail metal losses to become excessive. Hence we explored complex-coupled lasers with sizable asymmetries in the κ value seen by the two counter-propagating waves (say 20%), but far from cancelling one with respect to the other. Still, in a laser, small differences in coupling may have large influence on the lasing mode competition.

Spectra with clear single mode behavior and device-grade characteristics in terms of output power (> 5 mW) were observed. Some devices were tested with external sources to check the existing of a substantial anisotropy in reflection in conditions close to transparency. The first ORL tests were also reported in Ref. [START_REF] Brac De La Perrière | Electrically injected parity-time symmetric distributed feedback laser diodes (DFB) for telecom applications[END_REF], whereby the coherence collapse was tracked by adequate spectroscopic means under increasing return losses from a desktoplength fiber fitted with a variable loop-based attenuator. An interesting proportion of lasers had sufficient feedback immunity to pass the IEEE802.3 stringent level, one of them tolerating up to -10 dB return. However the data were somehow dispersed and hard to exploit in a clear-cut manner. Fig. 2(a) shows for instance that complex-coupled lasers provide very similar spectra on both sides, in spite of the grating asymmetry. For a weak κ value, and henceforth a weak κL value, this is somehow possible, because the two counter-propagating waves are large throughout the device. But it suggests that the extra role of the facets is predominant in these conditions.

A similar conclusion can be drawn from the set of L-I curves from lasers of a same bar, Fig. 2(b). These lasers were designs with increasing phases ϕ, by steps of π/4: this means relative shifts of the dielectric and metallic gratings (Fig. 1a) of Δz = Λ/8 (Δz ≈ 80 nm for third-order grating) between the two gratings, from one laser to the adjacent one in the bar. The two PTS lasers manage to get more light on the (conventional) SD side than on the other side denoted SU. But the overall trend is not marked to pinpoint a genuinely major role of PTS-related effects. ) applies to all data. The L over R preference for the PTS cases is not showing up strongly, in spite of the existing directionality in reflection [START_REF] Brac De La Perrière | Electrically injected parity-time symmetric distributed feedback laser diodes (DFB) for telecom applications[END_REF], with interplay of the laser facets being the likely cause.

In the rest of the paper, we attempt to explain these behaviors by invoking the non-trivial role of the facets in these configurations. The issue is somewhat reminiscent of how the surface cut of metamaterials with underlying periodicity influences the observed properties. The facet reflection (~0.5 in amplitude with antireflection [AR] or high-reflection [HR] coating) is of the same order of magnitude as the single pass reflection from the grating, hence its phase can play a large role.

Such physics already exists for simple DFB index-coupled lasers, but the issue was different: the modal structure of the gain grants two modes with similar gain thresholds on each side of the grating (for first order operation). Then the facets are mainly lifting the degeneracy, and devices are eventually sorted to discard those with unfit SMSR (secondary mode suppression ratio) in terms of acceptable products. Here, we have gratings with a more clearly asymmetric behavior from the start. The facets thus tend to spoil their quality. And we cannot play much the game of quarter-phase shifted DFB with large κL (that naturally grants quite some feedback immunity) because the metal losses would incur a substantial penalty in raising the threshold point. Hence we must clarify the role of facets to assess the viability of our approach in its present form, or inspire solutions that make it good enough for actual devices.

Modeling of various gratings and of the impact of variable facets

An important result of the study reported in [START_REF] D'amico | Nanoscale quantum optics[END_REF] is that despite the fact that the asymmetry properties of the complex profile grating are spoiled by the reflection of the facets, the combined action of the index-and loss-modulated gratings is nonetheless very beneficial for improving the single-frequency operation of DFB lasers. The reasons for this behavior can be understood by considering the threshold gain required to achieve a lasing emission for the different longitudinal modes inside the DFB laser. This threshold gain can be determined by modeling the DFB laser array as a 1D periodic stack of layers of thickness /8 whose refractive index matches the effective index of the Bragg waveguide at the corresponding position along the light propagation axis. This is illustrated in Fig. 3 on the example of a single grating period of length . Mirrors will be added in the next section, but now the termination are facets with air. The relevant longitudinal mode threshold gain g th can classically be determined from the condition of the observation of the first pole of the transmission matrix. To this end we consider realistic DFB lasers examples with Bragg grating modulation parameters corresponding to those determined from experimental data. We thus revisit classical results but with the extra trends dictated by these values. The distribution as function of wavelength of the transmission poles for 1-mm-long conventional DFB laser with modulation of only real part of the refractive index n Re = 0.002 is shown in Fig. 4a. The two lowest poles at g th = 32 dB/cm are corresponding to the onset of a laser emission at the upper and lower edges of the Bragg grating stop-band centered around 1.549 µm. It is obvious that such a dualfrequency operation does not comply with the requirements of telecommunications. Note that the poles located far from the stop band at g th FP =55 dB/cm are those of the Fabry-Perot cavity, i.e. the wave essentially averaging all grating effects, formed by the end facets with 28% reflection. The case of only loss modulated DFB grating with n Im = 0.0004 is illustrated in Fig. 4b. In contrast to the case of real refractive index modulation, this time a clear single frequency behavior with g th Bragg =78 dB/cm is observed. The improvement of the single frequency operation is obtained however at the price of an important penalty on the threshold gain. This is also visible from the gain required for reaching Fabry-Perot cavity poles, which is also strongly increased up to g th FP =122 dB/cm. The situation corresponding to the complex Bragg grating with n Re = 0.002 and n Im =0.0004 is illustrated in Fig. 4c. The Bragg grating pole level is now g th Bragg = 72 dB/cm, i.e. 6 dB/cm lower with respect to the case of only loss modulated grating. Note that the gain level of Fabry-Perot cavity poles is the same (g th FP =122 dB/cm).

An obvious question is whether it is possible to further reduce the Bragg-grating-related gain-level while keeping single-frequency behavior and high level of gain discrimination with respect to the poles of the Fabry-Perot cavity? The solution to this problem is to use the unidirectional reflection properties characteristic for Parity-Time symmetric gratings. The issue is that in our case on the one hand n Im ≪ n Re , as is clear from the discussion above further increase of n Im would be highly detrimental for the threshold gain, and on the other hand, due to the relatively small Bragg grating contrast (n Re = 0.002), the reflection from the grating is strongly spoiled by the facets reflection.

To solve these problems, it is necessary to introduce an additional asymmetry in the system by making one facet highly reflective (HR) and another one with anti-reflection (AR) coating while keeping the relevant product R HR R AR identical to that of an uncoated DFB. When the asymmetry of the facets reflection coincides with that of the Bragg grating (in the sense of favoring stronger reflection on the same side), this results in a strong decrease of the threshold gain. The corresponding situation is illustrated in Fig. 4d with a less-demanding value R AR = 0.08. As can be seen, g th Bragg is reduced down to 30 dB/cm and this pole is separated by at least 70 dB/cm from the nearest poles.

This result implies very stable single-frequency operation of the DFB laser, but it was obtained for a particular example of a 1-mm-long DFB laser with an integer number of grating periods. This choice could be far from neutral, as the interplay of facet reflectivity and grating reflectivity is strong, both having similar orders of magnitude. However, in practice, it would be very difficult to control the exact grating termination with the corresponding precision of /30. Therefore, it is important to evaluate the robustness of single-frequency operation compared to random phase termination of the ends of AR and HR coated Bragg gratings. This can be done by introducing the facet termination related phase defined as the ratio of the incomplete terminal period of length L with respect to the Bragg period  (see Fig. 5), then determining the lowest g th Bragg in the spectrum for each pair of values of AR and HR facet phases. The results of such a procedure are shown in Fig 6 .  In the case when the grating asymmetry coincides with that of complex Bragg grating (AR and obviously laser output on the low reflection side of the grating), there are large white color regions visible in Fig. 6a corresponding to the relatively low threshold gain g th Bragg < 40 dB/cm. The remarkable feature is that these regions are rather weakly sensitive to the exact phase (position) of the AR-coated end facet termination. In contrast, they are confined to specific areas of the HR coated end facet termination, namely phases around ±45° or ±50°. The total yield corresponding to the white colored regions can therefore be estimated to at least 25% for a fabrication with a random facet phase.

The situation is very different in the case when the asymmetry of AR and HR facets reflection is opposite with respect to that of the Bragg grating. This can be understood physically, as the wave amplitudes are not fitting those of the Bloch eigenmodes of the grating on account of the HR reflection in particular when the local mode would seek low reflection. The unidirectional properties similar to those of loss-modulated Bragg grating are thus showing up, despite the fact that the modulation of the imaginary part of index profile n Im = 0.0004 is five times lower with respect to the modulation of the real part n Re = 0.002. Hence, the interplay of facets and gratings cannot be devised in "rather additive" terms. But fortunately, the good solution still offers a tractable landscape allowing for random facet fabrication with acceptable yields. In addition to the low threshold gain level, AR/HR coated lasers also exhibit a high level of gain discrimination between the lowest and next pole. The corresponding situation is shown in Fig. 6b. These results ensure a stable single-frequency lasing with a high SMRS level and also provide the expectation of a better tolerance with respect to optical return loss. Fig. 6. Color maps showing variation of the threshold gain g th Bragg and threshold gain difference as a function of both the AR and the HR facets phase, when the asymmetry of AR and HR coatings reflectivity coincides with that of complex profile Bragg grating (HR on high reflective side of grating). (a) Gain threshold map with contours. Note the weak dependence on the AR facet for the low-threshold regime of HR phase around 0°/360° (for the opposite choice, AR and HR anti-coinciding, the gain threshold map is mostly above 85 dB/cm). (b) Map of the gain difference with their contours. The dashed line is the 40 dB gain-threshold contour of (a), delimiting the low threshold area. The secondary mode suppression ratio (SMSR) seems also good in those areas, essentially above 50 dB; (c) Gain threshold map in the case when the asymmetry of the PTS Bragg grating reflection is opposite to the asymmetry of AR/HR coating.

Conclusion

Laser diode with PTS symmetry in the DFB grating generalize the constrained complexcoupled DFB laser diodes that was first explored across the 90s. Continued interest for broader operating characteristics of real-life DFB lasers in networks of all kinds is a logical incentive to explore novel solutions such as PTS DFB laser diodes. We have specifically explored the role of grating cut by randomly placed facets in the non-trivial complex-coupled distributed-feedback lasers, based on our previous first experimental results. The interplay of facet and gratings was expected due to the order of magnitude of the reflection and the coherent nature of the interaction, but its physics was not hitherto reported. The issue bears some similarity with the generic issues of cutting into complex metamaterial with elaborate unit cells. We found relevant trends from an adapted version of the classical transmission matrix tool and threshold gain analysis.

The main result here is that a choice of relatively well mastered HR and AR facets can be put to good use to get low-threshold devices with a decent yield (>50%), with only a penalty of the HR facet position. A low threshold of about 32 dB/cm modal gain is predicted for 1mm-long PTS DFB laser diode devices having our measured real and complex grating modulation characteristics with only one fifth of complex modulation. A more complete study of SMSR and immunity to optical feedback should help defining optimal devices. With the trend found here and the rather clear physical basis for their understanding, the desired properties of large SMSR and strong immunity could see their emergence largely favored.

Fig. 1 .

 1 Fig. 1. (a) Sketch of the modulated ridge width that contributes to most of the real index modulation, and of the metallic grating, that is the only contributor the imaginary index modulation. Their relative positions are shifted as indicated. (b) Micrograph (tilted view) of a typical device (3 rd order grating, ca. 720 nm pitch) with the dielectric and metallic grating, and the contact on the ridge middle.

Fig. 2 .

 2 Fig. 2. (a) Typical spectra of a given complex-coupled grating device of length 1960 µm at 200 mA with cleaved facet, taken from both side. The DFB operation does provide single-mode operation and the weak difference shows that this is obtained with a relatively weak κL product (b) L-I data (power vs. intensity) of lasers from the same bar with phases increasing by steps of Λ/8 (equivalently ϕ increasing by steps of π/4); L and R designates the two sides. The data for different phases are shifted along the dashed line. The inner optical power axis (0-12 mW) applies to all data. The L over R preference for the PTS cases is not showing up strongly,

Fig. 3 .

 3 Fig. 3. Schematic of the complex index profile Bragg grating period modeled by the Abelès transmission matrix method.

Fig. 4 .

 4 Fig. 4. Maps of the position of the transmission matrix poles for different types of Bragg grating DFB lasers. a) dielectric grating DFB laser with refractive index modulation n Re =0.002. b) Metallic grating DFB laser with refractive index modulation n Im =0.0004. c) double grating DFB laser with refractive index modulation n Re =0.001 and n Im =0.0004,  Re =0,  Im =. D) AR and HR coated double grating DFB laser with refractive index modulation n Re =0.002 and n Im =0.0004,  Re =0,  Im =3/2. R AR =0.079, R HR =0.999.

Fig. 5 .

 5 Fig. 5. Schematic of the complex index profile Bragg grating modeled by the transmission matrix method. The phases of the two facets (AR and HR) are scanned throughout the 1-period long outer zones, and the study is carried out for either AR/HR and HR/AR choices, given the PTS grating asymmetry. The period is much smaller than the length (Λ≪L).

Acknowledgments

This work has partially been funded by the by the French Agence Nationale pour la Recherche (project "PARTISYMO", contract number N° ANR-18-CE24-0024-01), and partly supported by the French RENATECH network. We are very grateful to Q. Gaimard, G. Aubin, S. Guilet, L. Le Gratiet, K. Merghem, F. Raineri and C. Ulysse for their expert help in many ways. VB acknowledges financial support from EDOM "Ecole Doctorale" through a PhD grant.


Badr Tajini

Hugo Richard

Bertrand Thirion

Functional Magnetic Resonance Imaging data augmentation through conditional ICA

Keywords: Conditional ICA, Data generation, Decoding studies

Advances in computational cognitive neuroimaging research are related to the availability of large amounts of labeled brain imaging data, but such data are scarce and expensive to generate. While powerful data generation mechanisms, such as Generative Adversarial Networks (GANs), have been designed in the last decade for computer vision, such improvements have not yet carried over to brain imaging. A likely reason is that GANs training is ill-suited to the noisy, highdimensional and small-sample data available in functional neuroimaging. In this paper, we introduce Conditional Independent Components Analysis (Conditional ICA): a fast functional Magnetic Resonance Imaging (fMRI) data augmentation technique, that leverages abundant restingstate data to create images by sampling from an ICA decomposition. We then propose a mechanism to condition the generator on classes observed with few samples. We first show that the generative mechanism is successful at synthesizing data indistinguishable from observations, and that it yields gains in classification accuracy in brain decoding problems. In particular it outperforms GANs while being much easier to optimize and interpret. Lastly, Conditional ICA enhances classification accuracy in eight datasets without further parameters tuning.

Introduction

As a non-invasive brain imaging technique, task fMRI records brain activity while participants are performing specific cognitive tasks. Univariate statistical methods, such as general linear models (GLMs) [START_REF] Friston | Analysis of fmri time-series revisited[END_REF] have been successfully applied to identifying the brain regions involved in specific tasks. However such methods do not capture well correlations and interactions between brain-wide measurements. By contrast, classifiers trained to decode brain maps, i.e to discriminate between specific stimulus or task types [START_REF] Loula | Decoding fMRI activity in the time domain improves classification performance[END_REF][START_REF] Shirer | Decoding Subject-Driven Cognitive States with Whole-Brain Connectivity Patterns[END_REF][START_REF] Varoquaux | How machine learning is shaping cognitive neuroimaging[END_REF], take these correlations into account. The same framework is also popular for individual imaging-based diagnosis.

However, the large sample-complexity of these classifiers currently limits their accuracy. To tackle this problem, data generation is an attractive approach, as it could potentially compensate for the shortage of data. Generative Adversarial Networks (GANs) are promising generative models [START_REF] Goodfellow | Generative adversarial nets[END_REF]. However, GANs are ill-suited to the noisy, high-dimensional and small-sample data available in functional neuroimaging. Furthermore the training of GANs is notoriously unstable and there are many hyper-parameters to tune.

In this work, we introduce Conditional ICA: a novel data augmentation technique using ICA together with conditioning mechanisms to generate surrogate brain imaging data and improve image classification performance. Conditional ICA starts from a generative model of resting state data (unconditional model), that is fine-tuned into a conditional model that can generate task data. This way, the generative model for task data benefits from the abundant resting state data and can be trained with few labeled samples. We first show that the generative model of resting state data shipped in Conditional ICA produces samples that neither linear nor non-linear classifiers are able to distinguish. Then we benchmark Conditional ICA as a generative model of task data against various augmentation methods including GANs and conditional GANs on their ability to improve classification accuracy on a large task fMRI dataset. We find that Conditional ICA yields highest accuracy improvements. Lastly, we show on 8 different datasets that the use of Conditional ICA results in systematic improvements in classification accuracy ranging from 1% to 5%. Fig. 1: Conditional ICA approach. Our method aims to generate surrogate data from Task and Rest fMRI data by synthesizing statistical maps that qualitatively fit the distribution of the original maps. These can be used to improve the accuracy of machine learning models that identify contrasts from the corresponding brain activity patterns.

Methods

Notations We write matrices as bold capital letters, vectors as small bold letters. X † refers to the Moore-Penrose pseudo inverse of matrix X, , tr(X) refers to the trace of matrix X and I k refers to the identity matrix in R k,k . 0 k refers to the null vector in R k . Spatial Dimension reduction. The outline of the proposed approach is presented in Fig. 1. While brain maps are high-dimensional, they span a smaller space than that of the voxel grid. For the sake of tractability, we reduce the dimension of the data by projecting the voxel values on the high-resolution version of the Dictionaries of Functional Modes DiFuMo atlas [START_REF] Dadi | Fine-grain atlases of functional modes for fmri analysis[END_REF], i.e. with p = 1024 components. The choice of dimension reduction technique generally has an impact on the results. However we consider this question to be out of the scope of the current study and leave this to future work.

Unconditional generative models (resting state data generation) Given a largescale resting-state dataset X rest in R p,n where n is the number of images (samples) and p = 1024 the number of components in the atlas, let us consider how to learn its distribution. Assuming a Gaussian distribution is standard in this setting, yet, as shown later, it misses key distributional features. Moreover, we consider a model that subsumes the distribution of any type of fMRI data (task or rest): a linear mixture of k ≤ p independent temporal signals. We therefore use temporal ICA to learn a dimension reduction and unmixing matrix W rest ∈ R k,p such that the k sources i.e the k components of S rest = W rest X rest are as independent as possible.

A straightforward method to generate new rest data would be to independently sample them from the distribution of the sources. This is easy because such distribution has supposedly independent marginals. We apply an invertible quantile transform q rest to the sources S rest so that the distribution of z rest = q rest (s rest ) has standardized Gaussian marginals. Since the distribution of z rest has independent marginals, it is given by N (0 k , I k ) from which we can easily sample. As shown later, this approach fails: such samples are still separable from actual rest data.

We hypothesize that this is because independence does not hold, and thus a latent structure among the marginals of the source distribution has to be taken into account. Therefore we assume that the distribution of z rest is given by N (0 k , Λ rest ) where Λ rest is a definite positive matrix. Λ rest can easily be learned from a standard shrunk covariance estimator: Λ rest = Σ rest (1-α)+αtr(Σ rest )I k where α is given by the Ledoit-Wolf formula [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF] and Σ rest is the empirical covariance of Z rest .

Our encoding model for rest data is therefore given by Z rest = q rest (W rest X rest ) and we assume that the distribution of Z rest is N (0 k , Λ k ). The generative model is given by the pseudo inverse of the encoding model:

Xrest = (W rest ) † (q rest ) -1 (ε) where ε ∼ N (0 k , Λ rest ).
Conditional generative models (generative model for task data) While resting state datasets have a large number of samples (10 4 ∼ 10 5 ), task datasets have a small number of samples (10 ∼ 10 2 ). As a result, there are too few samples to learn high quality unmixing matrices. Therefore, using the unmixing matrix W rest learned from the resting state data, we rely on the following nonlinear generative model for brain maps in a certain class c:

x c = (W rest ) † q -1 (ε) (1)
with ε ∼ N (µ c , Λ).

In order to maximize the number of samples used to learn the parameters of the model, we assume that the quantile transform q and the latent covariance Λ do not depend on the class c. However, the mean µ c , that can be learned efficiently using just a few tens of samples, depends on class c. An overview of our generative method is shown in Fig. 2. Fig. 2: Conditional ICA approach in depth. The approach proceeds by learning a temporal ICA of rest data X rest ∈ R p,n , resulting in independent sources and unmixing matrix W rest ∈ R k,p . Applying the unmixing matrix to the task data, we obtain samples in the source space S task ∈ R k,n . Afterwards, we map S task to a normal distribution, yielding Z task ∈ R k,n . Then, we estimate the covariance Λ ∈ R k,k (all classes are assumed to have the same covariance) and the class-specific means µ 1 , . . . , µ C ∈ R k according to Ledoit-Wolf's method. For each class c, we can draw random samples Ztask c ∈ R k,n fakes from the resulting multivariate Gaussian distribution N (µ c , Λ) and obtain fake data Xc ∈ R p,n fakes by applying the inverse quantile transform and re-mixing the data using the pseudo inverse of the unmixing matrix. We append these synthetic data to the actual data to create our new augmented dataset on which we train classifiers.

Related work

In image processing, data augmentation is part of standard toolboxes and typically includes operations like cropping, rotation, translation. On fMRI data these methods do not make much sense as brain images are not invariant to such transformations. More advanced techniques [START_REF] Zhuang | Fmri data augmentation via synthesis[END_REF] are based on generative models such as GANs or variational auto-encoders [START_REF] Kingma | Auto-encoding variational bayes[END_REF]. Although GAN-based method are powerful they are slow and difficult to train [START_REF] Arjovsky | Wasserstein GAN[END_REF]. In appendix Table 7, we show that Conditional ICA is several orders of magnitude faster than GAN-based methods.

Our method is not an adversarial procedure. However it relates to other powerful generative models such as variational auto-encoders [START_REF] Kingma | Auto-encoding variational bayes[END_REF] with which it shares strong similarities. Indeed the analog of the encoding function in the variational auto-encoder is given by e(x) = Λ -1 2 q(W rest x) in our model and the analog to the decoding function in the variational auto-encoder is given by d(z) = (W rest ) † q -1 (Λ 1 2 z) in our model. As in the variational auto-encoder, e approximately maps the distribution of the data to a standardized Gaussian distribution, while the reconstruction error defined by the difference in l2 norm d(e(x)) -x 2 2 must remain small. Another classical generative model related to ours is given by normalizing flows [START_REF] Kobyzev | Normalizing flows: An introduction and review of current methods[END_REF]. When W rest is squared (no dimension reduction in ICA), the decoding operator d is invertible (its inverse is e) making our model an instance of normalizing flows. A great property is thus the simplicity and reduced cost of data generation.

Software tools

We use the FastICA algorithm [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] to perform independent component analysis and use k = 900 components. We use Nilearn [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF] for fMRI processing, Scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], Numpy [START_REF] Harris | Array programming with numpy[END_REF] and Scipy [START_REF] Virtanen | SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF] for data manipulation and machine learning tools.

Code availability Our code is available on the following git repository: https: //github.com/BTajini/augfmri.

Experiments

Dataset, data augmentation baselines and classifiers used

The unmixing matrices are learned on the rest HCP dataset [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF] using 200 subjects. These data were used after standard preprocessing, including linear detrending, band-pass filtering ([0.01, 0.1]Hz) and standardization of the time courses. The other 8 datasets [START_REF] Orfanos | The brainomics/localizer database[END_REF][START_REF] Pinel | Genetic and environmental contributions to brain activation during calculation[END_REF][START_REF] Pinel | Genetic and environmental contributions to brain activation during calculation[END_REF][START_REF] Pinel | The functional database of the archi project: Potential and perspectives[END_REF][START_REF] Pinel | Fast reproducible identification and large-scale databasing of individual functional cognitive networks[END_REF][START_REF] Poldrack | A phenome-wide examination of neural and cognitive function[END_REF][START_REF] Shafto | The cambridge centre for ageing and neuroscience (cam-can) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF][START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF] are obtained from the Neurovault repository [START_REF] Gorgolewski | Neurovault. org: a web-based repository for collecting and sharing unthresholded statistical maps of the human brain[END_REF]. The classes used in each dataset correspond to the activation maps related to the contrasts (such as "face vs tools") present in the set of tasks of each dataset. Details are available in appendix Table 3.

We consider 5 alternative augmentation methods: ICA, Covariance, ICA + Covariance, GANs and CGANs. When no augmentation method is applied we use the label Original.

The ICA covariance method applies ICA to X task to generate unmixing matrices W task and sources S task = W task X task . To generate a sample xc from class c, we sample independently from each source restricted to the samples of class c yielding stask The ICA + Covariance method combines the augmentation methods ICA and Covariance: samples are drawn following the ICA approach, but with some additive non-isotropic Gaussian noise. As in ICA we estimate W task and S task from X task via ICA. Then we consider R task = X task -W task S task and estimate the covariance Σ R of R task via LedoitWolf's method. We then generate a data sample xc from class c as with ICA and add Gaussian noise ñ ∼ N (0, Σ R ). Samples are thus generated as xc + ñ.

The GANs (respectively CGANs) method use a GANs (respectively CGANs) to generate data. The generator and discriminator have a mirrored architecture with 2 fully connected hidden layer of size (256 and 512). The number of epochs, batch size, momentum and learning rate are set to 20k, 16, 0.9, 0.01 and we use the Leaky RELU activation function.

We evaluate the performance of augmentation methods through the use of classifiers: logistic regression (LogReg), linear discriminant analysis with Ledoit-Wold estimate of covariance (LDA) perceptron with two hidden layers (MLP) and random forrests (RF). The hyper-parameters in each classifier are optimized through an internal 5-Fold cross validation. We set the number of iterations in each classifier so that convergence is reached. The exact specifications are available in appendix Table 4.

Distinguish fake from real HCP resting state data

This experiment is meant to assess the effectiveness of the data augmentation scheme in producing good samples. Data augmentation methods are trained on 200 subjects taken from HCP rest fMRI dataset which amounts to 960k samples (4800 per individual). Then synthetic data corresponding to 200 synthetic subjects are produced, yielding 960k fake samples and various classifiers are trained to distinguish fake from real data using 5-Fold cross validation. The cross-validated accuracy is shown in Table 1. Interestingly, we observe a dissociation between linear models (LogReg and LDA) that fail to discriminate between generated and actual data, and non-linear models (MLP and RF) that can discriminate samples from alternative augmentation methods. By contrast, all classifiers are at chance when Conditional ICA is used.

Comparing augmentation methods based on classification accuracy on task HCP dataset

In order to compare the different augmentation methods, we measure their relative benefit in the context of multi-class classification. We use 787 subjects from the HCP task dataset that contains 23 classes and randomly split the dataset into a train set that contains 100 subjects and a test set that contains 687 subjects. In each split we train augmentation methods on the train set to generate fake samples corresponding to 200 subjects. These samples are then appended to the train set, resulting in an augmented train set on which the classifiers are 2, show that Conditional ICA always yields a higher accuracy than when no augmentation method is applied. The gains are over 1% on all classifiers tested. By contrast, ICA+Covariance and ICA lead to a decrease in accuracy while the Covariance approach leads to non-significant gains. To further investigate the significance of differences between the proposed approach and other state-of-the-art methods, we perform a t-test for paired samples (see Table 5 in appendix). Most notably, the proposed method performs significantly better than other data augmentation techniques. Given the large size of the HCP task data evaluation set, this significance test would demonstrate that the gains are robust. Note that the random forest classifier is not used in this experiment as it leads to a much lower accuracy than other methods. For completeness, we display the results obtained with random forest in appendix Table 6.

Visualization of fake examples produced by GANs, CGANs and Conditional ICA are available in appendix Figure 5.

Gains in accuracy brought by conditional ICA on eight datasets.

In this experiment we assess the gains brought by Conditional ICA data augmentation on eight different task fMRI dataset. The number of subjects, classes and the size of the training and test sets differ between dataset and are reported in appendix Table 4. The rest of the experimental pipeline is exactly the same as with the HCP task dataset. We report in Fig. 3 the cross-validated accuracy of classifiers with and without augmentation. We notice that the effect of data augmentation is consistent across datasets, classifiers and splits, with 1% to 5% net gains. An additional experiment studying the sensitivity of Conditional ICA to the number of components used is described in appendix Figure 4.

Models

Discussion & Future work

In this work we introduced Conditional ICA a fast generative model for rest and task fMRI data. It produces samples that cannot be distinguished from true actual rest by linear as well as non-linear classifiers, showing that it captures higher-order statistics than naive ICA-based generators. When Conditional ICA is used as a data augmentation method, it yields consistent improvement in classification accuracy: on 8 tasks fMRI datasets, we observe increase in accuracy between 1% and 5% depending on the dataset and the classifier used. Importantly, this performance was obtained without any fine-tuning of the method, showing its reliability. One can also notice that our experiments cover datasets with different cardinalities, from tens to thousand, and different baseline prediction accuracy. It is noteworthy that Conditional ICA is essentially a linear generative model with pointwise non-linearity, which makes it cheap, easy to instantiate on new data, and to introspect.

The speed and simplicity of Conditional ICA as well as the systematic performance improvement it yields, also make it a promising candidate for data augmentation in a wide range of contexts. Future work may focus on its applicability to other decoding tasks such as the diagnosis of Autism Spectrum Disorder (ASD) [START_REF] Dvornek | Identifying autism from resting-state fmri using long short-term memory networks[END_REF][START_REF] Eslami | Asd-diagnet: a hybrid learning approach for detection of autism spectrum disorder using fmri data[END_REF][START_REF] Eslami | Auto-asd-network: a technique based on deep learning and support vector machines for diagnosing autism spectrum disorder using fmri data[END_REF] or Attention-Deficit/Hyperactivity Disorder detection (ADHD) [START_REF] Mao | Spatiotemporal deep learning method for adhd fmri classification[END_REF]. Other extensions of the present work concern the adaptation to individual feature prediction (e.g. age) where fMRI has shown some potential.

Dataset

Subjects, Classes Train/Test Neurovault collection hcp [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF] 787, 23 100/687 collection 4337 cam-can [START_REF] Shafto | The cambridge centre for ageing and neuroscience (cam-can) study protocol: a cross-sectional, lifespan, multidisciplinary examination of healthy cognitive ageing[END_REF] 605, 5 100/505 collection 4342 brainomics [START_REF] Orfanos | The brainomics/localizer database[END_REF] 94, 19 50/44 collection 4341 archi [START_REF] Pinel | The functional database of the archi project: Potential and perspectives[END_REF] 78, 30 40/38 collection 4339 la5c [START_REF] Poldrack | A phenome-wide examination of neural and cognitive function[END_REF] 191, 24 100/91 collection 4343 pinel2012archi [START_REF] Pinel | The functional database of the archi project: Potential and perspectives[END_REF] 76, 10 40/36 collection 1952 pinel2009twins [START_REF] Pinel | Genetic and environmental contributions to brain activation during calculation[END_REF] 65, 12 35/30 collection 1952 pinel2007fast [START_REF] Pinel | Fast reproducible identification and large-scale databasing of individual functional cognitive networks[END_REF] 133, 10 70/63 collection 1952

c

  and mix the data: xc = (W task ) † stask c . The Covariance method generates a new sample of synthetic data in class c by sampling from a Multivariate Gaussian with mean µ c and covariance Σ, where µ c is the class mean and Σ is the covariance of centered task data estimated using Ledoit-Wolf method. In brief, it assumes normality of the data per class.

Fig. 3 :

 3 Fig. 3: Accuracy of models for eight multi-contrast datasets. Cross validated accuracy of two linear (LDA and LogReg) and one non-linear classifier (MLP) with or without using data augmentation. The improvement yielded by data augmentation is displayed in red. Black error bars indicate standard deviation across splits while white error bars indicate standard deviation across splits with no augmentation.

Fig. 4 :

 4 Fig.4: Accuracy of augmented discriminative models when varying k. We use 100 train subjects from the HCP task dataset to train Conditional ICA with k components and generate 200 fake subjects. The classifiers are trained on the train and fake subjects and tested on the left-out 687 subjects. We repeat the procedure for various values of k using 5 random splits per value and report the mean accuracy across splits as a function of k. The dotted line represents the number of components that has been used in our experiments (k = 900).

Fig. 5 :

 5 Fig. 5: Data generation visualization. Visualization of real (Original) and synthetic brain maps from three generation methods: Conditional ICA (Ours), CGANs and GANS. Three cognitive tasks are shown (reward, punish and math).

Table 1 :

 1 Distinguish fake from real HCP resting state data We use HCP resting state data from n = 200 subjects (960k samples) and produce an equally large amount of fake data (960k samples) using data augmentation methods. The table shows the 5-fold cross validated accuracy obtained with various classifiers. When Conditional ICA is used, all classifiers are at chance.

	Models	LDA LogReg Random Forest MLP
	ICA	0.493 0.500	0.672	0.697
	Covariance	0.473 0.461	0.610	0.626
	ICA + Covariance 0.509 0.495	0.685	0.706
	GANs [8]	0.501 0.498	0.592	0.607
	CGANs [17]	0.498 0.493	0.579	0.604
	Conditional ICA 0.503 0.489	0.512	0.523

trained. Results, displayed in Table

Table 2 :

 2 Comparing augmentation methods based on classification accuracy on task HCP dataset We compare augmentation methods based on the classification accuracy (Acc), precision (Pre) and recall (Rec) obtained by 2 linear classifiers (LDA and LogReg) and one non-linear classifier trained on augmented datasets on HCP Task fMRI data. We report the mean accuracy across 5 splits.

			LDA			LogReg			MLP
		Acc	Pre	Rec	Acc	Pre	Rec	Acc	Pre	Rec
	Original	0.893 0.889 0.891 0.874 0.869 0.873 0.779 0.782 0.778
	ICA	0.814 0.809 0.813 0.840 0.836 0.839 0.803 0.805 0.802
	Covariance	0.895 0.894 0.895 0.876 0.877 0.875 0.819 0.823 0.820
	ICA + Covariance 0.816 0.811 0.812 0.840 0.839 0.840 0.815 0.819 0.814
	GANs [8]	0.877 0.871 0.870 0.863 0.865 0.864 0.771 0.779 0.775
	CGANs [17]	0.874 0.875 0.872 0.874 0.872 0.875 0.726 0.731 0.725
	Conditional ICA 0.901 0.903 0.905 0.890 0.888 0.890 0.832 0.835 0.831

Table 3 :

 3 Datasets used in the experiments. The table provides references to the datasets that were used for our experiments, with the number of subjects, the number of classes, the number of subjects in train and test set in each cross validation split and the collection number in Neurovault.

	Methods Optimizer	Hyper-parameters
	LogReg L-BFGS	inverse L 2 regularization strength
		(20 000 iterations)	in {0.0001, 0.001, 0.01, 0.1, 1}
	LDA	Least-squares solver	Estimation of covariance
			using Ledoit-Wolf's method
	RF	-	Default parameters in sklearn
	MLP	Adam	ReLU activation function, fully con-
		(20 000 iterations,	nected architecture with two hidden lay-
		momentum: 0.9,	ers both of size 1024, L2 penalty coeffi-
		batch size: 32,	cient: 10 -5
		learning rate: 0.0001)

Table 4 :

 4 Optimizers and hyper-parameters of classifiers For each classifier, we give the optimization method used as well as the value of hyper-parameters. ≤ 0.01 P ≤ 0.05 ICA + Covariance P ≤ 0.0001 P ≤ 0.0001 P ≤ 0.01 GANs P ≤ 0.001 P ≤ 0.001 P ≤ 0.0001 CGANs P ≤ 0.001 P ≤ 0.01 P ≤ 0.0001

	Models	LDA	LogR	MLP
	Original	P ≤ 0.01 P ≤ 0.05 P ≤ 0.0001
	ICA	P ≤ 0.0001 P ≤ 0.001 P ≤ 0.001
	Covariance	P ≤ 0.05 P

Table 5 :

 5 Statistical significance of the differences between Conditional ICA and other augmentation methods P-values obtained after a t-test for paired samples is performed on the data used to produce Table2.

	Methods	Original ICA COV. ICA + COV. GANs CGANs Cond. ICA
	Random Forest 0.782 0.778 0.780	0.780	0.780 0.779	0.783

Table 6 :

 6 Accuracy of Random Forrest Mean accuracy obtained from data used for producing the Table 2 using a Random Forest classifier according to the chosen augmentation method.

	Methods	Running time (secs)
	GANs	12948.2 (≈ 3,60 hr)
	CGANs	11015.1 (≈ 3,05 hr)
	Conditional ICA	62 s

Table 7 :

 7 Running time. We display the running time of three methods used to generate synthetic data. Conditional ICA is several orders of magnitude faster than GANs or CGANs. In practice the computational overhead induced by Conditional ICA is negligible.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 Framework Program for Research and Innovation under Grant Agreement No 945539 (Human Brain Project SGA3) and the KARAIB AI chair (ANR-20-CHIA-0025-01).


Mathieu Gardrat
email: mathieu.gardrat@entpe.fr

Pascal Pluvinet

Markov based mesoscopic simulation tool for urban freight: SIMTURB

Keywords: Urban freight transport, Model, SIMTURB, Multi-agent simulation, Markov process Avertissement

. With an architecture based on a Markov process, this model offers a complement and to some extent an alternative to multi-agent simulation models, since it makes possible to characterise precisely the routes of freight transport vehicles in a conurbation and characterise the movements of each agent (e.g. vehicle).

Introduction

Urban freight transport modelling has known for more than 20 years radical evolutions and contributions, from statistical models to optimisation techniques, from very basic spreadsheets to very complex multi-agent models. Today, although the vast lack of knowledge has been partially filled, the need for decision support tools is stronger than ever, given the growing complexity of city logistics and needs for relevant indicators. We present in this paper the contribution of Markov processes in this endeavour with a presentation of the SIMTURB model and its results.

The basis of the modelling process of SIMTURB is the disaggregation of a freight distribution matrix into routes via Markov chains [START_REF] Kemeny | Markov chains[END_REF], [START_REF] Janssens | The development of an adapted Markov chain modelling heuristic and simulation framework in the context of transportation research[END_REF]. In this case, the Markov process is used to build routes based on O-D matrices. Thus, routes are Markov chains of varying lengths for which transition matrices are in fact flow distribution matrices. To keep organisational consistency, the goods distribution matrix is subdivided into functionally homogeneous subsets that take into account the vehicle type (3 classes), the management mode (3 classes) and the length of the routes (7 classes). Coupled with geolocalised establishments and a road network, this approach makes it possible to spatially and temporally characterise the vehicle routes, without requiring complex multi-agent approach.

In the first section we present a literature review on urban freight transport modelling and Markov chains. The second section presents the methodology of the SIMTURB model. Section 3 will be dedicated to the confrontation of SIMTURB and results produced by the FRETURB freight generation model. The last section is dedicated to the conclusion and discussion.

Literature review

Urban freight modelling, compared to passenger traffic modelling, is often considered as more a complex and less mature topic in transportation sciences [START_REF] Ogden | Urban goods movement : a guide to policy and planning[END_REF]. Although the first models designed for urban freight were based on four-step models, their efficiency was quickly questioned [START_REF] Toilier | Freight transport modelling in urban areas: The French case of the FRETURB model[END_REF]. More recent researches have therefore led to alternative methods to the classical four-step models. Most urban freight models now integrate the notion of trips as an output to summarise the intricate structure of urban freight flows [START_REF] Comi | Urban freight transport demand modelling: A state of the art[END_REF]. This estimation of trips generated in a city and their distribution can be done using various methods. These methods are usually based on:

 commodity flows converted into delivery generation, which allows building trip sequences [START_REF] Boerkamps | GoodTrip-A new approach for modelling and evaluation of urban goods distribution[END_REF]- [START_REF] Nuzzolo | Urban freight demand forecasting: a mixed quantity/delivery/vehicle-based model[END_REF].  traffic counts, used to deduce vehicle flows or to calibrate the distribution of commodity or trip-based models. These models do not necessarily include freight trip generation but rather a structure for freight trip distribution [START_REF] Cascetta | Estimation of trip matrices from traffic counts and survey data: a generalized least squares estimator[END_REF], [START_REF] List | Estimating truck travel patterns in urban areas[END_REF]. This method implies that the characterisation of trips is limited to O-D couples and types of vehicles, but it fails to bring to light data on types of goods or organisational patterns [START_REF] Bera | Estimation of origin-destination matrix from traffic counts: the state of the art[END_REF].  establishment surveys in order to include trip generation as a primary step for the estimation of vehicle flows [START_REF] Slavin | Enhanced framework for modeling urban truck trips[END_REF]- [START_REF] Holguin-Veras | Freight generation, freight trip generation, and perils of using constant trip rates[END_REF].

Although trips are relevant to quantify the impact of freight flows in cities and provide numerous characteristics on their organisations, this perspective is not suited for all types of urban freight assessments. Indeed, many analyses require specific indicators such as the volume of vehicles or routes (which are indeed related to the number of trips but need further processing to convert trips or commodity flows into vehicles or routes).

To cope with these new needs, multi-agent models in recent years have known a strong development, specifically due to increasing computational capacities [START_REF] Anand | Ontology-based multi-agent system for urban freight transportation[END_REF]. These models are however extremely complex to run and require heavy consumptions of computational resources, furthermore these calculations are often made on a sample of the total amount of the existing agents in a real situation: in order to produce relevant results, multi agent models use a statistically significant sample of agents to represent the reality [START_REF] Ralha | A multi-agent model system for land-use change simulation[END_REF], [START_REF] Maggi | Understanding urban mobility and the impact of public policies: The role of the agent-based models[END_REF]. Effectively, what is considered as a relevant sample or what categories of actors are modelled, is usually built on strong hypotheses and archetypal behaviours. Some agent approaches usually work for aggregated freight demand and do not specify logistics organisations such as routes and are usually limited to some portions of urban activities, typically retailers and transport operators [START_REF] Van Duin | Towards an agent-based modelling approach for the evaluation of dynamic usage of urban distribution centres[END_REF], [START_REF] Marcucci | Urban freight transport and policy changes: Improving decision makers' awareness via an agent-specific approach[END_REF].

It is also important to note that multi-agent models do not necessarily estimate freight flows and routes, but are used for general policy interactions. Although very useful for stakeholders' interactions, this relative simplification of reality fails to take into account the vastness of logistics organisations and extreme but existing transport behaviours [START_REF] Schroeder | A computational framework for a multi-agent simulation of freight transport activities[END_REF]- [START_REF] De Oliveira | Multi-agent modelling approach for evaluating the city logistics dynamic in a vulnerability situation: An exploratory study in Belo Horizonte (Brazil)[END_REF], even if some are based on various types of surveys [START_REF] Le Pira | Towards a decision-support procedure to foster stakeholder involvement and acceptability of urban freight transport policies[END_REF].

If we focus on models simulating transportation behaviours, the main question faced is to build realistic routes. In these models, tours are built either using:

 iterative methods based on incremental trip chain or shipment combination [START_REF] Nuzzolo | A system of models for the simulation of urban freight restocking tours[END_REF], [START_REF] Thoen | Descriptive modeling of freight tour formation: A shipment-based approach[END_REF], [START_REF] Kim | Empirical comparison of tour-and trip-based truck travel demand models[END_REF]. This methods appear relevant but are quite complex to calibrate and implement, using heavy computational resources.  clustering methods which are less resource consuming but tend to over-optimise routes since they tend to concentrate deliveries based on spatial proximity rather than functional characteristics [START_REF] Outwater | Tour based and supply chain modeling for freight: integrated model demonstration in Chicago[END_REF].

Given this state of research, we noticed that very few models use Markov chains in transportation research [START_REF] Levinson | Paving new ground: a Markov chain model of the change in transportation networks and land use[END_REF], [START_REF] Ramezani | On the estimation of arterial route travel time distribution with Markov chains[END_REF], specifically in logistics, although such choice models offer relevant solutions. In fact, most of the studies using these methods are usually applied in the fields of (including, but not limited to) applied mathematics, biology and physics [START_REF] Lusseau | Effects of tour boats on the behavior of bottlenose dolphins: using Markov chains to model anthropogenic impacts[END_REF]- [START_REF] Randall | Analyzing Glauber dynamics by comparison of Markov chains[END_REF]. In this context, Markov chain modelling offers an interesting alternative to agent based modelling allowing efficient data processing for rich data sets production, as it will be explained in this paper

Methodology

The SIMUTRB model mixes various data sources such as pre-existing models and surveys and relies on several submodules to produce disaggregated results, which are detailed in Figure 1. The next subsections will describe the successive steps of the model, starting with the input data, which are generated by the FRETURB model but also specific surveys on urban freight transport. We then detail how the Markov process is carried out and detail the various functions that allow characterisation of the delivery tours.

Generating operations and OD matrices

The main data input of the model are O-D matrices used to initiate Markov processes. In order to produce these matrices, a generation of movements (step A) has to be processed beforehand. This specific step and further elements on FRETURB and the generation algorithm of the model are detailed in the appropriate references [START_REF] Toilier | Freight transport modelling in urban areas: The French case of the FRETURB model[END_REF], [START_REF] Routhier | FRETURB V3, a policy oriented software of modelling urban goods movement[END_REF].

In this model, the number of operations me of firm e is linked to a categorisation of the firms according to C non-empty subsets: E = {E1,E2,…Ec,…EC} such that:

e Є Ec => Me = φ(ac,wc,pc)

[1]
where a is the category of activity of the firm; w designates a class of workforce (based on the number of jobs in the firm e); p is an index that distinguishes different functions of the premises served (4 functions are distinguished: shops, warehouses, offices, head office).

Where ac, wc and pc are particular values of three variables for a firm of category c. As mentioned in the above references, the number of operations generated by zone z and by logistic category l is written as

G = M.L = gz,c,l [2]
With M a matrix indicating a number of operations by zone z and establishment category c and L a matrix of logistics structure defining the proportion of operations of logistics category l, in each establishment of category c. The logistics category l is characterised by:

the type of vehicle used (utility vehicles less than 3.5 tons, rigid trucks, semi-trailers); the management mode (own account or third party); organisation by "single leg delivery" (or "line haul") with a one stop trip (Figure 2), or by round with several trips. For the rounds, two types of operations are distinguished : o Main trips are the first and last trips of a tour o Connectors are the trips that link the establishments served by the tour to each other and between the two main trips. Indeed, in order to be paired, and for each logistics category l, these zones have to be distant from   l z, (when the trip is the first/last leg of a tour or is a line haul) or    l z, (when the trip is a connector). The pairing procedure is designed to take into account the capacity of each zone to generate operations of a category l. For example, the function    l z, indicates that an operation generated by a zone 𝑧 𝑖 and operated by a third party carrier in a LGV during a n-stop tour (defining the category l) will induce a given distance between two connectors. The target zones 𝑧 𝑗 are then selected when they are at a distance    l z, and generate an operation of the same category l. If there are no operations of category l in this zone 𝑧 𝑗 , the selection is extended to the neighbouring zones corresponding to a two sided confidence interval. The procedure is repeated for all the zones and all l categories of operations, until all operations are paired. 𝑧 𝑖 and 𝑧 𝑗 can be the same zone in the case of internal flows. The coherence is insured because the trips between a zone 𝑧 𝑖 and a zone 𝑧 𝑗 keep the same logistics characteristics (vehicle, management mode, organisation -single leg delivery or tour-and, when appropriate, size of the tour). According to a 25 class categorisation, the beginning of each l type trip which touches the zone 𝑧 𝑗 matches the movement of l type trips generated in 𝑧 𝑖 . The final O-D matrix 𝑇 𝑖𝑗 is the sum of the O-D 𝑙 𝑖𝑗 (𝑣, 𝑚, 𝑙) obtained for each vehicle class v, management mode m and classes of tour lengths c:

𝑻 𝒊𝒋 = ∑ [𝒍 𝒊𝒋 (𝒗, 𝒎, 𝒄)] (𝒗,𝒎,𝒍) [3]
This method is convergent as it is based on the consistency between the location of activities, types of operations, distances between two consecutive operations and provides us with the distribution matrix T which is used to process Markov chains.

Implementation of the Markov process

The core of the SIMTURB methodology is the implementation of a Markov process (step C). We use a discrete-time Markov chain with a finite number of states to simulate freight delivery tour chains. Markov chains respect the Markov property for which the future state of a process only depends on the present state of the same process. The chain describes a set of states noted S with (𝑋 𝑛 ) 𝑛∈𝑁 a sequence of random variables.

In this context, Markov property is therefore noted:

ℙ(𝑿 𝒏+𝟏 = 𝒙 𝒏+𝟏 | 𝑿 𝟎:𝒏 = 𝒙 𝟎:𝒏 ) = ℙ(𝑿 𝒏+𝟏 = 𝒙 𝒏+𝟏 | 𝑿 𝒏 = 𝒙 𝒏 ) [4]
A Markov chain needs the construction of a transition matrix, for which transition probabilities are noted p(x,y) the matrix being formalised as follows:

𝑷 = (𝒑(𝒙, 𝒚)) 𝒙,𝒚∈𝑺 [5]
With 𝒑(𝒙, 𝒚) ≥ 𝟎 This transition matrix is therefore a probability matrix with each row summing to 1 and being a measure of the probability S. In the case of SIMTURB, transition matrices are subsets of flows distributions matrices giving the number of operations going from zone i to zone j. Since we use a finite discrete time Markov chain to model sequences of delivery tours, it is necessary to know the actual length of tours depending on the type of tour. In our case a finite Markov chain is an efficient method to know all the states of the chain since each state is a zone delivered.

We therefore have to assign to each route a number of stops, corresponding in the Markov process to a chain length corresponding to the number of states on which we want information. This chain length cl attribution is carried out using probabilistic choice based on the type of vehicle v and management mode m.

𝒄𝒍 = ℙ(𝒗, 𝒎) [8] With ∑ 𝒑(𝒗, 𝒎) = 𝟏 𝒗,𝒎 [9]
And

∑ 𝒄𝒍 = 𝑻 𝒊𝒋 [10]
Once each route is affected with a tour length, the Markov process is carried out to build the tour to determine which zones are served. In order to keep the coherence in terms of transport behaviour, 63 transition matrices are built as subsets of the matrix T taking into account 3 main parameters: In order to balance the distribution of flows, each time a Markov chain is processed, the corresponding transition matrix is updated by subtracting the number of operations of the tour for each corresponding O-D generated in the Markov chain.

Geographical affectation -establishments

Once the zones visited by each tour are defined, the tours have to be more specifically defined with the points visited and their characteristics (e.g. the establishments). This phase (step D) corresponds to an affectation of establishments for each point of the routes simulated through the Markov chains following a two-stage sampling method. The first stage uses as primary sampling unit the establishment category e for each zone z and type of operation l in which the probability pe of being chosen is equivalent to:

𝒑𝒆 = 𝑴 𝒛,𝒍,𝒆 ∑ 𝑴 𝒛,𝒍,𝒆 𝒛,𝒍 [12]
Finally establishments are chosen within the strata defined in [START_REF] Toilier | Freight transport modelling in urban areas: The French case of the FRETURB model[END_REF] for each operation of tours depending on their category e and the number of operations of logistics category l. This specific process implemented without replacement guarantees having not the same establishment served twice in the same tour.

Each tour is therefore characterised by a set of successive points having geographical coordinates (longitude, latitude) permitting the estimation of distances thanks to Dijkstra algorithm implemented on the road network.

Following these steps, a Vehicle Routing Problem (step E) algorithm is then used to correct potential aberrations in terms of total distances travelled for the tours. This algorithm is used not to fully optimise the tours in terms of distances, but to render a faithful average distance for each tour in each category of route v,m,c. These average distances are derived from data produced from the French urban freight surveys [START_REF] Routhier | FRETURB V3, a policy oriented software of modelling urban goods movement[END_REF], [START_REF] Patier | How to improve the capture of urban goods movement data[END_REF], [START_REF] Bonnafous | French surveys of the delivery approach: from cross-section to diachronic analyses[END_REF].

Freight load affectation

Each of the delivery tour type compiles data about the initial weight loaded and their unloading profile depending on the type of activity delivered (step G). The weight loaded is a function of the number of delivery points in a tour and the type of vehicle. These elements are determined thanks to the French urban freight surveys describing 2200 routes of vehicles and 14000 stops [START_REF] Patier | How to improve the capture of urban goods movement data[END_REF]. During the delivery tour definition, three types of shape for the unloading pattern are applied (Figure 4):  the "regular" pattern that is defined by a regular unloaded weight all along the delivery tour, which equals to 50% of weight unloaded at mid-route,  the "hollowed" pattern is characterised by an unloaded weight superior to 50% of the initial weight at mid-route,  the "rounded" (or bulging) pattern is characterised by an unloaded weight inferior to 50% of the initial weight at mid-route.

These trends were determined thanks to the French surveys on urban goods movements [START_REF] Patier | How to improve the capture of urban goods movement data[END_REF], [START_REF] Bonnafous | French surveys of the delivery approach: from cross-section to diachronic analyses[END_REF]. In order to render these unloading pattern of vehicles, arithmetic functions were formalised using the main characteristics of the routes, such as the number of stops and the initial weight loaded.

The regular unloading function, corresponding to an affine function:

𝑳(𝒙) = - 𝑰 𝒏 𝒙 + 𝑰 [13]
The hollowed unloading curve can be modeled by using the following function:

𝑯(𝒙) = 𝟏 𝝆+𝛂𝒙 + 𝜷 𝒏 𝒙 [14]
The rounded unloading curve can be represented by the following function: The components being, α a coefficient determining the slope of the curve with 0<α<1 the nearer from 0 the more the curve takes a linear shape; 𝐼 the initial weight loaded in the vehicle; x the variable (referring to the stop number); n the total number of stops and In order to model loading/unloading patterns, this method needs the initial weight loaded and the number of stops. Thanks to this modelling process, it is possible to characterise a large variety of dimensions of urban freight routes, encompassing all types of vehicles, modes of organisation, modes of management, types of establishments… with a high level of spatial and temporal resolution (Figure 5, Table 1). We can consider that the level of precision of the model is mesoscopic (and not microscopic) since it is not simulated in real time with other agents and does not take into account behaviours at a vehicle level (such as gear changes, road trajectory, etc.).

The data produced in the model however offer the possibilities to work on a large variety of scenarios concerning public policies (Low Emissions Zones, Off-Hour Deliveries…) or organisational changes (Urban Consolidation Centers, fleet modifications…). It can also be used to calibrate large simulations for multi-agent models thanks to its high level of precision and vehicle based observation unit offering interoperable data input.

Confrontation to the FRETURB model

In this section, we discuss the use of Markov process to reconstitute urban freight movements and compare it to the data produced by FRETURB and data observed in the French Urban Freight Surveys. We will examine the organisational and spatial dimensions of the routes produced and verify their coherence to these data, which correspond to the state of the art. Simulations carried out to make these comparisons were made on Lyon's urban area (France) for the year 2018.

Organisational structure

As a first indicator of coherence, we compare the characteristics of the movements simulated.

The global number of movements distributed by SIMTURB, is relevant with the distribution operated by FRETURB. The number of movements in SIMTURB equals 98% of the initial FRETURB situation. However, if this first global figure is satisfying, the detail and structure shows some deviations. We indeed find some discrepancies specifically concerning the types of vehicles estimated.

First we can observe a strong convergence of SIMTURB for management modes (Table 2), reproducing the practices observed in the various urban freight surveys carried out in France [START_REF] Bonnafous | French surveys of the delivery approach: from cross-section to diachronic analyses[END_REF]. However, we find an overestimation of articulated trucks in opposition to rigid trucks (Table 3). This deviation is due to variation in the affectation of tour lengths for each type of vehicle. Fortunately these errors compensate themselves in the overall results, moreover these deviations are relatively weakly significant. The SIMTURB simulations are however much closer to the FRETURB results concerning modes of organisations (e.g. types of trips), since the attribution of trip length is entirely dependent of the 7 classes of tour lengths used in the model and guarantee limited deviations. This organisational structure has to be put into perspective with the distances generated by the model, which is a complementary indicator of coherence for the model. If we compare the results of the model to the data observed in freight surveys and generated through FRETURB, we observe globally similar trends.

Spatial structure

The previous sub-section showed strong organisational convergence between FRETURB and SIMTURB. We however noted hints of deviation concerning the spatial structure of flows specifically concerning the types of vehicles involved, which in turn have an influence on how the flows are distributed on a territory. We therefore start with an analysis of the distances generated by the model.

We observe that SIMTURB generates longer routes in some cases, specifically for vans and tours with a high number of points to deliver. This deviation is due to slightly longer inter-stop distances than usually observed which consequently has an influence on the total length of the route (Figure 7, Figure 6). These trends are identical for management modes and types of vehicles (although they tend to concentrate most on vans). The overestimation of total tour length is however compensated by shorter movements for linehaul and short tours than those measured in freight surveys. This difference is explained by the distribution of flows: movements generated in the model stay in the simulation area (the urban area) without exchanging with outside zones (regional, national and international flowswhich are taken into account in surveys although in limited proportions). Although typically representing long distance transport at regional, national or international levels, line haul trips in SIMTURB are limited to the area of study, consequently being shorter. We can effectively see that movements in shorter tours and line haul generate less distance than those observed in FRETURB and freight surveys.

Moreover short tours and line haul represent the most common types of movements (Figure 8). Since more than 50% of movements in the model are affected to line-haul and tours shorter than 10 movements, these deviations compensate themselves in the overall results. The comparison of the Origin-Destinations flows compared to the original O-D matrices generated by FRETURB is not entirely correspondent. In this section we analyse the deviation of SIMTURB when we use the Markov chains as they are originally produced (but with a strong deviation to urban freight surveys in terms of tour length) and routes optimised by a VRP algorithm. In order to measure the deviation of SIMTURB we carry out an analysis of the difference of flows according to the distances between each zone (Figure 9).

We can notice strong deviations for short O-D distance. With the original Markov process these deviations are however smaller than the VRP optimised Markov process. The optimised tours tend to overestimate flows (relatively to the FRETURB model) for O-D couples distant from 5 to 12 kilometers.

On the opposite, the model tends to slightly underestimate flows for more remote O-D couples (between 12 kilometers and more). On the overall, the errors are compensated and the models converge for longer O-D distances.

These deviations can be explained by two reasons. The first is the difference between FRETURB built O-D matrices, which are distributed on a continuous basis, and the Markov chains, which are based on a discrete and finite process. In this case, O-Ds with very low flows equal to very low probabilities in the Markov process, which tends to build chains on O-Ds with larger flows. On the contrary, O-Ds with less flows (which incidentally are usually more distant) are less likely to be drawn in the process.

The second reason is the use of a VRP, which necessarily modifies OD flows. The VRP emphasises deviations between the unoptimised tours produced through the Markov process and optimised tours: the algorithm keeps the number of points for each tour and the zones served, only rearranging the order in which the deliveries (points) are done in order to keep distances travelled consistent with observed tours in freight surveys. The number of movements therefore stays unchanged between unoptimised and optimised procedures, but the distances travelled and OD structures can evolve drastically (more than 300,000 daily kilometres travelled). We can however see that although the flows generated by non-optimised Markov processes have a better convergence with FRETURB flow structure, this process estimates longer distances, showing that the VRP is relevant, reducing the deviations in terms of distances generated, which is more relevant for evaluating pollutant emissions and congestion.

Conclusion

This paper presented an original method to simulate urban freight delivery tours using Markov chains at a mesoscopic level. After describing the general method, we analysed the deviations of the model to pre-existing data and models used for freight transport diagnosis and simulation (FRETURB and Urban Freight Surveys).

We noted that the overall organisational structure and volumes estimated by SIMTURB corresponded to FRETURB simulations and data observed in Urban Freight Surveys. However, the spatial structure of flows estimated by the model could show some deviations because of the transition between a continuous and discrete process and the use of VRP algorithm.

Although deviations to FRETURB are in some case rather significant, these differences do not necessarily mean that the model itself is not efficient in realistically reproducing urban freight tours. Therefore, the next steps to validate the model will be to confront the SIMTURB modelled data with samples of routes operated by freight carriers.

Thanks to its ability to account for the structures at stake in urban freight transport, the use of SIMTURB allows to work on a wide variety of scenarios useful for both local authorities and private operators. For instance, it can be used to examine the impact of LEZ, off-hour regulations, use of UCCs or the implementation of new types of vehicles [START_REF] Holguin-Veras | Economic impacts of staffed and unassisted off-hour deliveries in New York City[END_REF]- [START_REF] Marcucci | The potential demand for a urban freight consolidation centre[END_REF].

The model has indeed been used in various studies concerning the use of new vehicles [START_REF] Augros | Traffic simulation for optimizing an urban delivery truck -The CITYMOVE project[END_REF], Low Emission Zones, and other researches studying the impacts of UCCs and new urban distribution technologies (studies produced under non-disclosure agreement and not freely available).

For now, the model only takes into account B2B flows, but some modifications in the modelling process, coupled with behavioural data collected on household deliveries [START_REF] Gardrat | The impact of new practices for supplying households in urban goods movements: method and first results. An application for Lyon[END_REF], [START_REF] Gardrat | Méthodologie d'enquête: le découplage de l'achat et de la récupération des marchandises par les ménages[END_REF] will allow the integration of such flows. This specific calibration will however call for new typologies of movements and freight generation points (not only establishments, but also households and their characteristics).

Other enhancements will include commodity modelling based on establishments' generation.

The current loading functions are based on the initial load at the starting establishment, which is a rough estimate of the demand for goods. The model does not in fact produce data on the types of goods, packaging and the receiving establishments themselves do not generate the demand for goods. Incidentally, the model is limited to delivery tour generation and does not take into account pick-up tours or mixed organisations (alternating deliveries and pick-ups in the same tour). A modification of the modelling framework should allow taking into account pick-ups in conjunction with the goods generated by the points served by tours. This work however needs the modification of O-D matrix generation with a functional perspective on deliveries and pick-ups.

Figure 1 :

 1 Figure 1 : SIMTURB modelling framework

Figure 2 :

 2 Figure 2 : modes of organisations and types of routes taken into account in SIMTURBFrom these volumes of operations characterised by a category l, the distribution step carries on thanks to a distance matrix calculated between each pair of zones. It is therefore possible to determine the zones likely to be linked by freight trips. Based on the number of operations of the logistic category l, obtained for each vehicle class v, management mode m and classes of tour lengths c in zone z, the distribution follows (step B) an attribution process based on:

 3

 3 types vehicles (vans, rigid trucks, semi-trailers)  3 modes of management (own account as a shipper, own account as a consignee or third party)  7 classes of tour lengths ranging from line haul routes, with a one stop trip, or by tour with a minimum 3 operations and more This gives for each combination a transition matrix P, in which probabilities is a ratio of the number of operations for each O-D on the total of operations of a v,m,c subset v being the type of vehicle, m the management mode and c the class of length of the tour. 𝑷 = [ 𝒍 𝟏𝟏 (𝒗,𝒎,𝒄) ∑[𝒍 𝒊𝒋 (𝒗,𝒎,𝒄)] ⋯ 𝒍 𝟏𝒋 (𝒗,𝒎,𝒄) ∑[𝒍 𝒊𝒋 (𝒗,𝒎,𝒄)] ⋮ ⋱ ⋮ 𝒍 𝒊𝟏 (𝒗,𝒎,𝒄) ∑[𝒍 𝒊𝒋 (𝒗,𝒎,𝒄)] ⋯ 𝒍 𝒊𝒋 (𝒗,𝒎,𝒄) ∑[𝒍 𝒊𝒋 (𝒗,𝒎,𝒄)]] = [ 𝒑 𝟏,𝟏 … 𝒑 𝟏,𝒏 ⋮ ⋱ ⋮ 𝒑 𝒏,𝟏 … 𝒑 𝒏,𝒏 ]

Figure 3 :

 3 Figure 3 : Initial weight loaded (goodsin kg) according to the number of deliveries in a tour and type of vehicle

Figure 4 :

 4 Figure 4: Unloading patterns (example for 20 stops)

Figure 5 :

 5 Figure 5 : examples of routes from SIMTURB simulation results in Lyon's urban area Vehicle Management mode Nb points Length (km) Total traffic time (h) Total stopping time (h) Total load (kg) Distance between 2 stops (km)

Figure 6 :Figure 7 :

 67 Figure 6 : total distance of a tour according to the number of points for types of vehicles, envelopes for survey values

Figure 8 :

 8 Figure 8 : histogram of tour lengths (line haul -2 stops-are evicted for better readability)

Figure 9 :

 9 Figure 9: Deviation to Freturb in O-D flow structure according to O-D distances

Table 1 : examples of routes from SIMTURB simulation results in Lyon's urban area

 1

	Total tour
	duration (h)

Table 2 : Freturb -Simturb comparison for management mode structure

 2

	Vehicles	Freturb	Simturb	Difference Deviation
		flows	flows	(%)
	Vans	146816	146153	-663	-0.5%
	Articulated	29045	30906	1861	6.4%
	Rigid trucks	44864	43187	-1677	-3.7%

Table 3 : Freturb -Simturb comparison for vehicle structure

 3

	Type	of	Freturb	Simturb	Difference Deviation
	operation		flows	flows	(%)
	Connectors		117710	117536	-174	-0.1%
	Main trips		23222	23240	18	0.1%
	Line Haul		79793	79470	-323	-0.4%

Table 4 : Freturb -Simturb comparison for trip type structure

 4

Table 5 : Comparison of distances generated with Freturb -Simturb

 5