Kossigan Bernard Dedey
email: kossigan.dedey@inrae.fr

David Grenier
email: david.grenier@inrae.fr

Tiphaine Lucas
email: tiphaine.lucas@inrae.fr

Mathematical modelling of uniaxial extension of a heterogeneous gas cell wall in bread dough: Stress fields and stress concentration analysis relating to the proving and baking steps

Keywords: viscoelasticity, visco-hyperelasticity, rupture

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

Introduction

Bread dough can be viewed as a dispersion of gas cells in a continuous hydrated gluten-starch matrix [START_REF] Gan | The microstructure and gas retention of bread dough[END_REF]. The crumb structure of baked breads (specific volume, gas fraction and texture) mostly depends on dough preparation before baking i.e. mixing and proving [START_REF] Dobraszczyk | Dough rheology and breadmaking[END_REF][START_REF] Eliasson | Cereals in breadmaking: a molecular colloidal approach[END_REF] but also, to a great extent, on the way the walls that separate the gas cells (GCWs) rupture during baking [START_REF] Dobraszczyk | Dough rheology and breadmaking[END_REF][START_REF] Hayman | Factors controlling gas cell failure in bread dough[END_REF]. When early GCW rupture occurs i.e. below both gelatinization and protein coagulation temperatures, the crumb collapses and large cells form within it. When the opposite occurs, and only a few or none of the GCWs rupture during baking, there is shrinkage at the cooling stage that follows due to the decrease in pressure in the cells that have remained closed [START_REF] Kusunose | Role of starch granules in controlling expansion of dough during baking[END_REF]. At the end of fermentation and during baking, some of the thinnest GCWs are reduced to approximately the size of starch granules [START_REF] Bloksma | Dough structure, dough rheology and baking quality[END_REF][START_REF] Sandstedt | The microscopic structure of bread and dough[END_REF]. These are the GCWs most likely to rupture. When such a size is reached, gluten and starch granules must be considered as interacting phases to account for their structural heterogeneities. The individual mechanical properties of phases and the interactions between them are crucial to any proper description of stress concentrations and of the areas where rupture is most likely to be initiated within GCWs.

The mechanical properties of composite materials such as GCWs in bread dough result from the collective mechanical properties of the separate constituents present in dough and from the interactions between them. In bread dough the chief constituents are starch, gluten and water.

Numerous studies have been devoted to the characterization of the individual mechanical properties of gluten [START_REF] Dreese | Dynamic rheological properties of flour, gluten, and gluten starch doughs. I. Temperature-dependent changes during heating[END_REF][START_REF] Faubion | Dynamic rheological testing of wheat flour doughs[END_REF][START_REF] Janssen | Rheological behaviour of wheat glutens at small and large deformations. Comparison of two glutens differing in bread making potential[END_REF][START_REF] Kokelaar | Strain hardening properties and extensibility of flour and gluten doughs in relation to breadmaking performance[END_REF][START_REF] Ng | Linear to nonlinear rheology of bread dough and its constituents[END_REF][START_REF] Wesołowska-Trojanowska | Rheological properties of gluten obtained from polish wheat cultivars[END_REF] and of starch granules [START_REF] Chiotelli | Effect of small and large wheat starch granules on thermomechanical behavior of starch[END_REF][START_REF] Herrera | Rheology of starch nanoparticles as influenced by particle size, concentration and temperature[END_REF]. Studies have been performed at temperatures higher than 25°C on dough or gluten but these are generally carried out using dynamic tests [START_REF] Dreese | Dynamic rheological properties of flour, gluten, and gluten starch doughs. I. Temperature-dependent changes during heating[END_REF]. The applied strain and strain rates (> 10 -2 s -1 ) in these studies were greater than those that occur during GCW extension in the course of baking, where strain rates are of the order of 10 -3 s -1 [START_REF] Lucas | Inflation, squeezing and collapse in wheat flour dough during baking: Effects of flour quality and oven temperature[END_REF]. The conditions for most methods are such that the link between rheology and baking performance is not a straightforward one [START_REF] Dobraszczyk | Dough rheology and breadmaking[END_REF]. During baking, half of the initial extensibility of the gluten is lost at about 65°C and a proportion of the other half by the end of baking (90°C) [START_REF] Attenburrow | Rheological properties of wheat gluten[END_REF][START_REF] Grenier | Gas cell opening in bread dough during baking[END_REF]. In the case of starch alone, studies at temperatures above 25°C have been carried out in excess water [START_REF] Carrington | Microrheology of Swollen Starch Granules, Gums and Stabilisers for the Food Industry 9[END_REF][START_REF] Desse | Individual swollen starch granules under mechanical stress: evidence for deformation and volume loss[END_REF][START_REF] Fisher | Deformation mechanics of individual swollen starch granules[END_REF]. At 25°C starch granules are non-deformable at pressures of the order of magnitude found in bread baking. When heated at 65°C in excess water, starch granules soften and become easily deformable (strain > 1.4) under shearing and compression. The behaviour of starch is strongly dependent on water content and the mechanical properties identified in excess water have little or no application to starch at water contents relevant to bread dough during baking. To the best of our knowledge, no work exists that addresses the mechanical properties of starch and gluten separately in controlled hydrothermal conditions relevant to those affecting GCWs during bread dough baking. For this reason, it is hard to obtain the mechanical parameters relevant to the full temperature range encountered during bread dough baking (25-140°C). Only the orders of magnitude of these mechanical properties can be captured. [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF] identified values for starch and gluten using Young's moduli of 90 kPa and 10 kPa respectively (with relaxation times of 100 s for starch and 10 s for gluten) for bread dough at room temperature and demonstrated that there was close agreement between numerical simulations and the experimental data. This makes it possible to rely on the orders of magnitude of the moduli selected for the present study.

Various constitutive models have been used to describe the mechanical behaviour of dough constituents and compute the stress involved. One mechanical model suited to both gluten and starch during bread dough baking is the "gel" model proposed by [START_REF] Gabriele | A weak gel model for foods[END_REF] and applied by [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF], although others have considered starch to be visco-plastic [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF]. A "gel" turns from liquid to solid and makes it possible to mimic phase transitions such as the starch gelatinization and protein denaturation that occur during dough baking. This "gel" material can be modelled using time-dependent models such as the Lodge or Maxwell models using either infinitesimal or finite strain.

Views differ on the interactions between starch and gluten. Gluten and starch are known to be chemically incompatible materials and should therefore have limited interactions and should slide along each other [START_REF] Eliasson | Cereals in breadmaking: a molecular colloidal approach[END_REF]. Nevertheless, some authors have reported a number of interactions at the gluten-starch interface in dough [START_REF] He | Factors controlling gas retention in non-heated doughs[END_REF][START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF][START_REF] Petrofsky | Rheological properties of dough made with starch and gluten from several cereal sources[END_REF][START_REF] Van Vliet | Strain hardening of dough as a requirement for gas retention[END_REF], in bread crumb and starch-based biopolymer composites [START_REF] Guessasma | Influence of interphase properties on the effective behaviour of a starch-hemp composite[END_REF]. There is hence no definitive view on the nature of the interactions between gluten and starch in bread dough during baking. To enable the simulation of all possible interactions, [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF] used a thin elastic layer that made it possible to simulate either non-cohesive or cohesive interactions at the gluten-starch interface. Such a tool will be likewise be used and further described in this study.

There are few works on the numerical modelling of GCW extension at the scale of dough constituents. [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF] examined the average stress values in a large number of rheological tests using finite elements and periodic conditions at this scale. It is interesting to note that, unlike the literature and associated modelling for dough rheology, these authors considered dough to be a heterogeneous material composed of gluten, starch and the interface between the two. Unfortunately, stress concentrations and possible locations for rupture initiation in GCWs at strain rates relevant to bread during baking were beyond the scope of their study.

Following on from the work carried out by [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF], the overall objective in the current study was to extend the microscale approach adopted in dough rheology in terms of the types of law applicable to the mechanical behaviour of dough constituents and the starch fractions, strain-rate range, and variations in the mechanical properties of dough constituents associated with the baking process. Two modelling approaches for the mechanical behaviour of bread dough were compared: linear viscoelasticity (infinitesimal strain) and visco-hyperelasticity (finite strain). This work will subsequently be completed by an analysis of the sensitivity of the computed stress to some of the model's input parameters.

In the current study, the process of experimental verification has begun and has been taken as far as possible given the extreme scarcity of input data to feed the model. For this purpose, the average computed stress was compared with the very small amount of experimental data available from the literature for measurements taken at low strain rates on bread dough at ambient temperature. The simulations of the GCW extension were run in two configurations where the starch is stiffer than the gluten and two where the opposite is the case, being representative of the beginning and end of baking, before and after the starch gelatinizes and the proteins denature. The most likely locations for rupture to be initiated within the GCW for these temporally defined configurations, with and without gluten-starch interactions, were analysed through reference to stress fields using the stress triaxiality concept.

Materials and methods

Geometries

