domain: sdv info.info-mo spi

text:
Etienne Descamps
email: etienne.descamps@grignon.inra.fr

Nathalie Perrot
email: nathalie.perrot@grignon.inra.fr

Ioan-Cristian Trelea

Sebastien Gaucel
email: sebastien.gaucel@grignon.inra.fr

Alain Riaublanc
email: alain.riaublanc@nantes.inra.fr

Alan Mackie
email: alan.mackie@bbsrc.ac.uk

Evelyne Lutton
email: evelyne.lutton@grignon.inra.fr

Ioan Cristian Trelea
email: cristian.trelea@agroparistech.fr

Coupling deterministic and random sequential approaches for structure and texture prediction of a dairy oil-in-water emulsion

Keywords: milk gel model, dairy products, oil-water interface, non linear optimisation, model coupling, knowledge integration

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

Surface-active molecules, such as proteins, polymers, ionic and non-ionic surfactants play a major role in the stabilisation of dispersed system such as oil in water emulsions [START_REF] Dickinson | Milk protein interfacial layers and the relationship to emulsion stability and rheology[END_REF]. In food systems, an important class of emulsifiers are proteins: they are adsorbed on the oil droplet surface during homogenisation. Stabilisation is a consequence of the ability of the proteins to generate repulsive interactions (steric and electrostatic) between oil droplets [START_REF] Mcclements | Protein-stabilized emulsions[END_REF].

The emulsifying properties of milk proteins are excellent and justify their wide use in food processing [START_REF] Dickinson | Caseins in emulsions: interfacial properties and interactions[END_REF]. Milk proteins are divided in two major categories: caseins (as casein micelles (CM) or individual caseins) and native whey proteins (WP).

During most of milk gels processing, the milk proteins solution undergoes a heat treatment. Usually, a heat treatment denatures soluble proteins and aggregates of proteins appear, with a strong impact on the physicochemical properties of the final product. Above 70 • C, whey proteins are partly denatured and form aggregates (WPA) while casein micelles are less sensitive to heat treatment but can form complexes with WP (CMWP). This last reaction is in competition with the whey proteins aggregation [START_REF] Guyomarc'h | Formation of soluble and micelle-bound protein aggregates in heated milk[END_REF]. All these phenomena lead to four potential types of proteins in the solution: (CM), (WP), (WPA) and (CMWP).

Several studies on milk gels [START_REF] Dickinson | Milk protein interfacial layers and the relationship to emulsion stability and rheology[END_REF][START_REF] Murray | Interfacial rheology of food emulsifiers and proteins[END_REF][START_REF] Gaygadzhiev | Influence of the emulsion droplet type on the rheological characteristics and microstructure of rennet gels from reconstituted milk[END_REF][START_REF] Knudsen | Droplet surface properties and rheology of concentrated oil in water emulsions stabilized by heat-modified beta-lactoglobulin B[END_REF][START_REF] Dickinson | Mixed biopolymers at interfaces: Competitive adsorption and multilayer structures[END_REF] showed the importance of the dynamics and competition between these types of particles taking place at a fluid-fluid interface of lipid droplets during the emulsion. Although studies describe the structuring of pure whey protein aggregates (WPA) submitted to heat [START_REF] Rabe | Understanding protein adsorption phenomena at solid surfaces[END_REF]. Nevertheless less is known about the behaviour of complex aggregates made of casein micelles and whey proteins (CMWP) [START_REF] Morand | Increasing the hydrophobicity of the heatinduced whey protein complexes improves the acid gelation of skim milk[END_REF]. Moreover, the data and expertise collected on complex aggregates is generally difficult to integrate in existing models [START_REF] Foucquier | Toward an integrated modeling of the dairy product transformations, a review of the existing mathematical models[END_REF]. Mixed solutions are not thermodynamically controlled processes of competitive adsorption, and cannot be predicted from any classical model like the Langmuir one [START_REF] Rabe | Understanding protein adsorption phenomena at solid surfaces[END_REF][START_REF] Dickinson | Caseins in emulsions: interfacial properties and interactions[END_REF]. In this context, stochastic approaches simulating the layer of adsorbed protein at the nanooscopic level are relevant, but need important computing time if several types of particles are considered. The surface-active molecules, i.e. the CMWP aggregates, are in competition with the pure whey protein (WPA) aggregates at the oil in water interface and play a major role in the network created at a higher level. As highlighted by Dickinson [START_REF] Dickinson | Mixed biopolymers at interfaces: Competitive adsorption and multilayer structures[END_REF], the interpretation of surface composition in emulsion containing the full range of aggregated milk proteins is quite complex and certainly not yet fully understood. Modelling is a good candidate to help to understand such systems.

Experiments have been previously performed, generating complex unexpected behaviour, for a same total concentration of proteins, when using different mass ratio of (CM). A wide range of interfacial composition have been generated when changing this initial condition. The work presented in this paper starts from those observations.

We propose a computational approach to simulate the structure of a dairy emulsion, from a given composition and process condition. The computational approach is validated on the available dataset. The model uses a stochastic approach to generate at the nanoscopic scale a random spatial configuration of the physical elements of the emulsion. The purpose of the model is to predict the emergence of the macroscopic structure from the local organisation at the nanoscopic scale. This approach is strongly inspired by the Random Sequential Adsorption (RSA) model [START_REF] Feder | Random sequential adsorption[END_REF] and its derivatives [START_REF] Nasir | Sequential and competitive adsorption of bovine serum albumin and βlactoglobulin, and their resistance to exchange with α-lactalbumin and β-casein[END_REF], with the particularity that the interface is not represented as a plane, but as a 3D surface made of randomly distributed fat droplets in the space.

Additionally, our model is able to manage the competition of elements of different sizes, which is not considered in previous RSA models.