Estimation of starch fraction in the GCW at the beginning of baking

At the end of proving, the average GCW thickness in wheat dough is 240 μm [START_REF] Besbes | Characterizing the cellular structure of bread crumb and crust as affected by heating rate using X-ray microtomography[END_REF][START_REF] Turbin-Orger | Kinetics of bubble growth in wheat flour dough during proofing studied by computed X-ray micro-tomography[END_REF] with an average starch fraction of around 46% v.v [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF]. However, the present work focuses on the small proportion of GCWs with a thickness of around 10-15 µm (the average thickness of the largest starch granules) as these are the most likely to rupture [START_REF] Grenier | Gas cell opening in bread dough during baking[END_REF]. To the best of our knowledge very few reports of very thin GCWs in dough are available in the literature, and the present study has processed the microscopic images provided by [START_REF] Sandstedt | The microscopic structure of bread and dough[END_REF] that were taken at the end of baking (Fig. 1). The starch fraction was identified by measuring the areas occupied by starch and gluten in the images (Image J, National Institutes of Health, USA). The starch fraction, expressed in m 2 of starch per m 2 of GCW, was found to be 33% for image c in Fig. 1 and may slightly decrease as the GCW becomes thinner, attaining lower values locally where the starch granules are beginning to move apart (see, for example, those shown inside the green box in Fig. 2.a). Very low starch fractions (such as the 8% reference value selected for the present work) are extreme occurrences and are encountered only within strings or strands in bread crumb adjacent to locations where rupture occurs [START_REF] Stokes | In situ mechanical testing of dry and hydrated breadcrumb in the environmental scanning electron microscope (ESEM)[END_REF]. The choice of such a low starch fraction for the reference simulation was also intended to facilitate computation at the practical stages of the work and to allow a high strain level to be reached without impeding the location of stress concentrations in the final stages of the work. The fraction was then increased to 11%, 16% and 28% in order to take into account the range of thinner GCWs ready to rupture and to arrive at a closer approximation of the starch fraction reported in the literature in thin GCWs at the end of proving. In the geometries, the starch fraction was increased by adjusting the gluten dimensions (length and width).

Reference simulation

For the reference simulation, the starch fraction was set to 8% v/v because this is the fraction found at the thinnest locations within thin GCWs likely to rupture (Fig. 2.a). Where a GCW is so greatly extended the longest dimension of the starch granule is already aligned within the gluten strip [START_REF] Bloksma | Dough structure, dough rheology and baking quality[END_REF]. The single large and lenticular (A-type) starch granule in cross section measured 10 µm (long half-diameter) by 5 μm (short half-diameter). It was contained within a continuous gluten strip of 100 μm in length (L) and 20 μm wide (l) (Fig. 2.b).

Comparison with published experimental results

A second 2D geometry was used to compare the computed average stress values with published data at the scale of a continuous dough (section 3.1.3). The geometry was that of a cylinder of dough of 3.75 mm in height and 7.5 mm in diameter and corresponded to that used in Ng et al. (2006) (Fig. 2.c). Axisymmetry was also considered (Fig. 2

.d).

A third and last 3D geometry (not shown) containing an A-type starch granule (a = 5 μm, b = c =10 μm) within a gluten strip of 100 μm length, 40 μm width and 20 μm thickness was used only once to check whether the results in 2D were relevant to those in 3D (see section 3.1.1.).

Governing equations

Linear viscoelasticity

The momentum balance for the quasi-static mechanical equilibrium is given by Eq. 1.

. = 0 (1)
The total stress tensor is broken down into a purely elastic part and a viscoelastic part as follows (Eq. 2).

= +

(2)

where is the tensor (Eq. 3) which describes the time-independent elastic behaviour and the viscoelastic stress tensor which depends on the strain history of the material.

= :

(3)

where is the fourth-order tensor of elasticity and is the linear strain tensor (Eq. 4).

= ( + (

where is the displacement. The viscoelastic stress tensor is derived from the constitutive equation of the 1-element generalized Maxwell equation (Eq. 5).

+ λ = λ G (5)
where λ is the relaxation time (gluten or starch), is the infinitesimal shear strain rate tensor and G the shear elastic modulus (gluten or starch).

Visco-hyperelasticity

The momentum balance equation in the case of hyperelasticity is written as Eq. 6.

∇. ( = 0 (6) where = + is the deformation gradient where is the second-order unit tensor and is the displacement. is the second Piola-Kirchoff stress tensor (Eq. 7). It is derived from the strain energy density function (W).

S =2 = !" # + $ (7)
where = is the right Cauchy-Green deformation tensor (Eq. 8). The elastic stress tensor is given by Eq. 8.

!"

# = %&" # + '(& # = det( ,-( : (8)
where %&" # and '(& # are the volumetric and isochoric part, respectively [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]. The non-linear strain . is given by Eq. 9.

= (( + + ( . ) (9)
The stress in the 1-element viscoelastic generalised Maxwell model is obtained from Eq. 10.

$ + λ $ = λ β 0(& # (10)
The dimensionless coefficient β > 0 denotes the strain energy factor and is written as Eq. 11.

β = E %'(2& / E 45 (11)
where E %'(2& and E 45 denote the Young modulus in the viscous branch (time-dependent elasticty) and purely elastic branch (time-independent elasticity) respectively in the 2-branch Maxwell model.

In this study, the Neo-Hookean model is used for the strain energy density function (Eq. 12).

W = - 7 μ( -3 (12)
where μ is the Lamé coefficient (shear modulus) and = tr( is the first invariant of the left Cauchy-Green tensor.

Initial conditions and boundary conditions

It was assumed that there was initially no stress within either the gluten or the starch.

Conclusions on how the stress is concentrated will not be affected, even if that assumption does not stand, provided that discussion is limited to stress increases and to the identification of locations where these are greatest. The comparison to given yield stress cannot, however, be guaranteed.

Reference simulation

Strain rates of the order of 10 1 to 10 2 s -1 [START_REF] Weegels | Large and fast deformations crucial for the rheology of proofing dough[END_REF] have been estimated within gas cell membranes during failure in bread dough during proving at 34°C. These strain rates are relevant to the study of crack propagation and increase in size of a hole within the GCW once a crack has been initiated. In the present study, we only focused on stress concentration before rupture. The extension of the still un-cracked GCW is hence driven by the growth of the gas cells between which it is sandwiched. At the end of proving and during baking, GCWs, before rupture, typically undergo strain rates ranging from 10 ,< to 5 × 10 ,= > ,-and Hencky strain above 1 [START_REF] Dobraszczyk | Dough rheology and breadmaking[END_REF][START_REF] Eliasson | Cereals in breadmaking: a molecular colloidal approach[END_REF][START_REF] Turbin-Orger | Growth and setting of gas bubbles in a viscoelastic matrix imaged by X-ray microtomography: the evolution of cellular structures in fermenting wheat flour dough[END_REF]. The strain rate depends on the GCW location within the loaf and decreases during baking. For the reference simulation, a typical decrease in strain rate at the core of the dough [START_REF] Lucas | Inflation, squeezing and collapse in wheat flour dough during baking: Effects of flour quality and oven temperature[END_REF] from 3 × 10 ,= > ,-to 1.15 × 10 ,= > ,-was used. The displacement u Ax, ± E 7 F was applied accordingly at the upper and lower boundaries (Eq. 13).

u Ax, ±

E 7 F = E 7 ε (t = 0 t (13)
where ε (t = 0 is the initial strain rate (3 × 10 ,= s ,-and t is the current time. The left and right-hand sides of the GCW are free to move (Fig. 2.b).

Comparison with published experimental results

In order to compare the simulations with published experimental results for the cylindrical geometry [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF], a constant strain rate was applied. Upper and lower boundary displacements were accordingly applied to match the conditions used in the experiment (Eq.

14)

. The right-hand side of the cylinder was free to move and there was no displacement on the axis of symmetry (Fig. 2.d).

u Ax, ± E 7 F = E 7 (e I (JKL J -1 (14)

The gluten-starch interface: the Thin Elastic Layer (TEL)

Gluten-starch interaction at the gluten-starch interface was modelled using a Thin Elastic Layer (TEL) boundary condition. The thickness of a very thin layer can easily be defined without taking specific areas into account, the thickness M of the water phase at the glutenstarch interface being of the order of hundreds of nanometers at most. The TEL decouples the displacements on the two sides of the boundary. The interface is characterized by the two elastic constants k O and k J , in normal (Eq. 15) and tangential (Eq. 16) directions respectively.

k O = P QRS (-,T QRS !(-UT QRS (-,7T QRS (15) k J = V QRS ! (16)
where ν 'OJ is Poisson's ratio, E 'OJ is Young's modulus and G 'OJ is the shear modulus of the material at the interface. The interaction between gluten and starch was cohesive when E 'OJ was set to 10 6 Pa and non-cohesive when it was set to 10 -6 Pa.

Material properties

Reference properties

Both gluten and starch were considered to be viscoelastic. As neither gluten nor starch was assumed to exhibit much time-independent elasticity, the time-independent Young's modulus was set to 10 Pa. Little information on this time-independent elasticity is to be found in the literature, where discussion of the residual stress shown in the experimental data is limited. In the present work, the mechanical behaviours of both gluten and starch were mostly controlled by their time-dependent elasticity as it relaxed. At the beginning of baking, gluten is highly deformable while starch granules are non-deformable at the gas pressure values relevant to baking. As the temperature increases above the range 70-80°C, gluten denaturation and subsequent cross-linking are accompanied by an increase in gluten Young's modulus [START_REF] Grenier | Gas cell opening in bread dough during baking[END_REF] (Fig. 3). By contrast, the uptake of water by starch granules and the associated phase transitions of starch upon heating are accompanied by a decrease in the rigidity of the starch granules. In excess water, starch granules lose their integrity even under gentle shearing, the temperature at which this transition occurs varying with the starch type.

Mechanical testing of individual starch granules immediately after heating in excess water has shown that potato starch granules can become highly deformable at pressures of the order of hundreds of Pa [START_REF] Carrington | Microrheology of Swollen Starch Granules, Gums and Stabilisers for the Food Industry 9[END_REF][START_REF] Desse | Individual swollen starch granules under mechanical stress: evidence for deformation and volume loss[END_REF][START_REF] Fisher | Deformation mechanics of individual swollen starch granules[END_REF]. At the levels of hydration relevant to bread dough, the extent of granule softening remains uncertain.

Microscopic observations of GCWs revealed that the granules appeared to have flattened significantly in wheat bread crumb [START_REF] Sandstedt | The microscopic structure of bread and dough[END_REF] and to an extreme extent in bread prepared with low-amylose starch which is reputed to soften/disrupt wheat starch at an earlier point than usual during heating [START_REF] Kusunose | Role of starch granules in controlling expansion of dough during baking[END_REF]. Where baking occurs at atmospheric pressure, the release of carbon dioxide in gas cells becomes limited in the range of temperatures (above 60-70°C) at which granule softening is likely to occur. This makes it highly improbable that GCWs will reach any great degree of extension in these conditions and there is therefore no sense in studying high levels of extension for such cases. It is relevant to do so, though, for innovative baking processes such as partial vacuum baking [START_REF] Grenier | Enhanced aeration of part-baked bread using a novel combination of baking and partial vacuum[END_REF][START_REF] Lucas | Facility and method for treating food products, such as cellular products, in particular for expanding said products[END_REF][START_REF] Rondeau-Mouro | Characterization of gluten-free bread crumb baked at atmospheric and reduced pressures using TD-NMR[END_REF][START_REF] Şimşek | Evaluation of partial-vacuum baking for gluten-free bread: Effects on quality attributes and storage properties[END_REF] where extension is forced at mid-baking by the decrease in pressure in the oven's atmosphere. This decrease in pressure is accompanied by water ebullition and gas extension and, consequently, GCW extension is enhanced in these conditions.

Two sets of Young's moduli were therefore considered for the computations. The first relates to early baking (before 50-60°C is reached) when the starch (100 kPa) is more rigid than the gluten (10 kPa) [START_REF] Attenburrow | Rheological properties of wheat gluten[END_REF][START_REF] Dreese | Dynamic rheological properties of flour, gluten, and gluten starch doughs. I. Temperature-dependent changes during heating[END_REF][START_REF] Khatkar | Dynamic rheology of wheat flour dough. II. Assessment of dough strength and bread-making quality[END_REF][START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF], and the second reflects the conditions of advanced baking (beyond 70°C) when gluten has low deformability (100 kPa) [START_REF] Attenburrow | Rheological properties of wheat gluten[END_REF]) and starch deformability is high (10 kPa) [START_REF] Fisher | Deformation mechanics of individual swollen starch granules[END_REF] (Fig. 3). The orders of magnitude for the first set derive from the work of [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF] and are simply inverted to form the second set because no data for wheat starch at moderate levels of hydration and high temperatures is available in the literature. For the reference simulation, the same relaxation time (10 s) was used for both gluten and starch. These parameters are assumed to remain constant throughout the extension (Fig. 3). All parameters are gathered in Table 1.

Comparison with published experimental results

Since the dough was treated as uniform in the experiment used for validation, the mechanical properties were also considered to be uniform throughout the dough cylinder for that specific case. The average Young's modulus of the dough (100 kPa) was chosen. Two relaxation times (10 and 100 s) were also used in order to evaluate the sensitivity of the average stress to strain rate.

Digital modelling

Finite element computations were carried out using the COMSOL Multiphysics® v. 5.4 (COMSOL AB, Stockholm, Sweden) MUMPS solver . A mesh convergence test was run to find the best balance between simulation time/use of computer resources and the stability of the model. The average run time varied from 10 minutes to 1 day on an Xeon (R)

W-2155 CPU Intel processor at 3.31 GHz, with 256 Go RAM.

Analysis of stress fields and the triaxiality factor

For stress analysis, the von Mises criterion was applied (Eq. 17). This equivalent tensile stress (known as von Mises stress, X YZ ) is that most commonly used for the analysis of yielding in materials science.

X YZ = [ (\ ]] ,\ ^^ ^U(\ ^^U\ __ ^U(\ __ U\ ]] ^U`(\ ]^^U \ ^_^U \ _] ^ 7 (17)
The triaxiality factor (a. b., Eq. 18) was also used to identify the regions of uniaxial extension and shear. When the triaxiality factor is around 0.33 the material is mostly uniaxially extended and when it is equal to zero the material is mostly sheared. In the stress, the numbers refer to the direction of the space. The first number is the direction normal to the surface upon which the stress is applied and the second is the direction of the component of the stress.

a. b. = \ cd \ ef (18)
where

X 4g = - = (X) = - = (X --+ X 77 + X == )
3 Results and discussion

Model analysis and experimental validation

Average stress values computed using linear-viscoelasticity in 2D and 3D were compared in order to check that the agreement between them was sufficient for the entire study to be conducted using 2D geometries (section 3.1.1). Viscoelasticity and visco-hyperelasticity approaches were then compared to assess where visco-hyperelasticity was more suited to the task than linear viscoelasticity (section 3.1.2). Last, the average stress values computed using hyperelasticity were compared to previously published experimental data [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF] (section 3.1.3).

Comparison between 2D and 3D geometries

Fig. 4 shows the average values for von Mises stress vs Hencky strain in 2D and 3D geometries using the reference simulation for the two selected cases (starch granule more rigid than gluten/gluten more rigid than starch granule) (Fig. 4). Where the starch granule was more rigid than the gluten and up to a Hencky strain of 0.67, the difference when using linear visco-elasticity in 2D and 3D did not exceed 2%. Where the gluten was more rigid than the starch granule, this difference reached the order of 20% at most. These discrepancies are in line with those reported by [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF]. 2D rather than 3D geometry was therefore used to reduce computation times without too greatly distorting the results.

Comparison between linear viscoelasticity and visco-hyperelasticity

The average values for von Mises stress vs Hencky strain up to 0.67 are shown in Fig. 5 for both linear viscoelasticity and visco-hyperelasticity. There was good agreement between simulations up to a strain of 0.08. Thereafter, the average stress values diverged. It is notable in Fig. 5 that, when E (Jh 2i /E j"kJ!O =0.1, the viscoelasticity model barely registers the strain softening observed experimentally by [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF] at low strain rates. For this reason, we selected the visco-hyperelasticity model for the next steps of the study.

Comparison with published experimental results

The average von Mises stress values computed using the axisymmetric geometry were compared to published experimental results obtained at ambient temperature from dough [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF]. The objective here was solely to validate the structure of the model. It is impossible to achieve greater quantitative validation because the mechanical properties used in the model are unlikely to have been exactly those of the dough used in the experiment. Fig. 6 shows both simulated and experimental average values for stress vs Hencky strain at different strain rates (0.3 s -1 and 0.003 s -1 ). At the high strain rate (0.3 s -1 ), no good quantitative agreement could be found between the simulations in the present study and the experiments reported in the literature (Fig. 6.a). This is probably due to an overestimate of the Young's modulus. There is unfortunately no guarantee that the elasticity of the dough used by [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF] in their experiment was not lower than that used in the simulation. However, similarities in strain-hardening were observed between the simulation and the experiment performed by [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF] and the stress exhibited almost no sensitivity to relaxation time (Fig. 6.a). Elasticity was the main driver for the stress in this high strain rate configuration. At the low strain rate (0.003 s -1 ), when the model was less sensitive to elasticity and far more so to viscosity, a quite good quantitative and qualitative agreement was found between the simulations and the experiments (Fig. 6.b). As expected, the average von Mises stress values within the GCW were strongly dependent on the relaxation time at this low strain rate (Fig.

6.b)

. This indicates that the estimation of the relaxation time for dough (or dough constituents) has to be quite accurate when addressing GCW extension during proving or bread baking, since it is the viscous aspect of the dough that mostly drives the stress under such low strain rate conditions. Note that this statement is only true if the assumption of very low time-independent elasticity was true. At low strain rates, precisely the opposite conclusions would be drawn if there were a degree of time-independent elasticity in the dough. In such a case, only time-independent elasticity would be relevant because all possible time-dependent elasticity relaxation would have already occurred.

Average stress values

Reference starch fraction

Fig. 7 plots the average values of von Mises stress as a function of the Hencky strain for starch to gluten moduli ratios of 0.1 and 10. The starch-gluten moduli ratios were found to have significant impact on the average von Mises stress. Average stress values where starch was less rigid than gluten were tenfold higher than those where starch was more rigid than gluten. This was largely due to the high proportion of gluten, which accounted for 92% of the total volume.

The nature of the interactions between gluten and starch was found to have no effect on the average von Mises stress within the GCW for the reference starch fraction (8%). This does not fit with the reported results of [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF], who found that the cohesive/noncohesive nature of the gluten-starch interface affected average stress values in every case tested in their study. The difference is explained by their choice of starch fraction, which averaged around 46 % in bulk dough without gas cells [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF][START_REF] Tanner | Bread dough rheology: an improved damage function model[END_REF]. The following section will address the different starch fractions and will extend the comparison with the data reported by [START_REF] Mohammed | Mechanical characterization and micromechanical modeling of bread dough[END_REF].

Effect of starch fraction

Fig. 8 shows the average von Mises stress as a function of Hencky strain for the reference (blue color) starch fractions of 11, 16 and 28% and compares the cohesive and non-cohesive hypotheses concerning the interface between starch and gluten; only the case where the starch is more rigid than the gluten is considered here. The nature of the interaction increasingly affected average von Mises stress values as the starch fraction increased (Fig. 8). From this, it can be concluded that the nature of the gluten-starch interaction affected the average stress within the GCW as soon as the starch fraction exceeded a threshold located between 8% and 11%. As expected, the decrease in average stress values that accompanied the increase in starch fraction was greater when the interface was non-cohesive. Slightly greater average stress values were found when a cohesive interface was involved. Indeed, within the area of gluten closest to the starch granule, the increase in starch fraction was accompanied by a greater increase in stress values because the gluten in this area was under greater strain than it was when the interface was non-cohesive.

Stress fields and triaxiality: stress concentrations in early and late baking

The objective in this section is to refine the foregoing analysis by considering local stresses rather than average values. The intention is to learn more about the most likely locations for GCW rupture in the different configurations. From this point on, we will adopt the working assumption that the locations with the highest concentrations of stress are where the rupture of the material (the GCW in this instance) is most likely to be initiated. In other words, we assume that rupture occurs when the yield stress value is exceeded. In a further step, this analysis will make use of these findings to predict the phenomena most likely to be present in microscopic observations of the GCW in baked bread crumb. It should be remembered that the range of moduli ratios that have been tested in the course of the present work are also of significance to the baking process and to gradual changes in molecular conformation upon heating, as detailed in section 2.4.1.

Fig. 9 and Fig. 10 show von Mises stress fields for starch fractions of 8% (reference fraction) and 28% at 0.67 and 0.4 Hencky strains respectively, for different starch-gluten modulus ratios (0.1 and 10) and for cohesive/non-cohesive interactions. We next provide the refined analysis referred to above for the reference fraction. The locations of stress concentrations in the 28% starch fraction were found to be identical to those in the reference fraction. Only the magnitude of stress and strain (0.41) before modelled failure differed.

Early baking, l mnop q > l rs ntu

Fig. 9 a depicts results when the starch granule is more rigid than the gluten for the two gluten-starch interface types: cohesive (left) and non-cohesive (right).

Where it was assumed that the interaction at the gluten-starch interface is cohesive, the starch granule reinforced the gluten film and the stress spread through the entirety of the gluten and starch by means of the cohesive interaction (Fig. 9.a, on the left-hand side). The rupture was most probably initiated at two points of extreme stress within the gluten film; one lying a few micrometres away from the starch granule in the direction of the extension within the gluten and another immediately adjacent to the rim of the granule (see arrows in Fig. 9.a, left). The maximum stresses were 1.8 × 10 3 Pa and 2.08 × 10 4 Pa for the cohesive and non-cohesive gluten-starch interfaces respectively and were thus separated by an order of magnitude. It is worth remembering that yield stress magnitudes in dough for low strain rates (0.1 s -1 ) ranged from 0.5 to 1.1 × 10 3 Pa depending on the quality of the gluten network [START_REF] Attenburrow | The fracture behaviour of starch and gluten in the glassy state[END_REF][START_REF] Chin | Dough aeration and rheology: Part 2. Effects of flour type, mixing speed and total work input on aeration and rheology of bread dough[END_REF][START_REF] Dunnewind | The kieffer dough and gluten extensibility rig-an experimental evaluation[END_REF]. Such a thin GCW containing a still-rigid starch particle is likely to have ruptured at a lower Hencky strain than 0.67. This would be all the more probable if a certain amount of stress had been previously stored in the GCW at the beginning of the extension process, as can be expected to occur in proven dough at the beginning of baking.

Fig. 9 also shows the triaxiality factor fields that provide information on those areas where unidirectional extension (Fig. 9.c) and shear (Fig. 9.e) occur (cohesive hypothesis shown on left). In those cases where the gluten-starch interface was cohesive, the uniaxial extension within the gluten was strongest close to the rims of the starch granules and spread along the sides of the granules (Fig. 9.c, left). Shear was also found within the gluten in the upper and lower regions close to the gluten-starch interface. This shearing probably increases the likelihood that the interface will rupture at these locations. If the cohesive assumption is valid, these simulations tell us that microscopic observations should reveal gluten in close contact with the starch at the mid-point of the granule's length and gluten shreds near the edges of the granule.

By contrast, non-cohesive interactions at the gluten-starch interface caused a concentration of stress that was largely confined to the gluten lying in the direction of the extension close to the very edge of the starch granule (Fig. 9.a, right). This was where the greatest uniaxial extension of the gluten film occurred (Fig. 9.c, right). It was also the point at which the GCW was most likely to rupture. The starch granule did not reinforce the gluten as it had done in the previous case and stress within the granule was much lower (see Fig. 9.b vs 9.a). In this hypothesis, once initiated, rupture will spread within the gluten and is then expected to slide around the starch granule since there is no interaction at the gluten-starch interface. Stress values were also high in the gluten film that was close to the lateral edge of the starch granule (Fig 9.a,right) and this constitutes a second highly plausible location for rupture to occur (see arrow in Fig. 9). If the non-cohesive hypothesis is valid, these simulations tell us that, taking into account the uncertainty as to the degree of slippage, the starch granules located closest to the point of rupture might be partially covered with gluten but, unlike those in the cohesive model, they should be totally free of gluten shreds at their edges.

3.3.2 Late baking, l mnop q < l rs ntu Fig. 9.b shows von Mises stress fields at 0.67 Hencky strain when the starch granule was less rigid than the gluten for cohesive (left) and non-cohesive (right) gluten-starch interfaces.

When it was sufficiently soft, the starch granule became elongated in the direction of the extension and the stress became concentrated along the lateral surfaces of the starch granule because of the extreme extension of the gluten along both sides of the granule (Fig. 9.b and d, green arrows). The rupture of the gluten was most likely to occur along the surface of the granule where both strain and stress were at their greatest. Here again, with a non-cohesive gluten-starch interface, the surface of the starch granule should be free of gluten following rupture. With a cohesive gluten-starch interface, some gluten should remain on the lateral surfaces of the starch granule. In the direction of extension, the stress decreased with distance from the starch granules. It is possible that he gluten that had previously been strained in these areas when it had been less rigid than the starch (replicating the early stage of baking, see section 3.3.1) had now been attenuated by the elongation of the starch granule and the stress redistributed between the gluten and the starch. Such a change is likely to increase the extensive capacity of the GCW but only if starch granule softening occurs before the gluten has become too stiff. We should note that these two events are not well documented in the literature. Note also that if granule softening were to occur very early in the baking process it is possible that the GCWs might rupture far less frequently, leading to an incomplete opening of the porous structure of the bread's crumb. This would cause the baked dough to shrink during the cooling step. Again, this transition cannot be analyzed further without information on the evolution of the starch modulus during the heating process.

This analysis has also revealed that the main difference found between cohesive and noncohesive interactions at the gluten-starch interface was an element of shear on the side of the granule where cohesive interactions were involved but none where interactions were noncohesive (Fig. 9.f, left vs right). This explains why the stress spread throughout the starch and gluten and why the maximum stress remained lower than when the gluten-starch interactions were non-cohesive. In this case, the triaxiality factor showed that uniaxial extension and shear mostly decreased with distance from the granule, following the direction of the decrease in stress (Fig. 9.d and f).

About the incompressibilty of the starch granule

The conclusions drawn up to this point rely on the hypothesis of incompressibility of the starch granule. Unfortunately we found no literature dealing with Poisson's ratio of starch granule to sustain the hypothesis of incompressibily. In order to identify how far the non validity of the assumption of starch incompressibility could affect the results, simulations were run using Poisson's ratio of 0.3, which is an extreme value for polymers (see supplementary material 1). On the one hand, previous conclusions on the location where the stress concentrates the most were unchanged. On the other hand, the effect of cohesion/noncohesion between gluten and starch on the mean von Mises stress decreased with decreasing

Poisson's ratio. This effect was almost totally cancelled (remaining difference of 2.5%) in the late baking configuration and was quite reduced (-20%) in the early baking configuration, where the starch granule was more rigid than gluten.

Conclusion and future research

Both linear viscoelastic and visco-hyperelastic models were assessed to determine their suitability for the computation of average stress and of stress fields during uniaxial extension within gas cell walls (GCWs) in bread dough during baking. Comuptation of viscohyperelasticity involving finite strain was found to be more suited to the task of replicating GCW extension under strain and strain rates close to those encountered during baking. The visco-hyperelasticity model made it possible to describe strain-hardening in sufficient accordance with experimental results [START_REF] Dobraszczyk | Dough rheology and breadmaking[END_REF][START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF][START_REF] Van Vliet | Strain hardening of dough as a requirement for gas retention[END_REF] for a high strain rate (ε = 3.10 ,-s ,-. Strain softening was computed at low Hencky strain for a low strain rate (ε = 3.10 ,= s ,-) in line with the work of [START_REF] Ng | Linear to non-linear rheology of wheat flour dough[END_REF]. The simulations also demonstrated that proper identification of relaxation times is important for the analysis of low strain rates such as those encountered during baking provided that not too much timeindependent elasticity is involved.

It was found that the nature of gluten-starch interactions had no significant effect on the average stress for the low starch fraction used as a reference (8%). Whether the gluten-starch interaction was cohesive or non-cohesive did, however, become significant when the starch fraction was increased. For a starch granule fraction of 28%, which was the closest to that reported in the literature for thin GCWs and which was estimated in the present study, the mechanical interaction between starch and gluten was found to have an impact on average stress values, as already evidenced for bulk dough in the literature. The nature of the interactions between starch and gluten during baking is still a matter of debate in the cereal science community. Our simulations confirmed that information on the nature of the interaction at the gluten-starch interface such as that proposed by Jekle et al. ( 2016) is important for the appropriate modelling of the stress fields within GCWs. The rupture location will depend on the nature of this gluten-starch interaction which can be determined by the presence or absence of a number of gluten shreds along the surfaces of the starch granules located closest to the hole produced by the rupture. This evidence may provide a useful indicator for the future microscopic analysis of holes in the GCWs of bread crumb.

This result also confirms that the interaction between starch granules and gluten may affect dough performance and that gluten and starch should be considered in interaction [START_REF] Gao | Influence of gluten and starch granules interactions on dough mixing properties in wheat (Triticum aestivum L.)[END_REF] rather than separately, as is commonly practiced in the literature. The experimental verification of the model remains insufficient, requiring further experimental work to identify the mechanical properties of dough within the relevant temperature range for baking carried out at low strain rates and in a controlled gaseous environment to reflect real-life conditions for bread dough during baking. Work should first be carried out to estimate the timeindependent elasticity of the protein matrix and liquid lamella at the dough-gas interface as this might strongly affect strain-hardening at low strain rates. It is almost impossible to record the low strain rates relevant to dough proving and baking using conventional rheometers. A dedicated experimental device for bubble inflation that is capable of reproducing both low and high strain rates in an atmosphere closely resembling that encountered in gas cells during the baking of real dough must therefore be developed. Such a device would enable the elimination of doubt concerning the mechanical properties of dough at high temperatures and provide data to support the conclusions drawn in the present study. The development of this device will form the next step in our investigation. Second, the possible evolution of the interactions at the gluten-starch interface between the beginning and end of baking and the changes that are even more likely to occur between the end of baking and the end of cooling require analysis. This work is necessary for the appropriate modelling of the stress fields in GCWs during the last stages of bread making to be carried out.

Fig. 1 .

 1 Fig. 1. Microscopic views of GCW cross-sections at the end of fermentation from Sandstedt (1954): a. thick GCW; b., c. and d. thin GCWs.

Fig. 3 .

 3 Fig. 3. Mechanical properties for T < 50-60°C, E /E =10 (early in baking) and T > 70-80°C, E /E = 0.1 (later in baking). The GCW spatial organisation and morphology shown in the upper part of the figure has been adapted from Grenier et al. (2021).

Fig. 4 .

 4 Fig. 4. Average values for von Mises stress vs Hencky strain in the case of cohesive interaction at the gluten-starch interface.

Fig. 5 .

 5 Fig.5. Plotted average values for von Mises stress vs Hencky strain, comparing linear viscoelasticity with visco-hyperelasticity where there is cohesive interaction at the gluten-starch interface.

  Fig. 9. von Mises stress fields within the GCW at the end of extension (0.67 Hencky Strain) for the two references E /E =10 (a, c, e) and E /E = 0.1 (b, d, f). Triaxiality factor: uniaxial extension (TF = 0.33, c and d) and shear (TF = 0, e and f). Cohesive (left-hand image) and non-cohesive (right-hand image) gluten-starch interfaces are shown for each subplot a, b, c, d, e and f.

  Fig. 7. Reference starch fraction -average values for von Mises stress vs Hencky strain computed with the proposed model for both the cohesive and non-cohesive gluten-starch interface hypotheses and for two different starch-gluten modulus ratios E Fig. 8. Starch fractions -average values for von Mises stress vs Hencky strain computed with the proposed model, for different starch fractions (11%, 16% and 28%), the last being of relevance to GCWs in bread dough at the end of fermentation. All calculations were performed for the case

				a 3,000					0.3 s -1
				400,000					E starch /E gluten = 0.1
			von Mises stress [Pa]	1,000 1,500 2,000 2,500 2,500 350,000 3,000 100,000 150,000 200,000 250,000 300,000 Von Mises stress [Pa] 1,000 1,500 2,000 Von Mises stress [Pa] 50,000 500			E starch /E gluten =0.1 Linear viscoelasticity Visco-hyperelasticity Model Ng, McKinley (2006) t = 100 s t = 10 s Cohesion Non-cohesion Strong wheat flour dough E starch /E gluten = 10
					0.0	0.2	0.4	0.6	0.8	1.0	1.2
				500	0.0 0	0.1	0.2	E starch /E gluten =10 Hencky strain, Ɛ 0.3 0.4 0.5 H	0.6	0.7
										Hencky strain, Ԑ
				0
				0.0	0.1	0.2			0.3 0.003 s -1 0.4	0.5	0.6	0.7
					80,000			Hencky strain, Ԑ H
	E	/E	= 10.	Von Mises stress [Pa]	20,000 40,000 60,000		t = 100s Ng, McKinley (2006) Model	/E	= 0.1 and E	/E	= 10.
											t = 10s
					0
					0.0	0.2	0.4	0.6	0.8	1.0	1.2
									Hencky strain, Ԑ
	Fig. 6. Average values for von Mises stress vs. Hencky strain, comparing the proposed model with
	experimental data obtained at two different strain rates (0.3 s -1 and 0.003 s -1 ) (Ng et al., 2006).

H b H

Table 1 .

 1 Nomenclature and input parameters

	Unit	Value	Ref

Acknowledgments

The first author, Kossigan Bernard DEDEY, is the holder of a grant from Université Rennes 1 and INRAE.


Manon Bourdeau

MSc Nicolas Guibert

Dr Emmanuel Fort
email: emmanuel.fort@univ-lyon1.fr

Sebastien Boulogne

Dr Emmanuel Lagarde

PhD Barbara Charbotel

MSc Sébastien Boulogne

Medicine consumptions and occupational road risk

Keywords: Road safety, responsibility analysis, road traffic crash, professional drivers, medicines

Introduction

Several studies in France [START_REF] Hours | Diseases, consumption of medicines and responsibility for a road crash: a case-control study[END_REF] and around the world have been carried out about the link between the consumption of several medicines and the risk of road traffic crash [START_REF] Elvik | Risk of road accident associated with the use of drugs: A systematic review and meta-analysis of evidence from epidemiological studies[END_REF]. Some drugs are identified as strongly associated to an increased risk of road crash, especially benzodiazepines which represent the most studied class in the literature [START_REF] Orriols | Benzodiazepinelike hypnotics and the associated risk of road traffic accidents[END_REF][START_REF] Orriols | Trends in benzodiazepine anxiolytics and z-hypnotics use among French drivers involved in road traffic crashes from 2005 to 2015: a responsibility case-control study[END_REF][START_REF] Nevriana | New, Occasional, and Frequent Use of Zolpidem or Zopiclone (Alone and in Combination) and the Risk of Injurious Road Traffic Crashes in Older Adult Drivers: A Population-Based Case-Control and Case-Crossover Study[END_REF][START_REF] Engeland | Risk of road traffic accidents associated with the prescription of drugs: a registry-based cohort study[END_REF][START_REF] Lai | Long-term use of zolpidem increases the risk of major injury: a population-based cohort study[END_REF][START_REF] Drummer | Odds of culpability associated with use of impairing drugs in injured drivers in Victoria, Australia[END_REF][START_REF] Gjerde | Alcohol, psychoactive drugs and fatal road traffic accidents in Norway: A case-control study[END_REF].

Other drugs have been less studied, with sometimes inconsistent or contradictory results. This is the case, for example, with anti-epileptic drugs [START_REF] Orriols | Prescription of antiepileptics and the risk of road traffic crash[END_REF][START_REF] Sundelin | Epilepsy, antiepileptic drugs, and serious transport accidents: A nationwide cohort study[END_REF][START_REF] Bramness | An increased risk of road traffic accidents after prescriptions of lithium or valproate?[END_REF], analgesics [START_REF] Wickens | Driving under the influence of prescription opioids: Self-reported prevalence and association with collision risk in a large Canadian jurisdiction[END_REF][START_REF] Monárrez-Espino | Analgesics and road traffic crashes in senior drivers: An epidemiological review and explorative metaanalysis on opioids[END_REF], antidepressants [START_REF] Sims R V | Relations among Chronic Medical Conditions, Medications, and Automobile Crashes in the Elderly: A Population-based Cas-Control Study[END_REF][START_REF] Meuleners | Psychoactive medications and crash involvement requiring hospitalization for older drivers: a population-based study[END_REF][START_REF] Orriols | Risk of injurious road traffic crash after prescription of antidepressants[END_REF][START_REF] Chang | Psychotropic drugs and risk of motor vehicle accidents: a population-based case-control study[END_REF][START_REF] Leveille | Psychoactive medications and injurious motor vehicle collisions involving older drivers[END_REF], muscle relaxants [START_REF] Bramness | The risk of traffic accidents after prescriptions of carisoprodol[END_REF] and antihistamines [START_REF] Ray | Psychoactive drugs and the risk of injurious motor vehicle crashes in elderly drivers[END_REF][START_REF] Howard | Sleepiness, sleep-disordered breathing and accident risk factors in commercial vehicle drivers[END_REF][START_REF] Verster | Antihistamines and driving ability: evidence from on-the-road driving studies during normal traffic[END_REF][START_REF] Orriols | Road traffic crash risk associated with prescription of hydroxyzine and other sedating H1-antihistamines: A responsibility and case-crossover study[END_REF][START_REF] Woratanarat | Alcohol, illicit and non-illicit psychoactive drug use and road traffic injury in Thailand: a case-control study[END_REF]. In addition, there are several substances in some of these drugs classes that may have an impact on driving ability. This is the case of substances used as hypnotics or antiepileptics, which are likely to have a sedative effect; or some substances used for diabetes that can cause hypoglycemia [START_REF] Orriols | Prescription Medicines and the Risk of Road Traffic Crashes: A French Registry-Based Study[END_REF].

Studies on the effects of medicines on the risk of road traffic crash are often focused on the elderly [START_REF] Monárrez-Espino | Analgesics and road traffic crashes in senior drivers: An epidemiological review and explorative metaanalysis on opioids[END_REF][START_REF] Sims R V | Relations among Chronic Medical Conditions, Medications, and Automobile Crashes in the Elderly: A Population-based Cas-Control Study[END_REF][START_REF] Meuleners | Psychoactive medications and crash involvement requiring hospitalization for older drivers: a population-based study[END_REF][START_REF] Ray | Psychoactive drugs and the risk of injurious motor vehicle crashes in elderly drivers[END_REF][START_REF] Labat | Prevalence of psychoactive substances in truck drivers in the Nord-Pas-de-Calais region (France)[END_REF]). Studies that examined the relationship between drug use and the risk of work-related motor vehicle crashes are scarce, often inconsistent and mainly focused on truck drivers [START_REF] Labat | Prevalence of psychoactive substances in truck drivers in the Nord-Pas-de-Calais region (France)[END_REF][START_REF] Khoshakhlagh | The relationship between illnesses and medical drug consumption with the occurrence of traffic accidents among truck and bus drivers in Tehran, Iran[END_REF][START_REF] Girotto | Continued use of drugs and working conditions among truck drivers[END_REF].

However, it can be argued that drivers on work-related journeys are more likely to have a different coping behavior due to many factors [START_REF] Fort | Road accidents, an occupational risk[END_REF][START_REF] Fort | Working conditions and risk exposure of employees whose occupations require driving on public roads -Factorial analysis and classification[END_REF]. For example, they may be more exposed to stress and fatigue because they are often required to move quickly and efficiently and because of the pressure that hierarchy can exert. This could lead to different medicine consumption behaviors compared to non-professional drivers.

On the other hand, there is reason to believe that some drivers on occupational journeys may have lower consumption prevalence than drivers on private journeys. Indeed, some of them are subjects to specific regulations. Moreover, their journeys do not have the same stakes or objectives as drivers on private journeys, as driving may represent a (financial) livelihood.

To our knowledge, no study has compared the differences in consumption and road crash risk according to the type of journey, private or work-related (either commuting or mission).

The objective of this study was to determine whether the prevalence of drivers' medicines use differed between work-related journeys (commuting or mission) and non-work-related (private) for ten major classes of drugs and to describe the impact of these drugs on the risk of being responsible for a road crash according to the type of journey.

Methods

Data sources

The study is based on three French nationwide databases: police reports (PR), the national police database of injurious crashes (IC) and the national healthcare insurance database (HCI database).

The drivers were included in the study by their national ID number (NID, also called social security number) that was extracted from police reports (PRs) by an automatic procedure.

Next, drivers were linked to prescription medicines reimbursement data around the date of the crash by their NID. Then, PRs were linked to the National police database of injurious road traffic crashes (IC) records by a probabilistic linkage method [START_REF] Orriols | Prescription Medicines and the Risk of Road Traffic Crashes: A French Registry-Based Study[END_REF].

Police reports (PRs)

In France, when a road crash with injuries occurs, police forces are called and have to file a report. This report contains general information about the crash. It is then scanned and stored as an image file by TransPv. Almost 70,000 reports are completed each year. The NID is sometimes recorded in the PR. It has been estimated, through manual assessment on a small sample in a previous study, that this NID could be found in 28 % of the drivers concerned [START_REF] Orriols | Prescription Medicines and the Risk of Road Traffic Crashes: A French Registry-Based Study[END_REF]. This database contains details about injurious road traffic crashes like on the circumstances of the crash or the vehicles and persons involved.

National police database of injurious road traffic crashes

It can contain the severity of the driver's injuries: unhurt, slightly injured, seriously injured (hospitalized for more than 24 hours) or killed (died within the 30 days following the crash) because the police may have to go to the hospital in order to carry out further investigations.

It also contains data on alcohol consumption. Indeed, police forces test every driver involved in an injurious crash for the presence of alcohol, using a breathalyzer. In the case where the test is positive (≥0.5 g/l), the driver refuses or is incapable of doing the breath test, the blood alcohol level is measured.

When the breath test is negative, the driver is recorded as not being under the influence of alcohol.

National health-care insurance database (HCI database)

The HCI database (Système National d'Informations Inter Régimes de l'Assurance Maladie) covers the entire French population. Each time a medicine is delivered to an outpatient in a pharmacy, a record is added, containing different information like the NID, the day of dispensing and the seven-digit code that identifies the medicine prescription. The HCI database also records 30 long-term chronic diseases (for which patients' health care expenses are fully reimbursed) by their ICD-10 code (International Classification of Diseases, Tenth Revision) and by the dates of the beginning and the end of the disease.

Participant inclusion

In order to mitigate the impact of previous crashes on medicine exposure, drivers were censored at their first involvement in a road traffic crash if they were involved in more than one crash during the study period.

The drivers included in the study are those aged between 18 and 65 years old, who were involved in a road crash in France between July 2005 and December 2015. Subjects whose police reports did not contain their national ID or responsibility were not included.

Medicines and exposure periods

Daily exposure to drug prescription was estimated for ten medicines: drugs used in diabetes, antihypertensives, muscle relaxants, analgesics, antiepileptics, anti-parkinson drugs, psycholeptics (containing all benzodiazepines), psychoanaleptics, other nervous system drugs and antihistamines for systemic use.

Drivers for whom treatment had been delivered on the day of the accident were not considered to be exposed to this treatment because the prescription of this drug could be a consequence of the crash.

Concomitant medicine exposure

In France, a four-level risk classification system of medicines affecting driving abilities is established, ranging from level 0 (no or negligible risk) to level 3 (major risk) (32,33).

Comparisons were adjusted for exposure to other medicines classified in the highest risk levels (levels 2 and 3) used the day of the crash.

Determining driver crash responsibility

A method adapted from Robertson and Drummer [START_REF] Robertson | Responsibility analysis : A Methodology to Study the Effects of Drugs in Driving[END_REF] was used to determine responsibility levels related to the crash. This method is frequently used in the field of road safety. Six factors who are likely to reduce the responsibility of the driver were considered: road, vehicle and driving conditions, type of crash, traffic rule obedience, and difficulty of task involved. Each driver was then assigned a score for each of these factors. This score went from 1 (favorable for driving) to 4 (not favorable for driving). The scores were then added together to obtain a final score. If the sum was between 8 and 15, the driver was considered to be responsible or partially responsible and if it was greater than or equal to 15 then he was considered not responsible.

Then, the principle of responsibility analysis was to compare the probabilities of exposure on the day of the crash between non-responsible drivers and those deemed responsible for the crash. It allowed to include responsible and non-responsible drivers who were directly selected from the driver population.

Statistical analysis

Descriptive analysis

Drivers were classified according to the type of journey when the crash occurred (private, commuting, mission). Using proportions and Chi-2 squared tests, we compared age, gender, socioeconomic category, injury severity, responsibility, vehicle type, time of day, location, alcohol level, frequencies of exposure to medicines of interest and pictogram's levels.

Responsibility analysis

A cases-control study approach was developed, where cases were responsible drivers involved in road crashes and the controls were non-responsible drivers involved in road crashes with two assumptions: The first one is that non-responsible drivers are representative of the driving population (reference population). The second one is that responsible and non-responsible drivers were both driving a vehicle on the day of the crash so they had the same probability of being on the road at the same time (especially in the case were several vehicles were involved). It also allowed providing as much information on others exposures and risk factors for both categories.

The association between driver responsibility in the crash and exposure to drugs used in diabetes, antihypertensives, muscle relaxants, analgesics, antiepileptics, anti-parkinson drugs, psycholeptics, psychoanaleptics, other nervous system drugs (which include parasympathicomimetics, drug used in addiction disorders, antivertigo preparations and other nervous system drugs) and antihistamines for systemic use by the type of journey was then investigated using univariate and multivariate logistic regression models with a stratification on the type of journey.

Interaction tests (likelihood ratio tests) were performed to determine if the type of journey did significantly modify the effect of medicines on responsibility.

Multivariate models included adjustment variables identified in the literature as potential risk or confounding factors: socio-economic characteristics (age, gender, socioeconomic category), health characteristics (exposure to another drug of level 2 or level 3 on the day of the crash), alcohol level and characteristics of the crash (year, month, day of week, time of day (daytime vs nighttime), location (urban area vs nonurban area), vehicle type, blood alcohol level, and injury severity.

An adjustment on the existence of the long-term chronic disease "epilepsy" was also performed for the model about antiepileptics use in order to look at the potential confusion effect of the disease. Indeed, epileptic seizures could be the triggering factor for the crash, independently of the effects of the drugs while antiepileptics are widely used for other indications as neuropathic pain, bipolar disorder or migraine. On private journeys, 1,644 drivers were exposed to antiepileptic drugs. Among them, 138 had a long-term chronic disease related to epilepsy. It was the case for 43 drivers on commuting journeys (398 drivers were exposed to antiepileptics) and 15 drivers on mission journeys (147 drivers were exposed to antiepileptics).

Missing data only concerned blood alcohol concentration (n = 22,225, 12.4 %). We estimated that it was not possible to make imputations on these missing values and consequently created a specific category "missing data". Data were analysed using R® statistical software package, version 3.2.3 and 3.5.0 (R foundation for statistical computing, Vienna, Austria).

Ethical statement

Confidentiality was ensured by using the anonymization function of the Health Care Insurance system [START_REF] Trouessin | FOIN: A nominative information occultation function[END_REF]. The study was approved by the French data protection authority (Commission Nationale Informatique et Libertés).

Results

Study population

Police reports were available for 229,962 persons involved in road crashes with injuries during the study period. Among them, 22,196 were excluded because they were under 18 years of age or over 65 years of age and 28,465 NIDs were excluded because they were pedestrians or passengers. 179,301 drivers were then included.

Responsibility in the crash was established for 179,269 drivers: 85,684 were identified as responsible and 93,585 as non-responsible (figure 1).

Drivers aged between 35 and 45 years were more represented among drivers on commuting (25.5%) or mission (27.8 %) journeys compared to drivers on private journeys (21.7 %).

The proportion of men was significantly higher among drivers on mission journeys (84.3 %) than among drivers on commuting (66.1%) or private journeys (67.4 %) (table 1).

Exposure to Medicines

Fourteen percent of the drivers included in the study (n = 25,466) were exposed to at least one prescription medicine on the day of the crash.

This proportion was about fifteen percent (n = 17,746) for private journeys, about thirteen percent (n = 5,400) for commuting journeys and about eleven percent for mission journeys (n = 2,320).

Standardized incidences on sex and age groups were also different depending on the type of journey with an incidence higher for private journeys (15,942 for 100,000 road crashes) than for commuting [START_REF] Bramness | An increased risk of road traffic accidents after prescriptions of lithium or valproate?[END_REF]963 for 100,000 road crashes) or mission [START_REF] Sundelin | Epilepsy, antiepileptic drugs, and serious transport accidents: A nationwide cohort study[END_REF]636 for 100,000 road crashes) journeys.

Drivers on private journeys had consumption prevalences twice as high as drivers on mission journeys for four classes of medicines: antiepileptics (1.4 % on private journeys and 0.7 % on mission journeys), psycholeptics (6.4 % on private journeys and 2.8 % on mission journeys), psychoanaleptics (4.8 % on private journeys and 2.2 % on mission journeys) and other nervous system drugs (1.3 % on private journeys and 0.7 % on mission journeys).

Consumption prevalences were similar among drivers regardless the type of journey for five medicines: antihypertensives, muscle relaxants, analgesics, anti-parkinson drugs and antihistamines (table 2).

Almost sixty-four percent (n = 16,775) of drivers exposed to one of the tenth medicines studied on the day of the crash were exposed to a level 2 treatment and more than five percent (n = 1,412) to a level 3 treatment (table 3).

After adjusting for risk factors identified from the literature (age, gender, socioeconomic category, blood alcohol concentration, year, month, day of the week and time of day of the crash, localization, vehicle type and injury severity) and for exposure to other level 2 or 3 drugs on the day of the crash, the risk of being responsible for a road traffic crash was significantly increased for drivers on private journeys for four medicines: antiepileptics (N03), psycholeptics (N05), psychoanaleptics (N06) and other nervous system drugs (N07).

For drivers on commuting journeys positive associations were found for three medicines: antiepileptics (N03), psycholeptics (N05) and psychoanaleptics (N06) and for drivers on mission journeys for three treaments: antiepileptics (N03), psychoanaleptics (N06) and other nervous system drugs (N07) (table 4).

After taking account of the existence of a long-term chronic disease for epilepsy, the association of the responsibility and the exposure to antiepileptics were lower but still significant for private journeys (OR=1.39 [1.17-1.64]) and for commuting journeys (OR=1.43 [1.08-1.89]) whereas the association was no more significant for mission journeys (OR= 1.31 [0.81-2.11]). In the same time, the association with the existence of a long-term chronic disease for epilepsy were significant (OR= 1.54 [1.15-2.06] for private journeys, OR= 2.94 [1.73-4.99] for commuting journeys and OR= 3.49 [1.43-8.52] for mission journeys).

The results of the responsibility and interaction tests analysis revealed a significantly higher risk for drivers on commuting and mission journeys than for drivers on private journeys for four medicines: antiepileptics, psycholeptics and psychoanaleptics for drivers on commuting journeys and antiepileptics, psychoanaleptics and other nervous system drugs for drivers on mission journeys. Indeed, for antiepileptics, the risk was 1.47 [1.25-1.73] for drivers who were on private journeys, 1.63 [1.24-2.15] for drivers who were on commuting journeys and 1.59 [1.01-2.51] for drivers who were on mission journeys.

For psycholeptics low differences were observed between private and commuting journeys.

The risk reached 1.17 [1.08-1.26] on private journeys and 1.19 [1.03-1.39] on commuting journeys. The association was not significant for drivers on mission journeys.

For psychoanaleptics, the risk was lower for drivers on private journeys, 1.26 [1. 14-1.40] journeys, and similar among drivers on commuting (OR = 1.37 [1.17-1.60]) and mission (OR = 1.35 [1.02-1.78]) journeys. Surprisingly, for private journeys, drivers exposed to antihistamines were less likely to be responsible for a road crash than unexposed drivers (OR = 0.89 [0.81-0.97]). Exposures to drugs used in the treatment of diabetes or hypertension, muscle relaxants, analgesics or anti-parkinson drugs were not significantly associated with an increase or a decrease in the risk of being responsible for a road crash, regardless of the type of journey performed.

Discussion

Lower medicines consumption prevalences were identified for drivers who were on commuting or mission journeys. On the other hand, slightly higher risks of being responsible for the crash were observed for drivers who were on commuting or mission journeys for four medicines: antiepileptics, psycholeptics (only for commuting journeys, not significant for mission journeys), psychoanaleptics and other nervous system drugs (only for mission journeys, not significant for commuting journeys). These differences may be explained by several factors. First, it is possible that the effects highlighted for some medicines, particularly antiepileptic drugs, are more related to the expression of the pathology itself (epileptic seizures) than to the treatment taken. As shown before, results of responsibility analysis did suggest a confusion effect of epilepsy after adjusting on the existence of a long-term chronic disease (epilepsy) for drivers. Nevertheless, these results are based on low numbers (196 drivers under anti-epileptics were identified as having a long-term chronic disease related to epilepsy (138 on private journeys, 43 on commuting and 15 on mission).

Epilepsy is one of the pathologies that can contraindicate driving vehicles. These findings may suggest that it is probably not antiepileptics but the disease itself or related seizures that are likely to represent a risk factor for a car crash. Furthermore, it can possibly be because of the use of antiepileptics for other diseases than epilepsy [START_REF] Molero | Associations between gabapentinoids and suicidal behaviour, unintentional overdoses, injuries, road traffic incidents, and violent crime: Population based cohort study in Sweden[END_REF]. Further studies should be made in order to investigate this hypothesis.

Anyhow, the risk found in our study for antiepileptics is partially consistent with other studies. Indeed, one previous study performed within the framework of CESIR project [START_REF] Orriols | Prescription of antiepileptics and the risk of road traffic crash[END_REF] found a positive and significant association between the exposure to antiepileptics and the risk of road crash (OR = 1.74 [1.29-2.34]) whereas other studies did not find any significant association [START_REF] Sundelin | Epilepsy, antiepileptic drugs, and serious transport accidents: A nationwide cohort study[END_REF][START_REF] Bramness | An increased risk of road traffic accidents after prescriptions of lithium or valproate?[END_REF].

The slightly higher increase in risk for crashes in mission and commuting may reflect the fact that there may be more drivers with severe epileptic conditions driving for professional sake.

Drivers with a higher risk of seizure could be more easily be deterred to drive for private reasons than when facing the need to exercise a professional activity in order to have an income.

In this study, the most important association on mission journeys has been highlighted for other nervous system drugs (OR = 2.04 [1.35-3.07] on mission journeys against 1.43 [1.21-1.70] on private journeys and 1.33 [0.99-1.77] on commuting journeys). 87,5% of these drivers were males, 35,6% were between 25 and 34, 32,7% between 35 and 44 and 20,2% between 45 and 54. 31,7% were professional drivers. 26,9% drove a car, 22,1% a commercial vehicle and 17,3% a truck. Other nervous system drugs contain parasympathicomimetics, drug used in addiction disorders, antivertigo preparations and other nervous system drugs.

Those treatments are used in order to treat specific pathologies so we can make the assumption that the risk of being responsible for a car crash is, like for antiepileptics, maybe linked to the pathology itself more than to the treatment.

Medicine exposure was ascertained from computerized records of reimbursed prescriptions filled at the pharmacy. It allowed to have data that was not subject to underreporting, which is a real problem often encountered when people self-report medicine exposure information.

On the other hand, the HCI database does not provide information about compliance or selfmedication. Non-compliance can sometimes induce exposure misclassification that would lead to underestimate the risk. Another limitation with this database is that real consumption dates were not available and we had to impute them from dispensing dates.

Our exposure measure should therefore be considered as a surrogate for actual consumption.

The responsibility analysis is a real strength of the study because it allowed using nonresponsible drivers, who are sharing common characteristics with responsible drivers as control group (being a driver, road type, hour of the day…). Furthermore, the responsibility analysis is based on two main assumptions. The first one is that non-responsible drivers are representative of the driving population. The second one is that responsible and nonresponsible drivers were both driving a vehicle on the day of the crash so it also allowed providing the same information on other exposures for both categories (especially for accidents where several vehicles were involved). This would not have been the case in a case-control study where controls were only selected from medical or from driver's license databases.

Finally, it is a method frequently used in the field of road traffic safety. This method has been approved by an independent expert evaluation of responsibility in a previous study [START_REF] Laumon | Cannabis intoxication and fatal road crashes in France: population based case-control study[END_REF] performed on the impact of illegal drug consumption, using this method and the same national database (but limited to fatal crashes). However, this method does not capture the risk, for non-responsible drivers, of being unable to avoid a crash that could be related to medicine use. This would lead to risk under-estimation. However, as responsibility studies typically report the increased odds of responsibility associated with some factor and as non-responsible drivers are considered to be comparable to non-accidented drivers, this study design could lead to interpretational bias and over estimation of the odds ratio when the responsibility odds ratio are considered as direct road crash odds ratio [START_REF] Rogeberg | A meta-analysis of the crash risk of cannabis-positive drivers in culpability studies-Avoiding interpretational bias[END_REF].

Seriously injured drivers were more likely to be taken to the hospital and so to have their national id number registered in police reports. Consequently, there was an association between injury severity and being part of the study and drivers who were more seriously injured and implicated in more serious traffic crashes were slightly over-represented in our sample. Killed or uninjured drivers had lower inclusion rates than severely injured driver.

In addition, the databases used for this study only include drivers involved in road crashes resulting in personal injury (to the driver or to someone else). Drivers who were involved in non-injury road crashes are not included. It could have led to a selection bias.

Other studies used the responsibility analysis approach in order to look at the relation between the exposure to some medicines classes and the risk of being responsible for a road crash [START_REF] Orriols | Road traffic crash risk associated with prescription of hydroxyzine and other sedating H1-antihistamines: A responsibility and case-crossover study[END_REF][START_REF] Barbone | Association of road-traffic accidents with benzodiazepine use[END_REF][START_REF] Lowenstein | Drugs and traffic crash responsibility: A study of injured motorists in Colorado[END_REF]. Nevertheless, to our knowledge, no study assessed the consumptions and risks differences between drivers on private journeys and drivers on professional journeys. However, several studies assessed the consumptions and risks differences between professional and nonprofessional drivers. We found no association between the exposure to antihistamines during occupational journeys and the risk of being responsible for a road traffic crash unlike a previous study conducted among professional drivers (OR = 3.44 [1.06-11.16]). However, this latter study included truck drivers only.

The risks identified in our study for psycholeptics and psychoanaleptics are consistent with other previous studies not conducted but no data have been so far published on professional drivers in particular [START_REF] Leveille | Psychoactive medications and injurious motor vehicle collisions involving older drivers[END_REF][START_REF] Ray | Psychoactive drugs and the risk of injurious motor vehicle crashes in elderly drivers[END_REF][START_REF] Bramness | Minor increase in risk of road traffic accidents after prescriptions of antidepressants: A study of population registry data in Norway[END_REF].

Conclusion

This study brings original results about the problematic of medicines consumptions and driving for purposes of work. Firstly, in terms of the prevalence of the exposure that is higher for drivers on private journeys. Then, in terms of the risk of being responsible for a road traffic crash linked to those consumptions, with a risk higher for drivers on commuting or mission journeys compared to private journeys for four medicines: antiepileptics, psycholeptics (only for commuting journeys, not significant for mission journeys), psychoanaleptics and other nervous system drugs (only for mission journeys, not significant for commuting journeys).

This study brings new information to health professionals (general practitioners, occupational practitioners and practitioners approved for driving licenses), especially about antiepileptics, psycholeptics, psychoanaleptics and other nervous system drugs consumptions. For those four medicines, an increased-risk was observed for drivers on commuting or mission journeys.

Preventive measures could be implemented in order to make occupational practitioners and employees aware of the occupational road risks associated with the use of some medications.

Further studies to confirm these results would be useful, including studies on the time to onset of work-related road crashes after the start of a drug prescription. This would clarify the impact of drugs and improve prevention and recommendations currently implemented by, for example, adjusting the time to return to driving.

  Details collected by police forces on each injurious crash are stored in the IC database (Bulletins d'Analyse des Accidents Corporels [BAAC]). This standardized database is led by the National Interministerial Road Safety Observatory (Observatoire national interministériel de sécurité routière [ONISR]). It is filled after police reports.

Figure 1 :

 1 Figure 1: Flowchart of the inclusion procedure.

Table 1 : Baseline characteristics according to the type of journey, N = 179,269.

 1

	Police reports
	652,663 reports
	(image files)	Police national database of injurious crashes
		745,200 records (including 1,661,609 individuals)
	305,964 NID collected in police reports	-1,239,324 drivers -422,285 other users (pedestrian or passengers)
	229,962 NIDs matched
	(drivers or other users)
		-22,196 drivers aged between 18 and 65 years
		old
		-28,465 other users (pedestrians or passengers)
	179,301 NIDs matched
	(drivers)
		32 drivers for which responsibility
		could not be determined
	179,269 drivers responsibility has
	been determined
	85,684 responsible	93,585 non-responsible

Acknowledgements

Competing interests:

The authors declare that they have no competing interests.