The paper is organised as follows. After a description of the experimental data used to optimise and validate the model (section 2), the model is developed in section 3. Results are presented in section 4: after a sensitivity analysis, the parameters of the model are fitted using an evolutionary optimisation approach (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. Discussion, conclusions and future work are finally developed in sections 5 and 6.

Experimental data

Various emulsions were generated with a range of controlled interfacial compositions (weight ratio of casein micelles (CM) to whey proteins (WP) from 80:20 to 12:88). Two sets of experiments were made:

M1

M2

Model of denaturation and aggregation of whey proteins

M2a

Model of adsorption at interface

M2b

Model of link creation between droplets Both sets of experiments were made with the same preheat treatment and protein solutions, but with different experimental devices and volumes. The protein phases of both experimental datasets were made of the same powders, with the same ionic strength, pH and ion composition.

M2 parameters M1 parameters

The two emulsion processes were based on the same principle: the continuous phase of the emulsion was formed from milk proteins dissolved in permeate. These milk proteins were a mixture of caseins (Promilk 852B, IDI company, France with 5% moisture, 1.5% fat, 85.5% nitrogenous matter/dry matter, 8.5% mineral matter, 4% lactose, 81% nitrogenous matter (on powder), 92% casein micelle, 2.6% Ca, 1.5% P, 0.3% K, 0.1% Na and 0.1% of Mg) and native whey proteins (BiPro, DAVISCO company, Minnesota with 5% max of moisture, 95% min. of protein, dry basis, 1% max fat, 3% max ash, 1% max lactose, a pH between 6.7 and 7.5) with milk permeate powder (Armor protéines, France with a pH of 6.0 min, 3% max moisture, 3% min proteins, 1% max fat, 82% lactose, 8% ashes). The continuous phase was prepared the day before use, was stored at 4 • C and then heated at 80 • C. The dispersed phase of the emulsion is made of saturated lipid: anhydrous milk fat (AMF) heated at 60 • C to become liquid. These two phases were then homogenised in order to get an emulsion.

The processes are different for the homogenisation and volume of resulting emulsion. For Database 1, the blending was done with a rotor stator (Polytron, Heidolph Silent Crusher M), in a low pressure homogeniser (Stansted Fluid Power, Stansted, UK) at 50 bars, whereas experiments in Database 2 were made with a blender (BL450 series, KENWOOD) with a shearing cycle (30 seconds low speed, 30 seconds of rest, then 2 x 30 seconds high speed) and in a manual homogeniser (EmulsiFlex -B3, AVESTIN) using 6 passes at 20 x 200 PSI.

In Database 1, the emulsion mass was 70 g (49 g of continuous phase and 21 g of dispersed phase) and in Database 2, the emulsion mass was 260 g (182 g of continuous phase and 78 g of dispersed phase).

In order to evaluate the impact of initial conditions on the structure and texture of the emulsion, experiments were carried out with various initial conditions (Table 1). The following initial conditions were kept the same for every experiment:

• For Database 1: the pre-heat treatment temperature of milk proteins was 80 • C, the denaturation level, i.e. the proportion of denatured WP in the solution, was around 0.6 and the mass of lipid was 21 g, for a total emulsion volume of 70 mL.

• For Database 2, the pre-heat treatment temperature of milk proteins was 80 • C and the mass of lipid The following measurements were collected for characterising the emulsions at a micro/nanoscale.

• Diameter and size distribution of lipid droplets. Laser light scattering was used to measure the diameter of the lipid droplets in the emulsion and evaluate the size distribution. In Database 1, measurements were made using a Saturn DigiSizer 5200 (Micromeritics, Norcross, USA). These measurements enabled us to calculate the initial free lipid surface S 0 . In Database 2, measurements were made using an LS 13320 Laser Diffraction Particle Size Analyser (Beckman Coulter). The measurement error for all considered devices is around 10%.

• Diameter and size distribution of whey protein aggregates. Dynamic light scattering was used to measure the diameter of the whey protein aggregates in the emulsion and evaluate the size distribution.

In Databases 1 and 2, these measurements were made using a Nanosizer ZS (Malvern Instrument, UK).

• Interfacial concentration (Γ) and percentage of adsorbed caseins (w cas ads ). In Database 1, the study of the composition of the interface of lipid droplets was carried out with the Patton and Huston technique [START_REF] Patton | A method for isolation of milk-fat globules[END_REF] which is used to separate droplets. The protein concentration at the interface was then quantified by the Markwell method [START_REF] Markwell | Modification of Lowry procedure to simplify protein determination in membrane and lipoprotein samples[END_REF]. SDS-PAGE electrophoresis was then used to determine the concentration of each protein at the interface. In Database 2, SDS-PAGE electrophoresis was carried out and gels were purchased from Invitrogen Ltd. (Paisley, UK). The measurement error of the Patton techniques and electrophoresis were each around 10%.

• Sensory measurements. Two groups of sensory texture were described during experiments for each conditions: gelified or liquid. Measurements were performed [START_REF] Curt | Formalization of at-line human evaluations to monitor product changes during processing : the concept of sensory indicator[END_REF] on the emulsion 2 hours after the processing at 20 (2) Viscoelastic measurements performed with a controlled stain rheometer (ARES, TA Instruments, New Castle, Delaware, USA) using a plate-plate geometry (diameter 40 mm, gap 100µm) at 20 • C. A dynamic mechanical spectrum (G', storage modulus and G", loss modulus as a function of the angular frequency) was obtained in a frequency range of frequency from 0.01 to 100 rad/s, at 10 % shear strain.

Coupled models for a milk gel emulsion

The model presented in this paper is a coupling between a first order differential kinetic model of protein denaturation (M1) and a stochastic model of simulation of the fat droplets interface colonisation (M2a) and bridge creation between fat droplets (M2b) (Figure 1). Coupling M1 and M2 allows to simulate the emergence of a network at a mesoscale level from local droplet considerations. Both models are presented below.

M1: First order differential kinetic model of protein denaturation

The first model aims at simulating the denaturation process, in particular the competitive reaction that occurs during the thermal denaturation of whey proteins. According to the values of the parameters of M1, denatured WP are more or less associated with CM or aggregated in the form of WPA. The outputs of M1 are the mass concentration of WPAs and the mass concentration of CMWPs. The assumptions are the following:

ρ wpcm (∞) ρ wp (0)
Ratio of denatured WP that associates with CM

• According to [START_REF] Tolkach | Reaction kinetic pathway of reversible and irreversible thermal denaturation of β -lactoglobulin[END_REF], the denaturation reaction of the β-lactoglobulin is a kinetic reaction of order 1.5(equation 1). This reaction is known as following an Arrhenius law (equation 2) [START_REF] Tolkach | Optimization of thermal pretreatment conditions for the separation of native α-lactalbumin from whey protein concentrates by means of selective denaturation of β-lactoglobulin[END_REF]. All native whey proteins are supposed to have the same kinetic as β-lactoglobulin.

• The reactions of aggregation and denaturation have the same order of 1.5. The two reactions are merged in one [START_REF] Chantoiseau | Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an eulerian/lagrangian approach[END_REF] (equation 3).

• The association step is a first order reaction with respect to the concentrations of casein micelle and denatured whey protein that are not yet aggregated. A denatured WP can join either the surface of a CM or a WPA. The affinity of association between WPs and CM is considered as constant. The possible association between an aggregate and casein micelles is neglected (equation 4).

• The simulated process is simplified: the temperature is considered constant along the heat treatment; denaturation is supposed to stop when the experimental measured denaturation rate is reached.

The full model is then made of a set of three differential equations (Figure 2):

• The denaturation reaction given by [START_REF] Tolkach | Reaction kinetic pathway of reversible and irreversible thermal denaturation of β -lactoglobulin[END_REF]:

dρ wp* dt = k d (T ) × ρ wp 1.5 (1)
k d (T ) follows an Arrhenius law according to temperature with known parameters [START_REF] Tolkach | Optimization of thermal pretreatment conditions for the separation of native α-lactalbumin from whey protein concentrates by means of selective denaturation of β-lactoglobulin[END_REF]:

k d (T ) = k d 0 exp {-E a d/RT } (2)
In our experiments, T is fixed to 80 • C, k d is also fixed until the experimental denaturation rate is reached.

Denaturation of native whey proteins

Aggregation of denatured whey proteins

Association of casein micelles with denatured whey proteins

k d k k cm ag (1) (2)
(3) • For the aggregation reaction, a similar law is used with an order of reaction of 1.5 [START_REF] Tolkach | Reaction kinetic pathway of reversible and irreversible thermal denaturation of β -lactoglobulin[END_REF][START_REF] Chantoiseau | Coupling fluid flow, heat transfer and thermal denaturation-aggregation of beta-lactoglobulin using an eulerian/lagrangian approach[END_REF]:

(CMWP) (CM) (WPA) (WPA) (WP) (WP*) (WP*) (WP*)
dρ wpa dt (t) = k ag ρ 1.5 wp* (3)
• The association between casein micelles and non aggregated denatured whey proteins is modelled using a linear dependant law:

dρ wpcm dt (t) = k cm ρ cm ρ wp* (t) (4)
In this model, k d is used as a reference for the other reactions, k ag and k cm , and fixed to 1. Indeed, only equilibrium values predicted by M1 are used by M2, which allows to fix one of the reaction constants arbitrarily.

M2 : Adsorption at interface and bridging between fat droplets 3.2.1. Model structure

The model is based on a three dimensional space representation. This 3D space is a cube representing a random distribution of the fat droplets. It has been calibrated using d 3,2 and d 4,3 , measured for each experiments (Table 1) as explained in section 3.2.2 below. For the adsorption process, a competition between WPA, WP and CM (including CMWP) is considered. This competition depends on the current interfacial concentration of adsorbed native WP at the surface of the droplets.

The model is made of two compartments (Figure 1). The first compartment simulates the adsorption and gives the final interface configuration at the local level of each fat droplet (M2a in Figure 1). The second compartment uses the space configuration to estimate the ability of each adsorbed element to link with another fat droplet (M2b in Figure 1).

A stochastic approach using Brownian dynamics is used to determine the potential of bridging between fat droplets (M2b).

Implementation of the adsorption model (M2a)

The following assumptions have been made, in order to reduce the unknown parameter space and the complexity of the model:

• A 1 :
The desorption phenomenon is neglected. The adsorption of the WP and of some casein proteins are known to be irreversible [START_REF] Dickinson | Competitive adsorption and protein-surfactant interactions in oil-in-water emulsions[END_REF]. As a consequence, we consider the WP aggregation and MC adsorption as irreversible too.

• A 2 : The native WPs can move on the interface [START_REF] Brun | Some effects of heat on the competitive adsorption of caseins and whey proteins in oil-in-water emulsions[END_REF]. To take into account this phenomenon, the model processes the native WPs in a different way than the other particles. Indeed, it only takes into account the surface occupied by all adsorbed native WPs, and does not allocate them a fixed position at interface. All further adsorbed elements are supposed to be able to push the whey proteins on the interface. WPs are managed as a population contrarily to the other particles in competition (WPA, CM and CMWP), where interaction between individuals is taken into account. A droplet only contains the information of the mass of native WPs adsorbed at its surface. Due to the high concentration of ionic elements present in the solution, the electrostatic repulsion is neglected.

• A 3 : The coalescence phenomenon is simplified, and obeys to an inverse exponential law. The homogeniser is supposed to generate an initial population of fat droplets with a mean radius of 100 nm. This population then converges to a distribution that corresponds to the experimental values of d 4,3 and d 3,2 . Coalescence is represented by a mean field law.

• A 4 : The size of the aggregates of whey proteins has been measured during the experiments : it is considered as fixed with radius s wpa = 65 nm in the given conditions [START_REF] Surel | Composition and structure of interface impacts texture of o/w emulsions[END_REF]. It has been reported that the fractal dimension of these aggregates is D wpa = 2.4 [START_REF] Verheul | Power law behavior of structural properties of protein gels[END_REF] and that the size of a native whey protein is around s wp = 1.5 nm [START_REF] Aymard | The effect of temperature and ionic strength on the dimerisation of β-lactoglobulin[END_REF], which allows to compute the weight of each aggregate (equation 9).

• A 5 : The size of casein micelles have also been measured with a size of s cm = 75 nm and a mass of m cm = 5.77 • 10 -16 g.

• A 6 : Spreading is a deformation of the particle at the surface of the fat droplet after adsorption. It is modeled as an instantaneous phenomenon whose caracteristic parameter is the increase factor of the apparent radius of the particle.

The model M2 simulates a fraction of volume V of the emulsion. Inside this volume, a population of fat droplets is placed. The droplets sizes are randomly generated using the experimental measures (distribution of fat droplets and fat volume fraction). The size distribution is defined with two parameters, d 4,3 and d 3,2 :

d 4,3 = i d 4 i i d 3 i ; d 3,2 = i d 3 i i d 2 i ( 5
)
The parameters of a log normal distribution are straighforwardly calculated from these two factors. Suppose the diameter distribution follows a ln N(µ, σ 2 ) law, d k then follows a log normal law ln

N(k • µ, k 2 • σ 2 ).
d 4,3 = e 4•µ+4 2 • σ 2 2 e 3•µ+3 2 • σ 2 2 = e µ+ 7•σ 2 2 ( 6
)
d 3,2 = e 3•µ+3 2 • σ 2 2 e 2•µ+2 2 • σ 2 2 = e µ+ 5•σ 2 2 (7)
µ and σ are obtained:

σ = ln( d 4,3 d 3,2 ) ; µ = ln(d 3,2 ) - 5 • σ 2 2 (8)
Knowing V, [CM], ρ wpa (∞) and ρ wp (∞), the numbers of particles N cm , N wpa and N wp are computed as follows.

N wpa is determined using the measured fractal dimension:

N wpa = ρ wpa (∞) m wpa × V with m wpa = m wp × s wpa s wp D wpa (9)
The number of native whey proteins, and casein micelles are computed in a similar way:

N wp = ρ wp (∞) m wp × V ; N cm = ρ cm m cm × V (10)
Starting from a random fat droplet space distribution with an empty interface, the model consists in fixing elements from the solution to a randomly selected location on the interface. This loop proceeds as follows:

• A fat droplet d is chosen according to a random law: the probability is proportional to the droplet total surface.

• For coalescence, instead of simulating a population of fat droplets that are progressively merged, the droplets are supposed to be made of several sub-droplets of radius r sub . Each droplet generated by the random space distribution corresponds to a set of sub-droplets of radius r sub with same cumulated volume. The radius follows an inverse power law r d sub (t) = r 0 + (r dr 0 ) • e -t . The surface is given by

s d (t) = N d sub (t) × 4π r d sub (t) 2 with N d sub (t) = 3 × v d × 4π r d sub (t) 3 -1
the mean number of sub droplets. Coalescence is considered as uniform for all droplets.

t is initialised to 0, and increased by an increment ∆t at each iteration:

1 ∆t = F cm × N cm + F wpa × N wpa ×        d s d (t)        + F wp × N wp ×        d (s d (t) -s d occ )        (11)
• If the WPs cover the whole available surface of the fat droplet (equivalent to Γ > Γ wp = 3.2mg.m -2 ), no other particle can be adsorbed. Another droplet is chosen.

• The choice of a type of particle is governed by a random law where the probability to choose an element is proportional to the affinity parameter of each type (F cm , F wpa and F wp ).

For each S ∈ {cm, wpa, wp} : proba

S = N S × F S U∈{cm,wpa,wp} N U × F U (12)
The position of the particle on the fat droplet surface is fixed as follows:

-If the particle is a WPA or a CM (including CMWP): its position on the surface is randomly chosen. If there is no collision with elements of the list of already fixed particles, the particle is added to the list attached to the fat droplet. The collision test is a geometrical computation based on the position and radius (after spreading) of each particle of the list.

-If the particle is a native WP: the collision test is equivalently replaced by a Bernoulli process, whose parameter is the ratio between the available surface and the total surface of the droplet.

The previous loop ends when a maximum number of unsuccessful trials (i.e. without fixed element) is reached, or when all fat droplets are saturated. The flow diagram of the M2a algorithm is presented Figure 4.

Bridging between fat droplets (M2b)

The estimation of the links is based on a stochastic process that considers each element capable of binding two droplets (i.e. whey protein aggregates and activated casein micelles CMWP). Casein micelles are considered to be activated (i.e. CMWP) as soon as they have fixed enough whey proteins on their surfaces. Model M1 gives an average mass of denatured whey proteins associated with casein micelles. The model introduces a normal distribution of the average mass of WP associated with a CM in order to avoid a threshold effect between (CM) and (CMWP) states:

ρ cm (∞) × V N cm
Using the dispersion factor γ cmwp , it is possible to estimate how many CM's have more than µ cmwp × m CM associated WP and become activated.

For each linkable particle of radius r already adsorbed by a droplet, the following steps are performed:

• all droplets able to bind with the considered particle are sorted from the nearest to the farthest,

• for each droplet of this list, a random number x is chosen in (0, K LINK × r). If the distance from the fat droplet to the considered particle is less than x, a link is created.

Outputs of the model

The model allows the calculation of three outputs:

• Γ, the interfacial concentration, given by the total mass of adsorbed particles including the associated WP divided by the total surface of all fat droplets.

• w cm ads , the ratio between the adsorbed caseins and the total amount of adsorbed proteins.

• G link , the average number of links by droplet.

Results

This section presents simulations for the model M with tests on an experimental data set. A threshold phenomenon has been observed in the [20%, 30%] range of w cm 0 [START_REF] Surel | Heating milk proteins impact interface proteins adsorption and texture of emulsion : a way to innovative textures[END_REF][START_REF] Surel | Composition and structure of interface impacts texture of o/w emulsions[END_REF]: There actually exists a complex balance between competitive reactions occurring in thermal treatment and complex dynamics of colonisation at the fat droplet interface during homogenisation. Understanding this behaviour is an open question.

A sensitivity analysis has been first performed, then a parameter optimisation has been run for estimating the optimal parameters of M. If quantitative experimental data are available on interfacial concentration and casein mass ratio at interface, data about numbers of links between particles are not available. This is why sensitivity analysis and parameter optimisation have not been performed for the parameters and outputs of M2b. Nevertheless it has been shown that the number of links is directly related to rheological properties measured during the experiments [START_REF] Dickinson | Structure and rheology of simulated gels formed from aggregated colloidal particles[END_REF]. Two parameters out of 9 of the models M1 and M2a are fixed: r 0 , the radius of a fat droplet at the homogeniser exit, is set to 100 nm, in concordance with the experimental observations. Γ wp is fixed to 3.2 mg.m -2 , in accordance to the literature [START_REF] Hunt | Adsorption behaviour of whey protein isolate and caseinate in soya oil-inwater emulsions[END_REF], for single layer colonisation.

The sensitivity analysis and the parameter optimisation are made on parameters of M1 and M2a only. A logarithmic scale is used for the multiplicative parameters (that are here k {cm,ag} for M1, and F {cm,wp,wpa} for M2). For the sensitivity analysis, these parameters are arbitrarily bounded between e -7 and e 7 (approximatively 0.001 and 1000), but they are not limited during optimisation. The spreading parameters (Ksp cm , Ksp wp ) are bounded between 1 (no spreading) and 4 (the particle takes 16 times its initial surface, which is equivalent to an occupied surface of 100 m 2 .g -1 for a WPA and 150 m 2 .g -1 for a CM, to be compared to the values founded section 4.2.1).

Four inputs have been considered for the sensitivity analysis: w cm (between 0% and 100%), c prot (between 40 g.L -1 and 60 g.L -1 in the given experimental conditions), d 3.2 (between 0.5 and 1 µm) and (d 4.3d 3.2 ) (between 0 and 0.5 µm).

Sensitivity analysis and parameter optimisation take as outputs the interfacial concentration and the casein mass percentage.

Sensitivity analysis

Methodology

The model is assimilated to a stochastic function f : I × Ω → O with: • Ω the sample space of the random generator of the model,

• I = (P × U) ⊂ R
• O ⊂ R 2 the output space for the model.

|P| = 50 reference points have been chosen using a random Latin hypercube distribution in the space I. A stochastic exploration is then performed for each reference point and for each dimension of I. A squared uniform random law is used. 60 points are tested, so as we get 60 × 11 simulations for each reference point. As a whole, we get 60 × 11 × 50 simulations.

The chosen perturbation is x p i = x i + 2 × r 2 × X i max -X i min , with x p i a perturbation on x i and r a random number in [-1 ; 1].

The median of 8 repetitions for each output is considered in the sensitivity analysis. The data are then analysed using a three order polynomial least square regression. Each polynomial allows to evaluate one of the partial derivatives:

∂y j ∂x i
From the points of the latin square, three subsets S are considered: w cm 0 < 30%, 30% < w cm 0 < 60%, w cm 0 > 60%.

Dependence d i, j (S ) between input i and output j is evaluated by the sum of the partial derivatives on each subset S : 4: Impact of each parameter on each output of M2a, expressed as a percentage of the total perturbation: () < 5%, (+) < 10%, (++) < 15%, (+ + +) < 20%, (+ + ++) < 30%, (+ + + + +) > 30%

d i, j (S ) =         (x,y)∈S ∂y j ∂x i         ×         i (x,y)∈S ∂y j ∂x i         -1
The sensitivity analysis have been performed on the computer grid MIGALE of the INRA Jouy-en-Josas. It has taken 3 days using 100 parallel jobs.

Table 4 presents the influence of parameters on the model response. M1 has a stronger influence when the initial concentration of caseins is low. This is due to the high concentration of native WPs that could be denatured in this case. F cm and F wp have more impact on w cm ads when 30% < w cm 0 < 60%. The impact of F wpa is more important when w cm 0 < 30%, but remains quite low. A possible reason could be due to the boundaries chosen for M1 parameters, that could disadvantage WPA formation. Ksp cm and Ksp wpa are key parameters for the model. They have a strong impact on both Γ and w cm ads .

Identification of parameters by optimisation

Evolutionary optimisation

The problem of identifying an optimal parameter setting for a model according to an experimental dataset is a complex optimisation problem. The coupled model presented above has an obviously non linear and stochastic behaviour, this is a reason why a robust stochastic optimisation has been used. The algorithm used in this section is an Evolutionary Algorithm (EA), which corresponds to a large set of techniques that rely on the computer simulation of natural evolution mechanisms (Artificial Darwinism) [START_REF] Fraser | Simulation of genetic systems by automatic digital computers. i. introduction[END_REF][START_REF] Bremermann | Optimization through evolution and recombination, Self-Organizing Systems[END_REF][START_REF] Holland | Outline for a logical theory of adaptative systems[END_REF][START_REF] Rechenberg | Evolutionsstrategie : Optimierung Technicher System nach Prinzipien der Biologischen Evolution[END_REF][START_REF] Bäck | An overview of evolutionary algorithms for parameter optimization[END_REF].

EAs copy natural evolution with very simple stochastic operators. According to Darwin [START_REF] Darwin | n the Origin of Species by Means of Natural Selection, or the Preservation of Favoured Races in the Struggle for Life[END_REF], adaptation is based on three mechanisms: random variations, inheritance, and survival / reproduction of the fittest individuals. Used for computational optimisation, this scheme has the major advantage to make only few assumptions on the function to be optimised (there is no need to have a continuous or derivable function for instance). In short, Evolutionary Algorithms consider a population of potential solutions exactly as a population of individuals of a natural population that live, fight and reproduce. The environment pressure is replaced by an "optimisation" pressure: the function to be optimised, the fitness, is considered as a measure of the adaptation of the individual to its environment (fitness). In this way, individuals that reproduce are the best ones with respect to the problem to be solved.

Search space and fitness function

The parameter estimation is turned into an optimisation problem as follows: the optimal parameter set is searched in the space of all possible parameter sets. The function to be optimised, the fitness function, is an estimation of the quality of the model: a least square fitting of the predicted values (interfacial concentration and mass percentage of caseins at interface) on the learning sample (Database 1). As the model is stochastic, • S = (x i , y i ) 7 i=1 ⊂ U × O the data set (Database 1). • z i = median j M((p, x i ), s j ) the estimation of each point by simulating model M.

Algorithm

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] has been used. CMA-ES is considered as one of the best real-value optimizers for single-objective problems, delivering high-quality results in a very limited amount of time. Its population size was set to λ = 256, and log scale weights have been used for the update of µ (the number of offspring's generated at each generation). Search was initialised at the following point: (k cm = 1m 1.5 .g -0.5 , k ag = 1,

F cm = 1, F wp = 1,F wpa = 1, Ksp cm = 2, Ksp wpa = 2)
During model development, non-linear optimisation methods have been tested (including CMA-ES, NSGA-2 and Monte-Carlo random search), and it appeared that CMA-ES converged faster than the others to one of the best fitting parameter sets. Due to the heavy computationnal nature of the model (depending to the parameters it takes from 20 min to 5 hours to evaluate one set of parameters on a recent core), only two run of CMA-ES have been made for the last revision of the model.

Figure 5 shows the fitness evolution of the two run of CMA-ES that have been made with the latest revision of the model. The longest run have taken 3 days and 7 hours on a computational server of 64 cores (with 8 Xeon(R) CPU E7-8837 processor at 2.67GHz of 8 core) to make 54 generations. Due to heavy computational issues the program have been written using C++ and Fortran program languages. The CMA-ES implementation used is from [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

The best global fitness has been found for: k cm = 4.66m 1.5 .g, k ag = 7.6•10 -3 , F cm = 0.13, F wp = 7.0•10 -4 , F wpa = 0.16, Ksp cm = 1.64 and Ksp wpa = 2.14.

• The difference of magnitude between k cm and k ag is due to their units (Table 2). It is then difficult to make a comparison between these reactions.

• F wp is lower than the other F. Recall the probability to choose one particle (CMs, WPAs, WPs) is proportional to the product of their number and the corresponding F coefficients (equation 12). The small value of F wp , is compensated by the large number of WP.

• Surface occupied by CMs and WPAs are computed using the spreading parameters Ksp cm and Ksp wpa . These values can be compared to the spreading of native WPS calculated from litterature data ( [START_REF] Dalgleish | Adsorption of protein and the stability of emulsions[END_REF], [START_REF] Hunt | Adsorption behaviour of whey protein isolate and caseinate in soya oil-inwater emulsions[END_REF]): around 1500 312 = 2.19. Organized structures like casein micelles or whey protein aggregates have been shown to spread in a less extent than native whey proteins.

• Using the identified spreading parameters, the specific occupied surface are equivalent to 26 m 2 .g -1 and 29 m 2 .g -1 for the CMs and WPAs respectively. The maximum surface occupied by a single layer of native WPs is fixed in the model to 312 m 2 .g -1 [START_REF] Hunt | Adsorption behaviour of whey protein isolate and caseinate in soya oil-inwater emulsions[END_REF]. The calculated parameters are compared to data measured during the experiments on a Langmuir balance. Surface areas occupied by 1 g of the various protein species were estimated in real emulsion systems produced with pure proteins. The apparent surface area occupied by 1 g of the different protein species was estimated by the specific surface of the emulsion (6/D32) divided by the protein load at the interface of droplets measured according to the method described by Patton and Huston [START_REF] Patton | A method for isolation of milk-fat globules[END_REF]. It was done for solutions of casein micelles and whey proteins near the ratios manipulated in this paper and for a same mean fat droplet size but at an air in water interface. It gives an order of magnitude for a target of behavior, even if it cannot be used for a quantitative validation. The values found for s wp and s cm were equal respectively to 473 m 2 /g and 98 m 2 /g, which is in accordance with the order of magnitude of the parameters we found.

Validation of the model

The stochasticity of the model is presented in Figure 6. The variability induced by this stochasticity is obviously higher for the interfacial concentration than for the casein percentage at interface.

The variability of the interfacial concentration is particularly high for initial percentage of caseins between 20 and 50 %. It is in this range that the competition between particle is the highest and can mostly impact the amount of available surface on the fat droplets. This might be enhanced by the fact that, according to table 1, the difference between d 4.3 and d 3.2 is higher in this case. This induces a more dispersed distribution of droplet sizes.

The percentage of casein at interface is less impacted by the stochasticity of the model than the interfacial concentration.

Figure 7 and Table 4.3 show the behaviour of the model fitting to experimental data. Generally, the model has a good fit for the two outputs of the system and for the two databases tested (Database 1 for the parameter estimation and Database 2 for the test). A mean error is calculated :

• For the output w cm ads : 9.3 % and 7.3 % respectively for Database 1 and Database 2 This is comparable to the experimental measurement error of 10%. However two points on the first Database 1 are not well predicted: one for an interfacial concentration at 21% of initial casein rate and the other for the percentage of casein at 32% of initial casein rate. Moreover, a weakness of the model for low values of w cm 0 is observed. It might be due to a change of behaviour at the interface. This lack of fit reflects the complexity of the interfacial organisation. The data set shows a complex balance between competitive reactions occurring in thermal treatment and complex dynamics of colonisation at the fat droplet interface during homogenisation [START_REF] Dickinson | Milk protein interfacial layers and the relationship to emulsion stability and rheology[END_REF]. An example of prediction of the spatial structure of the gel by M2b is presented Figure 8. It is reached for a non fitted set of parameters (K link max = 5, µ cmwp = 0.6). It seems qualitatively coherent with the sensory measurements achieved during the experiments [START_REF] Surel | Heating milk proteins impact interface proteins adsorption and texture of emulsion : a way to innovative textures[END_REF][START_REF] Surel | Composition and structure of interface impacts texture of o/w emulsions[END_REF]. Indeed we know that the number of links is directly related to rheological properties ( [START_REF] Dickinson | Structure and rheology of simulated gels formed from aggregated colloidal particles[END_REF]). As a consequence, it seems possible to have an indirect measurement of sensory characteristics through an evaluation of the number and the nature (WPA or CMWP connection) of created links.

For low values of w cm 0 , (Figure 8a) solutions were evaluated as firm gelified structures. It is correlated with the number of links predicted by the model (here 842). Same results are observed for w cm 0 = 49% (Figure 8b), with two kinds of links: one created by CMWPs (here 859) and another by WPAs (here 244). The type of links involved in the structure has probably an impact on the texture. Nevertheless, it can not be depicted using a global sensory analysis. For w cm 0 = 80% (Figure 8c), the structure appears completely different, with a number of links radically reduced. It is coherent with liquid perception of the gel for this initial ratio.

Discussion

The model is able to fit relatively well the experimental data except for some of the points in Figure 7. In order to improve the results several ideas can be considered: • The organisation of the interface is currently assumed to be in a single layer. In reality it has been suggested that the interface might be more complex when several types of particles are in competition [START_REF] Dickinson | Milk protein interfacial layers and the relationship to emulsion stability and rheology[END_REF]. Nevertheless few quantitative knowledge is available.

• The stochastic representation of the competition between casein micelles and whey protein aggregates could be improved by adding knowledge on the colonisation phenomena.

• Another further development could be a model that simulates the coalescence between fat droplets as a stochastic process instead of the mean field approach.

Nevertheless the above improvements would increase the computation complexity of the model which is already high. The adequate model complexity is a question of balance between the computation burden and the physical reality.

Conclusion and future work

This combination of a deterministic and a stochastic model shows relevant results in replicating the behaviour of the complex artificial concentrated milk protein mix under study. Such a model can help us to have a better understanding of the system through simulation analysis. Further studies will be done along the directions described in the previous section.

Acknowledgements

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under the grant agreement FP7-222 654-DREAM.

Figure 2 :

 2 Figure 2: The three steps of a denaturation reaction.

Figure 3 :

 3 Figure 3: Coalescence and adsorption at the exit of the homogeniser

  d 4,3 and d 3,2 are determined from the mean values of the log normal distribution e µ+ σ 2 2 :

Figure 4 :

 4 Figure 4: Flow diagram of the algorithm for the adsorption part of the stochastic model (M2a)

Figure 5 :

 5 Figure 5: Fitness value during optimisation

Figure 6 :

 6 Figure 6: Model behaviour of 16 different runs that point out the stochasticity of the model

Figure 7 :Figure 8 :

 78 Figure 7: Comparison between experimental data (from Databases 1 and 2) with the model behaviour when using the best parameter set

Table 1 :

 1 Initial conditions and measurement results for Databases 1 and 2.

	d 3.2 : Surface area mean diameter of fat droplets,
	d 4.3 : Volume mean diameter,
	w cm 0 : Initial percentage of caseins in the solution,
	w cm ads : Percentage of adsorbed caseins,
	c prot : Protein concentration in water phase.
	Γ: Interfacial concentration.
	w cas 0 d3.2 d4.3 (%) (µm) (µm) (g.L -1 ) c prot	w cas ads (%)	Γ (mg.m -2 )
			Database 1
	13	0.5	0.8	48.4	9	7.1
	19	0.45	0.7	48.8		4.4
	21	0.5	0.8	48.9		3.9
	26	0.4	0.7	48.9		3.4
	32	0.45	1.1	49.4		5.7
	49	0.57	0.9	49.7		7.2
	80	0.8	1.0	50.4		6.1
			Database 2
	13	0.76 1.19	48.3	0	7.79
	31	0.94 1.55	47.5	4	5.48
	49	0.94 1.60	49.7		8.04
	80	0.86 1.43	50.4		7.22

Table 2 :

 2 Model M1 Glossary

	Name	Description	Unit
		Parameters
	k d	Denaturation rate of native WP	m 1.5 .g -0.5 .s -1
	k ag	Aggregation rate of denatured WP	m 1.5 .g -0.5 .s -1
	k cm	Association rate between denatured WP and CM	m 3 .g -1 .s -1
		Inputs
	η	Denaturation level of WP	%
	ρ wp (0) + ρ cm	Initial proteins mass concentration	g.m -3
	w cm 0	Initial mass percentage of casein	%
		State variables
	ρ cm	CM mass concentration	g.m -3
	ρ wpcm	Mass concentration of WP associated with CM	g.m -3
	ρ wpa	WP aggregates mass concentration	g.m -3
	ρ wp	WP natives mass concentration	g.m -3
	ρ wp*	Free denatured WP mass concentration	g.m -3
		Outputs of the model

Table 3 :

 3 Model M2 Glossary

	Name	Description	Unit
		Constants
	s wpa	WPA radius	nm
	s wp	Native WP radius	nm
	m wp	Native WP mass	g
	s cm	CM radius	nm
	m cm	CM mass	g
	D wpa	Fractal dimension of the WPAs	*
		Parameters of M2a
	r 0	Radius of the fat droplets at the exit of the homogeniser	nm
	Ksp cm Ksp wpa	Maximal spreading of CM Maximal spreading of WP aggregates	* *
	Γ wp	Maximal concentration of a one layer native WP interface	mg.m -2
	F cm	Adsorption rate of CM (including CMWP)	*
	F wp	Adsorption rate of natives WP	*
	F wpa	Adsorption rate of WP aggregates	*
		Parameters of M2b
	K lnk wpa	Link factor of WP aggregates	*
	K lnk cmwp	Link factor of activated CM	*
	µ cmwp	WP mass activation threshold of a CM	g.g -1
	γ cmwp	Variance of WP association on CM	g.g -1
		Inputs
	V	Volume of the space simulated by the model	m 3
	N cm	Number of CM	*
	N wpa	Number of WP aggregates	*
	N wp	Number of natives WP	*
	(G i ) k 1	Fat droplets population present in the space
	m cmwp	Mass mean of WP associated on CM	g.g -1
		State variables
	(P i ) p 1	Set of adsorbed elements including type and possible geometric
		information
	(L i ) q 1	Set of links between particles and fat droplets
		Outputs of the model
	Γ	Interfacial concentration	mg.m -2
	w cm ads	Mass percentage of casein adsorbed at the interface	%
	G link	Average number of links by droplet	*
	p lnk CMWP	Percentage of CMWP in the links set	%
	(*) dimensionless.

Table

			4.1.2. Results of the sensitivity analysis
	X	w cm 0 < 30%	30% < w cm 0 < 60%	w cm 0 > 60%
	Parameter w cm ads	Γ	w cm ads	Γ	w cm ads	Γ
	k ag	++	+ + +	+		+
	k cm	+ + + + + ++	++	+	++
	F cm	+ + +	++	+ + +	+	++	+
	F wp	++	++	+ + ++	++	++	+
	F wpa	++		+
	Ksp cm	++	+ + +	++	+ + + + + + + + + + + + + + +
	Ksp wpa	+ + +	+ + +	+ + ++	++	+ + ++	++

Table 5 :

 5 Model evaluation using the best fitted parameter of M1 and M2a.

	w cm 0		w cm ads	Γ
	exp exp model	exp	model
	%	%	%	mg.m -2 mg.m -2
			Database 1
	13	9	18	7.10	4.63
	19	13	26	4.40	4.29
	21	22	29	3.90	4.94
	26	41	35	3.40	4.07
	32	68	41	5.70	5.15
	49	61	59	7.20	6.35
	80	83	84	6.10	6.27
			Database 2
	13	0	18	7.79	6.06
	31	40	38	5.48	6.15
	49	53	56	8.04	6.25
	80	78	84	7.22	6.44

The authors would like to thank Claire Surel, Marc Anton from INRA BIA Nantes, and Marine Rouland from IFR Norwich, for their precious help in building of the model and providing experiments. Computational resources have been provided by ISC-PIF and INRA-Jouy-en-Josas (MIGALE grid).

• To build the model, a first database was collected from experiments conducted at the pilot plant of INRA BIA (Institut National de Recherche Agronomique, Biopolymères -Interactions -Assemblages, Nantes, France) [15].

• To validate the model, a second database was collected from experiments conducted in the , Norwich, England) [16].


Marion M Doyennette

Isabelle Deleris Déléris
email: isabelle.deleris@grignon.inra.fr

Gilles Feron
email: gilles.feron@dijon.inra.fr

Elisabeth E Guichard
email: elisabeth.guichard@dijon.inra.fr

Isabelle I Souchon
email: isabelle.souchon@grignon.inra.fr

Ioan-Cristian I.-C Trelea

G Féron

I C Tréléa

Main individual and product

Keywords: dynamic model, aroma compounds, food oral processing, physiology, mass transfer 1

Introduction

The mechanistic modelling of aroma compound release during food consumption allows one to calculate, from known physical laws, the amount of aroma compounds transferred over time in each anatomical compartment involved during food oral processing (mouth, nasal cavity, pharynx).

The first mechanistic models have focused on predicting the release of aroma compounds from a two-phase emulsion (water, oil) in contact with gas [START_REF] Harrison | Mathematical models of flavor release from liquid emulsions[END_REF]. They were based on physicochemical principles governing the release of volatile molecules from a food matrix: (i) the mass conservation of volatile compound, (ii) the mass transfer at the emulsion-gas interface (interfacial penetration theory), (iii) equilibrium properties at the emulsion-gas interface [START_REF] Overbosch | Flavor release in the mouth[END_REF].

First-order chemical kinetics have also been included in some models to describe reversible interactions between aroma compounds and non-volatile compounds such as macromolecules [START_REF] Harrison | Mathematical models of flavor release from liquid emulsions[END_REF]. However, these first models are not really representative of in vivo phenomena because the geometry of the system (surfaces and volumes) is assumed constant (which is not the case during food consumption) and they do not consider dynamic phenomena such as the dilution with salivary flow or the cyclic breathing of the subject. In addition, some of these models do not include any comparison with experimental data release [START_REF] Harrison | Mathematical models of flavor release from liquid emulsions[END_REF].

Further development of these pioneering models has led to a better representation of in vivo conditions occurring during food consumption, by including notably parameters related to individual physiology, such as salivary flow, periodic breath, etc. [START_REF] Normand | Modeling the Kinetics of Flavour Release during Drinking[END_REF][START_REF] Wright | Modelling flavour release from a chewed bolus in the mouth: Part II. The release kinetics[END_REF]. In addition, one of the major improvements in these models was to consider the aroma persistence phenomenon i.e. aroma release from bolus deposit covering the pharyngeal mucosa after swallowing. These models showed the relative roles of product and consumer characteristics and were validated against experimental data.

The first model including physiological data was proposed by Normand et al. [START_REF] Normand | Modeling the Kinetics of Flavour Release during Drinking[END_REF] in the case of liquid products and highlighted two main release regimes: (i) release due to equilibrium batch extraction (only pertinent for few breaths after swallowing) and (ii) release from lubricated mucosa (persistence phenomena). In the case of semi-liquid food, the most comprehensive model to date is the one developed by Tréléa et al. [START_REF] Tréléa | Mechanistic mathematical model for in vivo aroma release during eating of semi-liquid foods[END_REF] and further developed by Doyennette et al. [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF] coupling aroma release in the mouth and in the pharynx.

Mechanistic models describing aroma release from masticated foods are far less available in the literature, mostly because of the difficulty to understand the complex mechanisms which are involved during the consumption of those foods. Compared to liquid products, mass transfer occurs through several interfaces (product/saliva, saliva/air) in "solid" matrices (needing mastication) [START_REF] Harrison | Computer simulation of flavor release from solid foods in the mouth[END_REF][START_REF] De Roos | Non-equilibrium partition model for predicting flavour release in the mouth[END_REF]. Also, additional phenomena have to be considered: product dissolution and melting due to intraoral manipulations (chewing, saliva incorporation, warming) and the generation of a dynamic exchange interface between product and saliva.

Existing models describing aroma release from solid food are not complete since they do not take breathing into account [START_REF] Wright | Modelling flavour release from a chewed bolus in the mouth: Part II. The release kinetics[END_REF][START_REF] Hills | Two-film theory of flavour release from solids[END_REF]. Furthermore, the assumptions used are not always transferable to other types of food matrices and to associated food oral processes than those studied (e.g. candy) or to other consumption patterns (e.g. sucking mechanism studied by Wright et al. [START_REF] Wright | Modelling flavour release from a chewed bolus in the mouth: Part II. The release kinetics[END_REF]). [START_REF] Hills | Two-film theory of flavour release from solids[END_REF] also highlighted the importance of modelling the variation in the contact surface between the product and saliva, whose change in time plays a central role in the aroma compounds release. In their model, they suggest a multi-fragmentation theory to represent the chewing of a candy (assimilated to a cube). This approach provides a law of time change in the contact surface between product and saliva, using only two parameters: the number of bites and the duration of chewing. One might expect, however, that chewing real products, like cheeses for instance, will not always produce cubes and multiply the area by a factor of 2 after each bite. Other studies have also investigated the fragmentation of solids placed in the mouth under the action of chewing. The degree of fragmentation of a food can be determined by different types of experimental measurements such as the analysis of particle size distribution in the bolus [START_REF] Peyron | Particle Size Distribution of Food Boluses after Mastication of Six Natural Foods[END_REF]. These experimental data allow the determination of the laws of fragmentation depending on the type of food: deterministic analytic laws of product fragmentation over time [START_REF] Tobitsuka | Retention of a European Pear Aroma Model Mixture Using Different Types of Saccharides[END_REF], or probabilistic laws [START_REF] Harrison | Computer simulation of flavor release from solid foods in the mouth[END_REF][START_REF] Agrawal | Mechanical properties of foods responsible for resisting food breakdown in the human mouth[END_REF][START_REF] Van Der Bilt | A mathematical description of the comminution of food during mastication in man[END_REF][START_REF] Wright | Modelling flavour release from a chewed bolus in the mouth: Part I. Mastication[END_REF]. While these approaches can be very comprehensive, they have several limitations. For example, some of them require complex experimental protocols that are difficult to implement (bolus spitting after a variable number of bites). Other approaches need the determination of many parameters (such as the number of chews, the number or the size of food particles after each chew, etc…), usually unknown, and which may depend on the product and/or on the individual.

Hills and Harrison

From this literature review, it appears that modelling the release of aroma compounds during solid food consumption remains a challenging task and that the results of published studies are difficult to extrapolate to other experimental conditions. In this context, the present study proposes a model simulating the release of aroma compounds, applicable to different food matrices and different individuals with a wide range of physiological characteristics. With this mechanistic model, our main objective was to understand the mechanisms and parameters governing the release of aroma compounds during the intra-oral manipulation and swallowing of "solid" food (needing chewing).

Mathematical modelling of in vivo aroma compound release

To do so, simulations issued from the model were compared to in vivo release data of two aroma compounds, measured by atmospheric pressure chemical ionization mass spectrometry (APCI-MS) on ten panelists, during the consumption of four sorts of cheese, varying in composition and texture. Model assumptions are presented and their validity is discussed.

2.1.Principles of the model

The aroma release model described in this study applies to food products requiring mastication. It is an extension of the model developed for liquid and semi-liquid foods by Doyennette et al. [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF]. The new insight is the consideration of the mastication process. The global eating process involves several steps shown in figure 1: the initial state of the system (product introduction in mouth), the intra-oral manipulation of the product, which consists in several masticatory cycles, one or more swallowing events and the resting phase (release during the post-mastication stage) which occurs when there is no more product in the mouth.

The intra-oral manipulation phase usually lasts until the first swallow. However, in some cases (particularly for firm products), it can extend beyond the first swallow.

Similarly to the model of Doyennette et al. [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF], the present description of the swallowing step includes simultaneous contractions of the oral cavity and of the pharynx, leading to air and product expulsion, followed by relaxation and filling with fresh air. This will cause the expulsion of the bolus formed in the mouth into the pharynx. Each swallowing step leads to the deposit of a small part of the in-mouth liquid phase of the bolus on the pharyngeal walls.

A residual amount of the solid part of the bolus remains in the mouth and is chewed again and mixed with saliva in order to form a new food bolus suitable for swallowing.

A schematic representation of the four compartments involved in the model design, as well as their connections and the mechanisms responsible for flavour release, are given in figure 2.

All variables and parameters of the interconnected compartments required for the model simulation are specified in this figure .    Compared to the previous model [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF], two main differences can be observed: first, the presence of three instead of two sub-compartments in the mouth (the non-dissolved food product, which is the solid part of the bolus, the liquid phase of the bolus, made of saliva and dissolved food product and the air phase) and secondly, the opening of the velopharynx during chewing, allowing the transport of aroma compounds between the air phases of mouth and pharynx (figure 2).

Concerning the food bolus fragmentation under mastication, the particle size distribution in the bolus after each bite has already been described in literature [START_REF] Harrison | Computer simulation of flavor release from solid foods in the mouth[END_REF][START_REF] Peyron | Particle Size Distribution of Food Boluses after Mastication of Six Natural Foods[END_REF][START_REF] Agrawal | Mechanical properties of foods responsible for resisting food breakdown in the human mouth[END_REF][START_REF] Van Der Bilt | A mathematical description of the comminution of food during mastication in man[END_REF][START_REF] Wright | Modelling flavour release from a chewed bolus in the mouth: Part I. Mastication[END_REF]. However, these approaches require the knowledge of parameters such as food particle size after each chew, which is difficult to determine experimentally, particularly for pasty products like cheese used in the present study. The relevant parameter for aroma compound release is not the particle size itself, but the contact area between the solid and the liquid phase of the bolus.

In this study, we focus on the generation of this area, as previously done in literature [START_REF] De Loubens | Mechanistic model of in vitro salt release from model dairy gels based on standardized breakdown test simulating mastication[END_REF].

2.2.Mathematical description of mass transfer in mouth during the eating of masticated foods

In the following paragraphs, the main differences with the previously published model for liquid products [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF] are emphasized.

The solid food product placed in the mouth is broken down during intra-oral manipulation.

Two concurrent phenomena can occur:

• the transfer of aroma compounds from the non-dissolved product into the liquid phase of the bolus,

• the melting (dissolution) of the product in the liquid phase of the bolus due to the combined action of the mastication, the saliva incorporation and the warming of the product in the mouth. This also leads to the release of aroma compounds contained in the dissolved product towards the liquid phase of the bolus.

In practice, it is very difficult to distinguish between the relative contributions of each mechanism. Moreover, the phenomena of transfer and of dissolution can be described by very similar equations. Considering the studied products (cheeses), we arbitrarily chose to assign the release of aroma compounds into saliva to a single mechanism (dissolution) while being aware that the transfer also contributes to this process.

Air/bolus interfacial conditions in the oral cavity

Similarly to Doyennette et al. [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF], the transfer resistance on the air side 1/𝑘 𝑂𝐴 was assumed to be negligible when compared to the transfer resistance on the product side 1/𝑘 𝑂𝐿 (the orders of magnitude of the mass transfer coefficients, representing the reciprocal of the resistances, are 10 -2 and 10 -6 m/s respectively [START_REF] Cussler | Mass Transfer in Fluid Systems[END_REF].

Therefore, the interfacial aroma compound concentration on the liquid bolus side, using the partition conditions at the interface, is given by:

𝐶 𝑂𝐴𝐿 * (𝑡) = 𝐶 𝑂𝐴 (𝑡) 𝐾 𝑂𝐴𝐿 (𝑡) (1)
The volatile mass flux between the air and the liquid bolus in the oral cavity φ OAL is mainly determined by the resistance located on the bolus side and is given by the difference between the liquid bolus concentration (C OL ) and the interfacial concentration

(C * OAL 𝜙 𝑂𝐴𝐿 (𝑡) = 𝑘 𝑂𝐿 (𝑡) × 𝐴 𝑂𝐴𝐿 (𝑡) × (𝐶 𝑂𝐿 (𝑡) -𝐶 𝑂𝐴𝐿 * (𝑡)) (2)
):

Air in the oral cavity

In addition to the aroma compound flux from the liquid bolus, the air in the mouth can also exchange aroma compounds with the air in the pharynx. Jaw movements during mastication induce velopharynx opening and cyclic air flow between the pharynx and the mouth [START_REF] Matsuo | Effects of Respiration on Soft Palate Movement in Feeding[END_REF].

The variation of aroma concentration in the air in the oral cavity C OA is thus due to the volatile flux from the liquid food bolus and from the air coming from the pharynx (Q OA 𝑉 𝑂𝐴 (𝑡) ×

𝑑𝐶 𝑂𝐴 (𝑡) 𝑑𝑡 = 𝜙 𝑂𝐴𝐿 (𝑡) + � 𝑄 𝑂𝐴 (𝑡) × (𝐶 𝐹𝐴 (𝑡) -𝐶 𝑂𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑂𝐴 (𝑡) ≥ 0 0 𝑖𝑓𝑄 𝑂𝐴 (𝑡) < 0 (3)

≥0

means that the air flows in the direction shown by the arrow in figure 2):

Little is known on the real change in V OA . It was assumed that the masticatory movements create a cyclic variation of the air volume V OA around a mean value V OA mean

𝑉 𝑂𝐴 (𝑡) = 𝑉 𝑂𝐴 𝑚𝑒𝑎𝑛 + Δ𝑉 𝑂𝐴 × sin ( 2 × 𝜋 × 𝑓𝑟 𝑜𝑝𝑒𝑛𝑖𝑛𝑔 × 𝑡) (4)
as it has been observed by Matsuo et al. [START_REF] Matsuo | Effects of Respiration on Soft Palate Movement in Feeding[END_REF]. In this case:

Therefore, the air flow rate coming from the mouth

Q OA 𝑄 𝑂𝐴 (𝑡) = 𝑑𝑉 𝑂𝐴 (𝑡) 𝑑𝑡 = 2 × 𝜋 × 𝑓𝑟 𝑜𝑝𝑒𝑛𝑖𝑛𝑔 × Δ𝑉 𝑂𝐴 × cos ( 2 × 𝜋 × 𝑓𝑟 𝑜𝑝𝑒𝑛𝑖𝑛𝑔 × 𝑡) (5)
is calculated as follows:

It is expected that Δ𝑉 𝑂𝐴 can be highly variable among individuals, but no quantitative data was found in the literature. A fair estimation of Δ𝑉 𝑂𝐴 seems to be 20% of 𝑉 𝑂𝐴 𝑚𝑒𝑎𝑛 , value which was used in our model. Based on observations from Matsuo et al. [START_REF] Matsuo | Effects of Respiration on Soft Palate Movement in Feeding[END_REF], we assumed that fr opening can vary among individuals, within a defined range. We supposed that the highest opening frequency was coordinated with the masticatory frequency (fr masticatory ), and that the lowest opening frequency was coordinated with the respiratory frequency (F R

Product in the oral cavity

). This assumption will be further discussed in the results and discussion section.

In line with the assumption of the dissolution rather than transfer mechanism discussed above, the aroma compound concentration in the solid (undissolved) food product fraction remains constant in time:

𝐶 𝑂𝑃 (𝑡) = 𝐶 𝑂𝑃 (0) (6)
Due to dissolution and fragmentation processes, the volume of the solid food product decreases over time, while its contact area with the liquid phase of the bolus increases and then abruptly decreases due to swallowing.

The dissolution of the product at a rate v gives the following equation:

𝑑𝑉 𝑂𝑃 (𝑡) 𝑑𝑡 = -𝑣 × 𝐴 𝑂𝐿𝑃 (𝑡) (7)
Due to fragmentation induced by mastication, the contact area between the solid product present in the mouth and the liquid bolus increases in time. The exact rate of change of this contact area is not known for pasty products, like cheeses used in this study. In absence of more detailed information we assumed that, as long as some solid product is present in the mouth and a regular mastication behaviour were considered, its contact area with the liquid bolus evolves linearly over time. This simplifying hypothesis might be refined in the future based on more detailed studies, however. With this assumption:

𝑑𝐴 𝑂𝐿𝑃 (𝑡) 𝑑𝑡 = � 𝐴 𝑂𝐿𝑃𝑑𝑒𝑔 -𝐴 𝑂𝐿𝑃𝑖𝑛𝑖 𝑡 𝑑𝑒𝑔 -𝑡 0 if 𝑉 𝑂𝑃 (𝒕) > 𝟎 0 if 𝑉 𝑂𝑃 (𝒕) = 𝟎 (8)
with the index "deg" meaning "at the current deglutition moment", and the index "ini" meaning "at food product introduction in mouth" (at t 0

Liquid bolus in the oral cavity

) or "just after the previous deglutition".

The liquid bolus compartment has a composition which evolves over time. It is initially composed of pure saliva, and is progressively flavoured by the addition of dissolved product.

Its volume increases with the addition of saliva (salivary flow) and with the incorporation of dissolved product, and periodically decreases after swallowing.

The volume of the bolus 𝑉 𝑂𝐿 (𝑡) can be thus divided into two parts:

𝑉 𝑂𝐿 (𝑡) = 𝑉 𝑂𝑆 (𝑡) + 𝑉 𝑂𝑃𝐷 (𝑡) (9)
In that case, we have:

𝑑𝑉 𝑂𝑆 (𝑡) 𝑑𝑡 = 𝑄 𝑂𝑆 (10)
with

𝑑𝑉 𝑂𝑃𝐷 (𝑡)
𝑑𝑡 the product dissolution rate, defined as

𝑑𝑉 𝑂𝑃𝐷 (𝑡) 𝑑𝑡 = 𝑣 × 𝐴 𝑂𝐿𝑃 (𝑡) (11)
The mass balance for a given aroma compound in the bolus leads to the following equation:

𝑑𝑉 𝑂𝐿 (𝑡)×𝐶 𝑂𝐿 (𝑡) 𝑑𝑡 = 𝜙 𝑂𝐿𝑃 (𝑡) -𝜙 𝑂𝐴𝐿 (𝑡) (12)
with the volatile mass flux φ OLP

𝜙 𝑂𝐿𝑃 (𝑡) = 𝑣 × 𝐴 𝑂𝐿𝑃 (𝑡) × 𝐶 𝑂𝑃 (𝑡) (13)
coming from product dissolution:

The mass flux φ OAL is given by the difference between the aroma concentration in the liquid phase of the bolus (C OL ) and the interfacial concentration (C *

OAL

The properties of the liquid bolus relevant for aroma compound transfer, namely the air/bolus partition coefficient (K ) (Eq. 2).

OAL , Eq. 1) and the mass transfer coefficient (k OL 2.3.Mathematical description of aroma release during the pharyngeal step , Eq. 2) change with the relative fraction of the saliva and dissolved product in the liquid bolus. This dependence was previously established [START_REF] Doyennette | The dynamics of aroma compound transfer properties in cheeses during simulated eating conditions[END_REF] and included in the model simulations.

Bolus in the pharynx

Phenomena governing flavour release from the bolus in the pharynx are similar to the ones described for liquid products [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF].

Air in the pharynx

The air in the pharynx exchanges aroma compounds with the bolus in the pharynx, and with the other compartments: the mouth (air flow 𝑄 𝑂𝐴 (𝑡)), the nose (air flow 𝑄 𝑁𝐴 (𝑡)) and the trachea (air flow 𝑄 𝑇𝐴 (𝑡)). Compared to the model of Doyennette et al. [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF], air flows to and from the mouth have been modified (due to the velopharynx opening during mastication). By convention, the air flow rates indicated in figure 2 are positive if they follow the direction of the arrow. The air balance in the pharynx at any time gives the following relationship:

𝑄 𝑁𝐴 (𝑡) = -𝑄 𝑇𝐴 (𝑡) + 𝑄 𝑂𝐴 (𝑡) (14)
For the considered aroma compound, the mass balance in the air of the pharynx gives the equation:

𝑉 𝐹𝐴 × 𝑑𝐶 𝐹𝐴 (𝑡) 𝑑𝑡 = 𝜙 𝐹𝐴𝐿 (𝑡) + � -𝑄 𝑂𝐴 (𝑡) × (𝐶 𝑂𝐴 (𝑡) -𝐶 𝐹𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑂𝐴 (𝑡) < 0 𝑄 𝑁𝐴 (𝑡) × (𝐶 𝑁𝐴 (𝑡) -𝐶 𝐹𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑁𝐴 (𝑡) ≥ 0 𝑄 𝑇𝐴 (𝑡) × (𝐶 𝑇𝐴 (𝑡) -𝐶 𝐹𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑇𝐴 (𝑡) ≥ 0 (15)

2.4.Conditions of the system after each deglutition

The deglutition step is very short compared to product residence time in mouth. It is thus described as quick contraction and relaxation events (figure 1).

Product and bolus in the oral cavity after deglutition

Aroma compound concentrations in the product and in the bolus are unchanged during deglutition. The volume of product in the oral cavity just after deglutition (V OPdeg+ ) decreases and corresponds to the volume of product just before swallowing (V OPdeg-) multiplied by the fraction of liquid bolus remaining in the mouth after deglutition (r L V ), i.e. :

OPdeg+ = V OPdeg-×r L
Similarly, the liquid bolus/product contact area in the oral cavity after swallowing is calculated as:

A OLPdeg+ = A OLPdeg-×r L

The volumes of saliva (V ( 17)

OSdeg+ ) and of dissolved product (V OPDdeg+ V 𝑂𝑆𝑑𝑒𝑔+ = V 𝑂𝑆𝑑𝑒𝑔-× V 𝑆𝑎𝑙𝑖𝑣𝑎𝑑𝑒𝑔+ V 𝑂𝑆𝑑𝑒𝑔-+V 𝑂𝑃𝐷𝑑𝑒𝑔-
, and

V 𝑂𝑃𝐷𝑑𝑒𝑔+ = V 𝑂𝑃𝐷𝑑𝑒𝑔-× V 𝑆𝑎𝑙𝑖𝑣𝑎𝑑𝑒𝑔+ V 𝑂𝑆𝑑𝑒𝑔-+V 𝑂𝑃𝐷𝑑𝑒𝑔- (18)
) in the oral cavity decrease after swallowing. We assume that the bolus liquid part (saliva and dissolved product) is uniformly swallowed. Therefore, we have:

with V Salivadeg+ being the volume of saliva usually present in the oral cavity after swallowing.

Aroma compound retention by lubricated mucosa

To take the reservoir effect of lubricated mucosa into account, additional compartments were included into the model. We assumed that air in the mouth, in the pharynx and in the nose was in contact with lubricated mucosa layers within the corresponding compartment. For example, in the nose, the volatile mass flux φ NAM

𝜙 𝑁𝐴𝑀 (𝑡) = 𝑘 𝑁𝑀 × 𝐴 𝑁𝐴𝑀 × (𝐶 𝑁𝑀 (𝑡) - 𝐶 𝑁𝐴 (𝑡) 𝐾 𝑁𝐴𝑀 ) (19)
between the air and the mucosa is given by:

The mass balance of aroma compound in air contained in nasal cavity leads to:

𝑉 𝑁𝐴 × 𝑑𝐶 𝐹𝐴 (𝑡) 𝑑𝑡 = 𝜙 𝑁𝐴𝑀 (𝑡) + � 𝑄 𝑁𝐴 (𝑡) × (0 -𝐶 𝑁𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑁𝐴 (𝑡) ≥ 0 𝑄 𝑇𝐴 (𝑡) × (𝐶 𝐹𝐴 (𝑡) -𝐶 𝑇𝐴 (𝑡)) 𝑖𝑓 𝑄 𝑇𝐴 (𝑡) ≥ 0 (20)
and in nasal mucosa:

𝑉 𝑁𝑀 × 𝑑𝐶 𝑁𝑀 (𝑡) 𝑑𝑡 = -𝜙 𝑁𝐴𝑀 (𝑡) (21)
Here, V NM 𝑉 𝑁𝑀 = 𝑒 𝑁𝑀 × 𝐴 𝑁𝐴𝑀 [START_REF] Shojaei | Buccal mucosa as a route for systemic drug delivery: a review[END_REF] is the volume of lubricated mucosa which is involved in the interaction with the aroma compound. It can be expressed as Similar equations were added in the mouth and in the pharynx to take the effect of the lubricated mucosa present in these compartments into consideration. In the absence of more detailed information, the contact area between lubricated mucosa and air was arbitrarily set to half of the contact area between the bolus and the air in the mouth, as well as in the pharynx.

Little information is available in the literature concerning values to be considered in the aroma retention model. The contact area between air and lubricated mucosa in the nose A NAM was set to 160 cm² [START_REF] Levitzky | Pulmonary physiology[END_REF]. In literature, values comprised between 5.6×10 -5 and 4.8×10 -1 were mentioned for air/mucus partition coefficient of butanol and octanol in bullfrog [START_REF] Hornung | Olfactory Mucosa-Air Partitioning of Odorants[END_REF]. Thus, a typical value of 1×10 -3 was selected here. Concerning the mucosa layer thickness, Shojaei [START_REF] Shojaei | Buccal mucosa as a route for systemic drug delivery: a review[END_REF] gave values between 500 and 800 µm in mouth and 100 to 200 µm for the gingival mucosa. It is expected, however, that the aroma compound will not necessarily have time to diffuse in the whole epithelium thickness, and these values are considered as upper limits for mucosa layer thickness involved in aroma retention.

Material and Methods

Flavoured cheese products

Four industrial cheese products (melt-cheese technology) varying in composition and texture (two fat levels and two firmness levels) and flavoured with ethyl propanoate and 2-nonanone were studied [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF]. They were coded according to their characteristics, with a combination of letters as following: S or F for Soft or Firm respectively and l or h for low-fat or high-fat respectively.

Panelist selection and their physiological characterization at rest

Ten individuals were selected from a large panel composed of 44 subjects: they were verified to be representative of the whole panel (mean and standard deviation) concerning physiological data, masticatory behavior and APCI-MS release kinetics [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF].

The volumes of oral, nasal and pharyngeal cavities of subjects were measured with the Eccovision Acoustic Rhinopharyngometer from Eccovision (Sleep Group Solutions, North Miami Beach, FL 33162). A software was developed to calculate automatically the air/product areas of oral and pharyngeal cavities for each individual [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF].

The tidal volume of each individual was measured with a spirometer (Pulmo System II, MSR, Rungis, France) [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF]. The respiratory frequency F R 3.2. In-nose measurements of aroma release by Atmospheric Pressure Chemical Ionization-Mass Spectrometry (APCI-MS) and data processing used in the model was calculated directly from the acetone signal measured by APCI-MS.

Aroma release was measured using APCI-MS [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF]. For each product, each panelist and for the three replicates, a mean curve and an envelope curve (representing the standard deviation of the replicates) were determined based on the peak lines (curve linking the maxima of aroma release profile ) [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF]. Release kinetics were divided in two different phases: "phase 1" corresponded to the chewing phase before the first deglutition and "phase 2" corresponded to the rest of the chewing phase (if present) and to the resting phase. The end of phase 2 was set when less than 10% of the maximal intensity was reached. Areas under curve (AUC) were determined for each phase.

Physiological characterization of individuals during cheese consumption

Chewing activity, bolus saliva content and mouth coating were determined for each panelist as described in Repoux et al. [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF].

It has been demonstrated from various studies [START_REF] Watanabe | A Comparison of the Effects of Tasting and Chewing Foods On the Flow-rate of Whole Saliva In Man[END_REF][START_REF] Gaviao | Chewing behavior and salivary secretion[END_REF] that salivary flow during food consumption is much higher than the salivary flow at rest or artificially stimulated. Therefore, an average rate of saliva incorporation in the bolus during food product consumption (Q OS

𝑄 𝑂𝑆 = 𝑝𝑒𝑟𝑐𝑒𝑛𝑡𝑎𝑔𝑒 𝑜𝑓 𝑖𝑛𝑐𝑜𝑟𝑝𝑜𝑟𝑎𝑡𝑒𝑑 𝑠𝑎𝑙𝑖𝑣𝑎 ×𝑎𝑚𝑜𝑢𝑛𝑡 𝑜𝑓 𝑓𝑜𝑜𝑑 𝑖𝑛𝑔𝑒𝑠𝑡𝑒𝑑 100 ×𝑑𝑢𝑟𝑎𝑡𝑖𝑜𝑛 𝑜𝑓 𝑝ℎ𝑎𝑠𝑒 1 (23)
)

was estimated with the following equation:

with the amount of food ingested = 6g, and phase 1 being the chewing phase before the first swallow. The percentage of incorporated saliva (relative to the product) was calculated as follows, based on moisture content of the initial cheese and of the bolus:

𝑝𝑒𝑟𝑐𝑒𝑛𝑡𝑎𝑔𝑒 𝑜𝑓 𝑖𝑛𝑐𝑜𝑟𝑝𝑜𝑟𝑎𝑡𝑒𝑑 𝑠𝑎𝑙𝑖𝑣𝑎 = 𝐻𝑀 𝑏𝑜𝑙𝑢𝑠 100-𝐻𝑀 𝑏𝑜𝑙𝑢𝑠 × (100 -𝐻𝑀 𝑐ℎ𝑒𝑒𝑠𝑒 ) -𝐻𝑀 𝑐ℎ𝑒𝑒𝑠𝑒 (24)
The percentage of incorporated saliva and the duration of phase 1 were determined during two separate experimental sessions, which can induce biased calculation (the salivary flow rate is known to vary daily and with the physiological state of the individual). Therefore, the average rate of saliva incorporation in the bolus during food consumption calculated with Eq. 23 was

only indicative and we preferred to use this parameter as a degree of freedom of the model for simulations.

In this study, neither the frequency nor the displacement amplitude of the soft palate during food consumption could be experimentally measured in a non-invasive way. Due to the uncertainty on the actual values of these parameters and to the fact that they have a similar effect on the air flow rate (Eq. 5), the volume variation of the oral cavity during the consumption of cheese was arbitrarily fixed to a reference value of 20% of the mean volume and the opening frequency of the velopharynx was a degree of freedom of the model.

The contact area between the liquid bolus and the product just before swallowing (A OLPdeglu

A 𝑂𝐿𝑃𝑑𝑒𝑔𝑙𝑢 = 3×V OPini R (25)
)

was calculated with the assumption of spherical particles:

The particle radius R was estimated with a compression device designed by the Laboratoire de Rhéologie (Grenoble, FRANCE) [START_REF] Yven | Consequences of individuals chewing strategies on bolus rheological properties at the swalllowing treshold[END_REF].

Air/bolus partition and mass transfer coefficients of aroma compounds in cheese bolus

The modification of the air/bolus partition coefficients of ethyl propanoate and 2-nonanone, due to modification of the bolus composition, was determined by the static phase ratio variation method (PRV) [START_REF] Ettre | Determination of gas-liquid partition coefficients by automatic equilibrium headspace-gas chromatography utilizing the phase ratio variation method[END_REF]. The mass transfer coefficients of aroma compounds into the bolus were obtained by non-linear regression from dynamic headspace experiments performed with Proton Transfer Reaction-Mass Spectrometry (PTR-MS) [START_REF] Doyennette | Mechanisms explaining the role of viscosity and post-deglutitive pharyngeal residue on in vivo aroma release: A combined experimental and modeling study[END_REF]. From those in vitro measurements, the relationships describing the modifications of air/bolus partition and the mass transfer coefficients of ethyl propanoate and 2-nonanone in bolus as a function of cheese/saliva mass ratio were integrated into the model. The cheese/saliva mass ratio r cs

𝑟 𝑐𝑠 = 𝑉 𝑂𝑃𝐷 𝑉 𝑂𝑃𝐷 +𝑉 𝑂𝑆 (26)
used in the equations is given by:

Estimation of cheese dissolution rate in artificial saliva

Considering the previously mentioned assumption of product dissolution, the dissolution rate was estimated from measurements of the release of ionic species into warm artificial saliva, using a conductivity probe. This measurement was chosen for its simplicity and for the possibility to perform on-line measurements. The release of salts is considered here as a marker of the matrix dissolution. The conductivity probe was calibrated at 35°C, with aqueous NaCl solutions prepared with deionised water. Concentrations were measured in g/L NaCl equivalent. A 6g-cylinder of cheese matrix, cut with a 24-mm punch, was placed in a beaker containing 400 ml of artificial saliva (including 0.216 g of mucin (from porcine stomach type II, SIGMA-ALDRICH) per 100g) and a magnetic stirrer. The monitoring of salt release (actually all species contributing to the conductivity signal) was made during two minutes using a conductivity meter (MPC HEITOLAB 350) and a probe (DCP 4). Product dissolution rate was determined by nonlinear regression using Matlab 7 (Natick, MA), using Eq. 27.

𝜙 𝑠𝑎𝑙𝑡 (𝑡) = 𝑣 × 𝐴(𝑡) × (𝐶 𝑃 𝑠𝑎𝑙𝑡 -𝐶 𝑆 𝑠𝑎𝑙𝑡 (𝑡)) (27)
where φ salt is the salt mass flux between the product and saliva, A(t) the product/saliva contact area and C Psalt and C Ssalt

Measurements of aroma compound retention by lubricated mucosa using Proton Transfer

Reaction-Mass Spectrometry (PTR-MS) the salt concentrations in the solid product and in artificial saliva respectively.

Flavoured solutions were prepared with mineral water (Evian), ethyl propanoate and 2nonanone (final concentration of 4.8mg/L for each). 134.8cm 3 Five panelists were recruited for this specific measurement and were instructed to suck up a mouthful of vial headspace with a straw, to swallow it and to continue to breathe normally through the nose. A minimum of three replicates was performed for each individual. During this experiment, the dynamic release of aroma compounds was measured online using Proton Rankings of samples were then obtained by the non-parametric test of multiple comparisons (Bonferroni method, level of significance set at 5%). For the correlation test, the Spearman test was performed, with a level of significance set at 5% (www.Anastats.fr

Results and Discussion

). Then, aroma compounds in the pharyngeal deposit are gradually released to the air phase.

During secondary swallowing events, the deposit is renewed by the food bolus coming from the oral cavity. However, the new layer is already depleted in aroma compounds, which leads decrease to zero.

4.2.Comparison of model predictions with experimental data for ethyl propanoate

As no interaction with lubricated mucosa was assumed for ethyl propanoate, simulations were performed using the model without lubricated mucosa compartment. To compare model simulations with experimental in vivo release data for ethyl propanoate, the air flow rate resulting from the frequency of velopharynx opening was adjusted to fit the simulated release kinetics before the first swallowing. The rate of saliva incorporation into the food bolus was also adjusted so that the simulated release kinetics after the first deglutition fitted the decay phase of the experimental curve.

All simulations have been satisfactorily fitted to experimental data (the average error of the model, for all cheeses and all individuals, was only 5.88%±2.81). The two unknown parameters of the model (the average rate of saliva incorporation into the bolus and the frequency of velopharyngeal opening) could be estimated.

Figure 4 presents some comparisons of simulations and experimental data for panelists S001 and S101. We can observe that these two individuals present diversified consumption and aroma release behaviours: for example, panellist S001 exhibits a consistently longer mastication time before the first swallow than panellist S101, and both panellists have longer mastication times for firm matrices, as one could expect. In spite of this diversity, the model correctly described the in vivo release profiles. Some slight deviations can be noticed, as for instance on figures 4.a and 4.b at about 10s and 30s respectively. A possible explanation could be that individuals have more complex behaviour during food consumption than currently assumed in the model; for example, the amplitude of jaw movement related to

parameter ΔV OA
The mean rate of saliva incorporation into the bolus Q in Eq. 5 may vary in the course of the mastication process. This cannot be correctly represented by the model with constant values of physiological parameters. One should also bear in mind that panellists' behaviour is not fully reproducible and the greyed area in Figure 4 represents an interval of ±1 standard deviation of 3 replicate experiments, which is statistically expected to contain less than 58% of the data.

OS (for all cheeses and all panelists), experimentally determined using the percentage of saliva incorporated in the bolus and the average chewing time, is 7.88±4.39 mL/min (figure 5, bar B). Its value remains similar to the one of Q OS Overall, the average rate of saliva incorporation into the bolus, either measured directly or estimated by fitting simulations to experimental data, was statistically much higher than the salivary flow measured with mechanical stimulation (9.00±3.65 mL/min vs. 2.65±0.87 estimated by fitting simulations to experimental data (9.00 ± 3.65 mL/min, bar C in figure 5). These values are also comparable (mean and standard deviation) to the salivary flow values measured experimentally during food consumption by other authors [START_REF] Mioche | The relationship between chewing activity and food bolus properties obtained from different meat textures[END_REF][START_REF] Patarin | Mechanical characterization of cheese food bolus: a new device[END_REF][START_REF] Anderson | The effects of bolus hardness on masticatory kinematics[END_REF].

mL/min, bar A in Figure 5) as already observed by Gaviao et al. [START_REF] Gaviao | Chewing behavior and salivary secretion[END_REF] (7.82±4.53 mL/min vs.

1.40±0.67 mL/min), meaning that the salivary flow obtained with a mechanical stimulus is not representative of the actual one during product consumption.

When products are compared (table 1), it appears that the average rate of saliva incorporation into the bolus for the cheese matrix Fh (firm with a high-fat content) is statistically lower than for the others cheese matrices (2-fold factor between cheese matrix Fh and cheese matrices Sh and FFl), probably due to a longest consumption duration of cheese matrix Fh, as already observed in previous studies [START_REF] Repoux | Combined effect of cheese characteristics and food oral processing on in vivo aroma release[END_REF]. Indeed, the rate of incorporation of saliva into the bolus can be expressed as:

𝑄 𝑂𝑆 = 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑣𝑜𝑙𝑢𝑚𝑒 𝑜𝑓 𝑠𝑎𝑙𝑖𝑣𝑎 𝑑𝑢𝑟𝑎𝑡𝑖𝑜𝑛 𝑜𝑓 𝑐𝑜𝑛𝑠𝑢𝑚𝑝𝑡𝑖𝑜𝑛 𝑏𝑒𝑓𝑜𝑟𝑒 𝑠𝑤𝑎𝑙𝑙𝑜𝑤𝑖𝑛𝑔
+ 𝑠𝑎𝑙𝑖𝑣𝑎𝑟𝑦 𝑓𝑙𝑜𝑤 𝑟𝑎𝑡𝑒, with the initial volume of saliva in mouth being the same whatever the consumption duration. If we assume that the salivary flow rate is constant during time, it follows that the longer the consumption duration, the smaller Q OS 4.3. Comparison of model predictions with experimental data for 2-nonanone will be. The presence of fat can also induce a decrease in salivary flow, as already observed during the consumption of model cheese matrices [START_REF] Tarrega | Aroma release and chewing activity during eating different model cheeses[END_REF]. In the present study, it is likely that the low value of the average rate of saliva incorporation into the bolus for the cheese matrix Fh is the result of the two phenomena mentioned above.

In vivo observations highlighted that ethyl propanoate and 2-nonanone exhibit quite different behavior in terms of persistence (figure 6): the release amount for 2-nonanone during phase 2 was systematically higher for all individuals and all cheese matrices than for ethyl propanaote (significant difference between the areas under curve of the two aroma compounds, Friedman and Bonferroni tests, p<0.05, not shown). We assumed that this phenomenon was the result of an interaction between this molecule and lubricated mucosa. A specific experiment, based on swallowing of aromatized air, was performed to check this assumption (see Material and Method section). An example of standardized release curve is shown in figure 7 for one panelist. Statistical analysis confirmed that 2-nonanone (solid line) persists much longer than ethyl propanoate (dotted line) in the exhaled air. This was observed for all panelists and all replicates. Measurements by PTR-MS of the in vitro release of these two molecules by diluting the headspace of a vial containing the same aqueous flavoured solution excluded a retention effect of the PTR-MS transfer line. All these results confirmed the existence of a significant retention of the 2-nonanone by oral, pharyngeal and/or nasal lubricated mucosa.

This phenomenon has already been highlighted for other volatile compounds such as ethanol or menthol [START_REF] Normand | Modeling the Kinetics of Flavour Release during Drinking[END_REF][START_REF] Van Ruth | Volatile flavour retention in food technology and during consumption: Juice and custard examples[END_REF].

Simulations performed with the model including lubricated mucosa compartments were compared to in vivo release data for 2-nonanone for cheese matrices FFl and Sl. No data was available for the two others cheese matrices since 2-nonanone was not detected because of its high hydrophobicity and of retention by fat. The frequency of velopharynx opening and the rate of saliva incorporation into the food bolus which were previously found for ethyl propanoate were kept identical. Results showed that all simulations satisfactorily fitted experimental data if a mucosa layer thickness of 100 µm and a transfer coefficient of 1×10 -1 4.4. Velopharynx opening cm/s were fixed for all individuals (figure 8). The average error of the model was 13.16%±6.80, for the two cheese matrices tested and for the ten individuals. These results will have to be refined in the future, through dedicated experiments to measure aroma compound retention by lubricated mucosa.

The results of the model fitting suggested that velopharynx opening can be different between individuals. We found that in 53% of cases, the velopharynx opening was synchronized with chewing frequency; in 37% of cases, it was synchronized with respiratory frequency and in the 10% remaining cases, subjects had an intermediate behavior. Although little data is available on this topic in the literature, these results are consistent with observations from Matsuo et al. [START_REF] Matsuo | Effects of Respiration on Soft Palate Movement in Feeding[END_REF], who studied the velopharynx opening of 9 persons by videofluorography during the consumption of various foods (banana, cookie, meat). They found that the opening frequency varied widely among subjects, but that the opening time was related to jaw movements. A link between breathing and velopharynx movements was also highlighted.

Analysis of model sensitivity

An analysis of model sensitivity was performed to evaluate the effect of some physiological or physicochemical parameters on the variation of nasal concentration with time. A standard nasal concentration curve was determined using mean values for panellist characteristics.

Then, each parameter was multiplied or divided by a factor of two (which is representative of the typical variability of physiological parameters).

To assess the effect of each parameter on release intensity, each simulated kinetic was scaled by the maximum of the standard nasal concentration curve. Three parameters were found to have a strong positive effect on release intensity: the product dissolution rate in the mouth, the mass transfer coefficient of aroma compound in the bolus and the air-bolus contact area in the mouth. The respiratory frequency had a significant negative effect on nasal aroma concentration by increasing aroma compound removal.

To assess the effect of each parameter on release dynamics, independently from global intensity, each simulated kinetic was scaled by its own maximum intensity. Results showed [START_REF] Gaviao | Chewing behavior and salivary secretion[END_REF]), E) Average rate of saliva incorporation in the bolus reported in [START_REF] Salles | In-mouth mechanisms leading to flavor release and perception[END_REF]) and F) Average rate of saliva incorporation in the bolus reported in [START_REF] Watanabe | A Comparison of the Effects of Tasting and Chewing Foods On the Flow-rate of Whole Saliva In Man[END_REF]). Table 1. Comparison of the average rates of saliva incorporation into the bolus (in mL/min) for the 4 cheese matrices. The letters indicate the classification group of the products (Friedman and Bonferroni tests are significant at the level 5%).

Table captions

4. 1 .Figure 3

 13 Figure 3 shows the time variation of 9 out of the 15 model variables for the release of ethyl

  to a general decrease in C FL . The air in the pharynx (C FA ) receives aroma compounds from the pharyngeal deposit and from the air of the oral cavity at each velopharynx opening (minor peaks observed on the curve) (figure 3.h). Overall, the kinetics of concentration in the pharynx C FA follow the one in oral cavity C OA , up to a dilution factor. The concentration of aroma compounds in the air in the nasal cavity C NA evolves similarly to the one of the air in the pharynx C FA (figure3.i). However, each inspiration of fresh air, which is aroma-free, makes the signal of the concentration in the nasal cavity C NA As illustrated in figure3, modelling provided insight in the change in non measured variables.An important outcome was an improved understanding of mechanisms governing aroma release kinetics and of the underlying assumptions.

of saliva in the bolus 3 Figure 1 .

 31 Figure 1. Schematic representation of the chronological steps of the consumption of a solid

Figure 2 .

 2 Figure 2. Schematic representation of the interconnected compartments and of the

Figure 3 .

 3 Figure 3. Time variation of 9 model variables for the release of aroma compounds during the

Figure 4 .

 4 Figure 4. Comparison of simulated and experimental release profiles for ethyl propanoate for

Figure 5 .

 5 Figure 5. Comparison of the average rate of saliva incorporation into the bolus estimated by

Figure 6 .

 6 Figure 6. Comparison of in vivo release profiles of ethyl propanoate (thin envelope curve and

Figure 7 .

 7 Figure 7. Comparison of in vivo release of ethyl propanoate (dotted line) and 2-nonanone

Figure 8 .

 8 Figure 8. Comparison of simulated and experimental release profiles for 2-nonanone for

  Figure 1.

Figure 2 .

 2 Figure 2.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

  Figure 5.

  Air was sampled from the subject's nose at a flow rate of 35 cm 3 /min via two inlets of a stainless nosepiece placed in both of the assessor's nostrils. The PTR-MS instrument drift tube was thermally controlled (T drift =60°C) and operated at P drift =200 Pa with a voltage set of

	U drift =600 V. Measurements were performed with the multiple ion detection mode on specific
	masses with a dwell time per mass of 50 ms. Ethyl propanoate and 2-nonanone were
	respectively monitored at m/z 103 (molecular ion) and m/z 143 (molecular ion). In addition,
	m/z 59 (acetone) was monitored as a breath marker and m/z 21 (signal for H 3	18 O + ) and m/z 37
	(signal for water clusters H 2 O-H 3 O + ) were monitored to check the instrument performances
	and cluster ion formation.
	3.7. Statistical analysis

Transfer Reaction Mass Spectrometry (PTR-MS, Ionicon Analytik, Innsbruck, Austria). The PTR-MS inlet was connected to the subject's nose via a 1/16" PEEK tube maintained at 60°C. Schott vials, equipped with caps fitted with two valves (Interchim SCAT or InterchimOmnifit, England), were half-filled with flavoured solutions and left at 35°C for 2 h to allow thermal equilibrium.

Due to the small size and the non-normal characteristic of the dataset, classical analysis of variance was not appropriate (univariate procedure using SAS/Stat® software). Instead, Friedman tests were performed with an Excel program (available at www.Anastats.fr).

Table 1 .

 1

	5. Table1
	Rate of saliva incorporation	FFl	Sh	Sl	Fh
	Mean value (mL/min)	10.50	10.80	9.43	5.27
	Standard deviation (mL/min)	5.70	4.97	6.16	8.08
	Group	A	A	AB	B

Acknowledgements