Modeling in food across the scales: towards a universal mass transfer simulator of small molecules in food

Multiscale modeling in food is the cutting-edge strategy to revisit food structure and food composition to meet specific targets such as bioavailability, oral perception, or to evaluate the contamination of food by chemicals. A special implementation of Langevin dynamics is proposed to describe mass transfer in structured food. The concepts of random walks over discrete times and physicochemical interactions are connected via an exact solution of the Fokker–Planck equation across interfaces. The methodology is illustrated on the calculation of effective diffusivities of small solutes in emulsions in relationship with their polydispersity, the volume fraction of dispersed phase d = [0.1, 0.4], the ratio of diffusion coefficients between the two phases, rD = [10−2, 102], and the partition coefficients between the continuous and disperse phases, K = [10−2, + ∞[. Simulated diffusion paths are detailed in 2D emulsions and the effective diffusivities compared with the core–shell model of Kalnin and Kotomin (J Phys A Math Gen 31(35):7227–7234, 1998). The same effects are finally tabulated for 3D emulsions covering the full range of food applications. The methodology is comprehensive enough to enable various extensions such as chemisorption, adsorption in the surfactant layer, local flows, flocculation/creaming.


Introduction
Mass transfer of small molecules plays a significant role in evaluating the shelf-life, the flavor, the mouthfeel, and the intestinal absorption of food constituents. The presented approach focuses on organic solutes with sensorial, nutritional, or sanitary interest, such as flavors, antioxidants or chemical contaminants originating from food contact materials. Similar descriptions could also apply to water vapor in dry foods. The state of the art offers a broad range of calculation and simulation methods to estimate diffusive fluxes in food [1] with a particular food geometry or to evaluate diffusivities [2][3][4][5][6] or thermodynamic properties [7][8][9][10][11][12] from molecular structures. There is, however, a lack of methodologies to describe the effect of food structure and the spatial organization of phases at intermediate scales above the molecular scale and below the scale where the effective medium approximation can be assumed. Although Fickian diffusion is strictly limited to ideal mixtures [13], continuous approximations averaged over a representative volume are physically acceptable at the cost of proper identification of driving forces and transport coefficients [14]. Such descriptions are commonly used to simulate mass transfer of water and organic solutes in complex food products or materials [15][16][17][18][19], but they fail to capture the effects of temperature (thermodynamic effects), the evolution of the microstructure during processing, ingestion or digestion, and more generally for any fractal arrangement of phases. As examples, they are poorly adapted when the dispersed phase is colloidal and is moving to the continuous phase by Brownian motion, sedimentation or aggregation, or when the mean free path of molecules is comparable to the length scale of the details, or when the slowdown of random walks is controlled by rare traps [20]. Particle methods and other meshless methods offer an alternative by avoiding space discretization and Lagrangian formulation. For fluids, smoothed particle hydrodynamics [21][22][23][24] and moving particle semi-implicit scheme [25] are the most popular. Alternative techniques of integration of kinematic equations on intricate spatial scales comprise high order spectral approximations [26], lattice gas models [27] and kinetic Monte-Carlo [28,29]. Describing mass transfer of solutes instead of advection requires, however, defining particles commensurable to molecules. They carry properties of real molecules, including mass, charge, position, and velocity. Representing all molecules of a food portion will remain inaccessible for several decades, but efficient and accurate descriptions can be recovered by following the trajectories of a smaller collection of solutes governed by Langevin dynamics (LD). The approach is very generic, and this study emphasizes on the correct implementation of thermodynamical concepts at interfaces such as local thermodynamical equilibrium, detailed balance, time and microscopic reversibility. A distinction between all concepts has been reviewed in Gorban [30]. In LD simulations, interactions with surrounding molecules are captured via an implicit solvent model and a proper stochastic equation, enabling the incorporation of temperature effects.
In an infinite medium, the representation is equivalent to the Fokker-Planck equation for the evolution of the probability density of a population of molecules and hence similar to the convection-diffusion equation. No bias is therefore expected. This work illustrates how such a general methodology can be devised and implemented. With this respect and without a loss of generality, a simple casestudy is presented. It aims at calculating (in few minutes or hours on a multicore processor) the effective solute diffusivities in 2D and 3D emulsions near their percolation threshold. The chosen problem captures enough complexity-the coupling between transport and thermodynamics around complex interfaces while enabling visualization of 2D trajectories. The formulation shows how properties broadly available in the literature, such as intrinsic diffusion coefficients, activity or partition coefficients, association/dissociation coefficients with macromolecules, can be mapped onto LD simulations. There is additionally significant literature on the subject to be used for validation and discussion. Earlier results, including the Maxwell-Garnett mixing rule, have been recapitulated by Crank [31]. Subsequent works [32,33] considered obstacle effects (e.g., impenetrable globules), but neglected the partitioning of the diffusing species between the continuous phase and the globules. Ronis [34] studied the effect of partitioning for a continuous phase separating two globules and evidenced the role of trapping in globules on effective diffusivities. The repeated ring approximation failed, however, to provide a practical solution for solutes with a low chemical affinity for the continuous phase. The core-shell (or coated sphere) model of Kalnin and Kotomin [35], discussed further in Kalnin et al. [36], provided the first extension of the Maxwell-Garnett equation when partitioning occurs. The free geometry parameter, Φ , introduced in Eq. (29) in Kalnin et al. [36] has however no physical meaning and needs to be adjusted according to the arrangement and morphology of the dispersion.
LD simulations solve all previous issues by enabling an averaging over various arrangements and packaging of globules. The principle of simulations and the theoretical treatment of mass transfer across interfaces are developed in section two. In particular, this work generalizes the kinetic approximation of Lejay and Maire [37] at interfaces, by introducing a proper treatment of sorption isotherms. The methodology to determine the effective diffusivities from long-term dynamics are detailed in section three. The following section recapitulates the scaled results obtained for a wide range of conditions, including the effect of diffusion coefficients ratio, partition coefficients, volume fraction of the dispersed phase, and emulsion morphology. The reliability of our results is finally compared to the analytical model of Kalnin and Kotomin [35].

Langevin dynamics (LD) formulation in structured foods
Duderstadt and Martin [38], transport theory in condensed matter differs from the approaches encountered in classical physics because it involves particles, not a continuum theory of matter, such as fluid dynamics. In particle transport, the random nature of particle interaction events requires replacing densities and mass fluxes by a field of probability densities or distribution functions. The exact number of particles in a specific region of space and time cannot be predicted. Concentrations and particle densities used at macroscopic scale represent, indeed, only the average number of particles per elemental volume unit within a certain period of time. Though it is broadly used at macroscopic scale, the first moments (averages of particle densities or fluxes) are not sufficient to describe diluted systems at microscopic scale, and two classes of problem, so-called direct or indirect problems, address the issue. The direct problem consists in calculating the distribution of particles in the medium given the composition and geometry of the host medium and the location of eventual sources of particles in the medium. The inverse problem is more difficult and aims at evaluating the effective characteristics of the medium through which the particles have propagated. The underlying random walk process has been described either as a succession of random steps on discrete positions or times or as a continuous process in time or space. Continuous representations are preferred to describe anomalous diffusion. In this work, solute transport is described as a special case of skewed Brownian motion random walks. Following the same ideas as Walsh [39] and Yor [40], the diffusion process is generalized by replacing random variates of cartesian coordinates by polar ones, where each particle moves along rays emanating from its previous position. This scheme is well adapted to describe random walks with discontinuous diffusion coefficients [41]. In short, particle positions are continuous, but they evolve over discrete times.

Equations of motion
Solutes dispersed at low concentration are described as self-moving particles due to thermal fluctuations and without interacting with themselves. This transport mechanism homogenizing solute density in any homogeneous and condensed phase is known as tracer diffusion or molecular diffusion. At macroscopic scale, the distribution of the diffusing species obeys to the general Fokker-Planck equation [42]. At microscopic scale, the motions of small solutes, e.g. not massive particles, are essentially controlled by viscous forces, so that their positions result from the random collisions with surrounding molecules. The energy of the motion is randomly pumped in or pumped out the system via a stochastic formulation of equations of motions. By assuming the conservation of momentum and angular momentum, the n-dimensional stochastic trajectory is captured by a special case of Langevin equation: where r(t) is the position vector in an Eulerian reference frame. A is a n × n matrix associated with low-frequency deterministic motions, whose eigenvalues i i=1..n represent the corresponding reciprocal correlation times. ξ(t) is a random vector, whose components are normal variates, with zero means and unitary standard deviations. The last term on the right-hand side is related to the symmetric positive-definite diffusion tensor, D, via its unique Cholesky decomposition: where Λ is a diagonal matrix of eigenvalues and where the columns of B are the corresponding unit eigenvectors. In practice, Eq. (1) is straightforwardly integrated with a fixed time step, t , via an Itô-Taylor expansion [43]: where I is the identity matrix.
When δt is chosen much smaller than 1/α 1 , where α 1 is the largest eigenvalue of A, the visited distance by molecular diffusion is much larger than the distance visited by any other deterministic cause. t (t) is the random increment of the Wiener process associated to trace diffusion. Since all increments increase mainly as √ t the integration scheme (3) has a low global order of 1/2. The convergence is slow and longterm iterations are required to sample the trajectory in space and frequency. Since only the second moments of the stochastic trajectory are desired, the choice of the initial starting positions do not modify the final result. Hence a parallel integration of Eq. (1) over m independent trajectories including each n steps offer considerable acceleration by several magnitude orders.
In this work, Eq. (1) applies to both the continuous phase, denoted "c", and the dispersed phase in it, denoted "d". The corresponding diffusion tensors are noted D c and D d , respectively. For simplicity, diffusivities were considered isotropic and diagonal matrices replaced diffusion tensors with uniform elements, D c and D d with D c and D d representing isotropic scalar diffusivities with SI units in m 2 s −1 .

An intuitive description of the free energy landscape on diffusivities
Diffusion on a free-energy landscape is not intuitive but plays a significant role in dense and multiphasic media such as food. It introduces two important consequences on trajectories: it can make some positions more probable and it modifies waiting times before a translation can be observed (i.e., a change in macroscopic state). For instance, heterogeneous dynamics are responsible for the broad mass dependence of solute diffusivities observed in solid polymers [3]. A simplified interpretation of the role of free energy landscape on diffusivities is presented in this section. The random walk along a periodic two-well potential involving two macrostates (or positions) denoted "c" (continuous phase) and "d" (dispersed phase) separated by a distance l∕2 is presented in Fig. 1, and corresponds to the In this 1D case and the transport across interfaces, it is essential to note that the waiting times are random not the positions. The thermodynamical equilibrium between both states is determined by the partition coefficient between the two macrostates, K , defined as: where p k eq k=c,d and c ex k=c,d are respectively the solute probability density at equilibrium and its excess chemical potentials in phase k . A more general definition of the local thermodynamical equilibrium at the interface between the macrostates/positions c and d can be envisioned within the framework of transition state theory (TST) modified by Kramers [44] and Mel'Nikov [45] and experimentally evidenced for transport of isolated molecules by Vestergaard et al. [46]. According to Kramers and Random walk along a periodic double-well potential in one dimension corresponding to Eqs. (6) and (7) in the limits of weak friction (i.e. no inertia), the escape rate of a Brownian particle out of a deep potential well is subjected to many recrossings. Additional discussions on symmetrizable (reversible) jump processes can be found in Chen [47]. At the exact position of the interface Ω , the probabilities to move forward or backward are the same. In other words, the flux of probabilities of transitions from c to d and in the opposite direction should be equal: The trajectories of solutes at the interface are not ballistic and are induced by the thermal fluctuations of the total free energy of the entire system at constant volume, F . A translation (whatever the direction: from A to B or from B to A) occurs when the system acquires an energy larger than the critical value F ‡ separating the two wells. The height of this energy barrier is greater than the free energies associated with the sorption states c and d, denoted F k k=A,B . By considering that the thermal fluctuations are Poissonian events (independent and stationary), the escape rates from A to B are hence exponentially distributed with an expectation, in the canonical ensemble, given by: where h is the Plank constant and k B the Boltzmann constant.
By noticing that the Wiener increments are identical for the sequences d → c → d and c → d → c, the effective diffusivity on the double-well surface energy reads: min p c eq , p d eq . Equation (7) is illustrative of a pure thermodynamic control on diffusivities, which should not be neglected in simulation and require proper implementation. Effective diffusivities, D, deviate from the homogeneous medium as soon as large K values are met. More discussions on the effect of surface energies can be found in Vitrac and Hayert [29]. In food systems, partitioning between phases is very likely due to specific interactions with the host medium at molecular scale [10,11] or due to the difference in size between the host molecules in phases c and d [8,12].

Interface conditions for discrete times: partition function for a single hopping event
Equation (7) demonstrates that that changes in chemical affinity near interfaces increase correlation times (see the interpretation of the thermodynamic factor in Eq. 10 of Fang and Vitrac [4]). Integrating Eq. (1) over discrete times while enforcing local thermal equilibrium via Eq. (5) and variable residence times is particularly challenging. Classical Riemann-Stieltjes integral of microscopic displacements along time cannot be used across interfaces, because almost all paths of Brownian motion fail to have bounded variation. Itô calculus already applied to derive Eq. (3), provides an alternative by defining the trajectory as a limit in probability of Riemann sums, but the probability density function must be twice continuously differentiable in space. This requirement is again not verified close to any interface ∂Ω, where abrupt changes in diffusion coefficients or chemical affinity occur. To reach a uniform probability at equilibrium, particle tracking methods used in porous or fractured media (i.e. D c ∕ D d > 1 and K = 1) rely on semi-reflecting barriers as reviewed by Delay et al. [48] and Salamon et al. [49]. Reflection is considered, however, only on the side of the high dispersive medium [50] with a r e f l e c t i o n r a t e . One corollarily is that a different number of particles cross the interface in a given amount of time [51] from c to d, and from d to c in violation of the principles of local thermodynamical equilibrium. The verification at any time and in a probabilistic sense of mass conservation and of a local partitioning of particles, as envisioned by Kramers, is challenging and not subjected to abundant literature.
First-passage properties offer an alternative to deriving a thermodynamically consistent numerical scheme at the microscopic scale (i.e. at both microscopic space and time scales). First-passage and first-escape times are defined as the time to hit the interface and to leave the current phase, respectively. By assuming that the local direction of translation is not affected by the proximity with interfaces (i.e. no long-range interactions such as Coulombic ones), the probability to hit the interface within a given time is given by the distribution of the particle along the line connecting the current position (coordinate along the line x = x 0 , at t = t 0 ) and the interface Ω ( x = 0 ). By indexing 1 the original half-plane x ≥ 0 and 2, its complementary part on the other side of the interface, the probability density function to find a particle at position x along the line after the dimensionless time Fo = D 1 (t−t0) is given for the two half-segments on both sides of Ω as: In domain 1, the particle may either diffuse or being reflected. Domain 2 contains the fraction of transmitted particle (i.e. crossing Ω ). Constants A , B and C are inferred from the condition of local equilibrium Fo , the distribution of probability reads piecewisely across the interface: Equation (9) enables to calculate the partition function associated with the considered hopping event. By interpreting it as the ratio between the fraction of remaining and transmitted solutes for a particle initially located at x 0 from Ω after a dimensionless duration Fo , one gets: On finite time scales, remaining solutes (i.e. in a probabilistic sense as we discuss only one event) include solutes which did not have "enough" time to cross and reflected ones. Both populations are indiscernible. On very long time scales, all remaining solutes have been already transmitted at least once and were back. The previous result of TST (Eq. 5) is now extended to semi-infinite geometries (i.e. the system is aperiodic). The new partition function between the remaining (i.e. reflected) and crossing probabilities depends not only on K but also on the square root ratio of diffusion coefficients (i.e. ratio of transmission rates) between the two half-spaces separated by the boundary Ω . The thermodynamical equilibrium is exactly and only verified at the interface Ω (i.e. the ratio of probabilities is set by K ), but it cannot be reached macroscopically since both media are semi-infinite. When particles crossed media with different properties (diffusivity or chemical affinity), they are not only dispersed, but they are also refracted and even diffracted when they cross multiple interfaces.

Choice of time steps and generalization to 2D and 3D structures
As discussed by Balescu [52, see Chapter 9] the sampling of random walks in inhomegenous media is significantly improved (i.e., better ergodicity) when ordinary dynamics such as Eq. (1) are transformed into statistical dynamics. In this perspective, the concept of "sampled time" should be envisioned as a "time-ordered path" to explore the configuration space accessible to particles rather than the true physical time (see Chapter 4 in Landau and Binder [53]). A not necessary but sufficient condition to ensure a Gibbsian equilibrium state is to verify at any time the detailed balance condition, also known as microreversibility, across ∂Ω [54]. The term of Gibbsian equilibrium is used because statistical mechanics is associated since its foundation with two partly incompatible definitions of equilibrium. Chemical and food engineers are familiar with the Gibbsian equilibrium. The equilibrium is defined from the stationary distribution of ensembles maximizing entropy. The Boltzmannian equilibrium imposes only the macro-variables to be almost constant on the equilibrium region, such as considered previously for semi-infinite media (see for example the discussion in Werndl and Frigg [55]). In shorts, Eq. (10) applied to each "jumping event" along a line in n = 1, 2, 3 dimensions enforces Gibbsian equilibrium irrespectively the implementation of the interface condition. Since the proposed algorithm of first passage near interfaces is analogous to the Fermat's principle in optics, any step of duration t intercepting the interface was split into two times t 1 + t 2 = t . t 1 was the ballistic time required to hit the interface and t 2 the complimentary time to complete partial specular reflection or partial refraction. Because only one possibility was considered during t 2 (particles were not duplicated), reflection was chosen when a random number picked uniformly between 0 and 1 was lower than pr 1 → 1| Ω, Fo, x = x 0 and refraction otherwise. Because refractive indices were associated with the phase velocity of the diffusion front, they were chosen proportionally to √ D 1 ∕ D 2 . Implementating reflection in a billard ball fashion (i.e. with equal indicent and reflected angles across the local normal) and reflection according to the Snell-Descartes law guaranteed that the normal probability flux was conserved across the interface for any initial incidence of the particle in 2D and 3D.
On time scales commensurable to t , any trajectory looks ballistic (jump process) with a random direction. Since an analytical equation was used, increasing t did not provoke any loss of information, while interfaces could be assumed plane and separating two semi-infinite halfspaces. The distribution of probability in a convex region is obtained by sampling several recrossings from different directions. For convex objects such as globules, the integration time step must be chosen small enough so that the probability of leaving the smallest object (e.g., globule) from its center and without making a turn is almost zero. The escape probability for small objects increases indeed dramatically with highly curved interfaces. Figure 2 compares exact analytical solutions and the corresponding escape rates for a single macroscopic jump between two infinite and semi-infinite spaces [Eq. (9) in 1D] and for spherical globules (numerical simulation). The comparable condition corresponds to Both descriptions succeed in reproducing the diffraction ring (two local maxima): one corresponding to the initial position and an other emanating from the outside of the interface where the chemical affinity and diffusion are higher. The escape rates of the half-space (one exit point corresponding to one single ballistic event in the simulation) and of the entire object (many exit points) do not coincide and should not be confused together (Fig. 2b). Only the domain associated to low Fo values is sought in simulations. Additionally, the calculations illustrate that diffusion has broken apparent symmetry around fractal and convex objects with poor penetrability. Mass diffusion evolves rapidly to surface diffusion around impermeable crowded or percolating objects. The complications associated with the arrangement of spheres was discussed early in Chapter 12 of Crank's reference textbook [31].

Food description
The presented case-study addresses generically the diffusion, sorption/desorption, and permeation in and across a food including a dispersed phase, as shown in Fig. 3. The continuous phase is assumed to be dense, viscous and associated with homogeneous diffusion properties. The discontinuous phase consists of globules with fixed positions. The eventual surfactant layer is neglected. Solutes (aroma, contaminants, etc.) can diffuse freely within both phases. The effective diffusivities are averaged on several 2D and 3D periodic systems sampling the real size distribution, arrangement, and packing of globules. A dimensionless formulation applicable to many emulsionned food products (dairy products, seasonings, etc.) is proposed. In its present form, application to foams is discarded, and all phases need to be dense. The fundamental reason is that the particle velocity should be comparable between all phases: few µm/min or less in liquids and solids against several km/min in gases.

Scaled emulsions
2D and 3D emulsions were represented as dispersed globules mapped onto a scaled periodic cell over Periodic conditions enabled to reproduce an infinite medium with a periodicity l 0 . The size of each globule was sampled from a log-normal distribution with arithmetic average, ∅ , and standard deviation, ∅ . Globules were inserted randomly according to a generic algorithm independent of the shape of the object to be inserted. The positions of all phases remained fixed during the simulation. To facilitate the sampling of the free space, a M × M or M × M × M fine grid (with M > 5 ⋅ 10 3 ) was used to mark sub-cell as occupied or free space. For each new globule of diameter ∅ to be inserted, a non-centered coarse grid with a step-size equal to 2 ⋅ ∅ was superposed to the fine one.  All sub-cells of the coarse grid that contained enough free space were randomly numbered and iteratively tested. Within each sub-cell, all possible centers laying within the relative positions −∅∕ 2, ∅∕ 2 × −∅∕ 2, ∅∕ 2 or −∅∕ 2, ∅∕ 2 × −∅∕ 2, ∅∕ 2 × −∅∕ 2, ∅∕ 2 (the subcell center was assumed to be located at 0,0 or 0,0,0) were checked for insertion so that the domain enclosed by the globule did not overlap any occupied site on the fine grid. When no feasible position was encountered on the coarse grid, the coarse grid was randomly translated up to 200 times. When no insertion was possible, a new globule was randomly picked. The whole process stopped when the set density in globules, defined by the fraction of total occupied sites, was achieved within a tolerance about 1%. All the structures that did not meet the target were discarded.
The tested conditions are summarized in Table 1. Two extreme morphologies were studied: a monodisperse emulsion, denoted E0, and a polydisperse one, denoted E1. The effect of tortuosity was tested by varying the volume fraction of dispersed phase, d, between 0.1 and 0.4. To reduce possible bias due to the arrangement of globules, all trajectories were averaged over 10-20 different structures. Examples of emulsions and their corresponding size distributions are given in Fig. 4. According to the set density in globules, each morphology included between 70 and 800 globules in 2D. Calculations in 3D comprised more than 10 4 globules.

Tested range of physico-chemical properties
The diffusion coefficient in the dispersed phase, D d , was chosen as reference as it enables to emphasize the effect of entropic trapping. Its value was set lower than l 2 0 4M 2 t −1 so that more than 95% of elemental trajectories during t had a length lower than the step of the fine grid. This condition ensured that even small globules were correctly sampled. The final granularity was far beyond the possibilities order spectral methods, such as finite elements. Simulations were performed for ratios r D = D c D d ranging between 10 −2 and 10 2 . The range covered most of the food applications (oil in water, water in oil, yogurt, sauce, etc.). Partition coefficients K , defined as the continuous-to-dispersed phase ratio of concentration at thermodynamical equilibrium were sampled over the whole scale: from 0 (no affinity for the continuous phase) to +∞ (no affinity for the dispersed phase). The interface conditions previously described was implemented using a search procedure, which did not presuppose any shape for the dispersed phase (e.g., disc or sphere). Interfaces and their surface normal were detected with the following tier approach: tier (1) coarse search on the fine grid, tier (2) continuous search by finding the intersection of trajectory with the parametric equation of the object contour, tier (3) recursive search of reflections using steps 1 and 2. An example of trajectory (one particle) is detailed in Fig. 5. The whole scheme can manage up to ten successive reflections and diffractions during t if needed. Simulations are performed at constant temperature.

Identification of diffusion coefficients
The effective diffusion tensor D was estimated from the asymptotic rate of increase of the covariance matrix of displacements with diffusion time: where ⟨⟩ t,i=1..m is the ensemble average operator over all simulated times, trajectories and sampled emulsions. It is worth noticing that t and are not equivalent. The first is indeed the simulation time (i.e., the physical time counted from the origin), whereas the latter is so-called diffusion time and represents a lag. Since only stationary trajectories vs are studied (i.e., disregarding the origin of times), it is a common practice to denote lags t in figures when it represents . The isotropic diffusivity is finally given by the expectation of the quadratic form: where tr() is the trace operator and n is the number of dimensions. msd( ) represents the average mean-square displacement (variance assuming zero mean) on time scale . The asymptotic regime is reached when the displacements are time-independent and with a white noise spectrum. In the absence of correlations, msd( ) increases as a power law , with an exponent close to unity. Confinement inside or between packed globules increases considerably the time required to reach independent displacements. For the particular value K = 0 , confinement occurs over all time scales, and the asymptotic diffusion coefficient is zero since lim →∞ msd( ) = dl 2 0 .
Simulation lengths were adapted to get msd values s times larger than dl 2 0 , with s ≥ 2 . According to Eqs. (11) and (13), the necessary condition reads: s values up to 20 and m values up to 500 were applied in this study. All simulations were performed for products mt longer than 10 7 δt. When m is large, starting from a distribution of particles reproducing an equilibrium distribution accelerates convergence.

Trajectories
The effect of r D and K on the trajectory of m = 5 diffusants is illustrated in Fig. 6 for r D = 2 and K = {0.1, 1, 10} for a polydisperse emulsion with d = 0.3 . The corresponding cumulated diffusion time was set to m ⋅ t = 5 ⋅ 10 6 ⋅ t to get a partly ergodic system but without reaching equilibrium (i.e., a single particle could not visit all globules; as a result, the trajectories can be distinguished). Corresponding residence times averaged over all diffusants were calculated on a 50 × 50 grid and plotted as isocontours.
Raising the affinity for the dispersed phase ( K → 0 ) confined the trajectories inside the globules (Fig. 6a) and increased the residence time in the globules (Fig. 6b). For K = 1 , the trajectories spread over all the simulated without concentrating in any particular phase (Fig. 6c). For all conditions and irrespectively of the K value, the residence times were longer inside the globules than in the continuous phase (Fig. 6d), in agreement with the twice lower diffusion coefficient in the disperse phase. In the very long term, residence times were not uniform between globules. As exemplified in Fig. 2b, c, dwelling times in large objects (slow escape rates) than in small ones (high escape rates).
The condition K = 10 filled more uniformly the space between globules (Fig. 6e) while visiting less frequently the dispersed phase (Fig. 6f ). In regions where globules were concentrated, residence time increased due to more tortuous paths.
From previous examples, trajectories exhibited strong correlations either in globules or between globules. Locally low diffusivities, tightly packed globules, and high chemical affinity for the dispersed phase caused confinement alone or in combination. Persistent entropic trapping associated with negative correlations is known to lead to significant subdiffusive regimes on intermediate time scales (i.e., before displacements become decorrelated).

Mean-square displacements
Dimensionless mean-square displacements (msd) versus dimensionless diffusion time ( tD d l −2 0 ) are plotted for a wide range of combinations r D × K × d × E in Fig. 7. Msd values of each curve were averaged over the trajectories of 100 diffusants. Since longterm diffusion can propagate only through the continuous phase, the asymptotic diffusion of the continuous phase, r D , is shown as reference. On a log-log scale, the asymptotic regime is identified by a unitary slope ( ≈ 1 ) and the effective diffusivity by its intercept. Negative correlations (subdiffusive regime associated with (t) ≤ 1 ) slowed down the increase of msd and translated the curve to the right.
In decreasing order of significance, r D and K had the main effects. Increasing r D translated msd curves in the same fashion (Fig. 7d-f ), whereas the effect of K (for K ≠ 1 ) was not monotonic and modified the local slope (t) on intermediate time scales (Fig. 7a-c). Large K values increased diffusion rates on all time scales. By contrast, the effect of small K values was not symmetric and affected more significantly long-term diffusivities with drops down to two magnitude orders in tested conditions.
In details, K and r D could have either cooperative or antagonist effects according to the considered condition. This interpretation was already suggested by Eq. (10). When the diffusion is slow in the continuous phase, the diffusion through the dispersed phase compensates the slowdown partly. This effect is very likely when the concentration in globules is high and when the values of K are low but not too small. By contrast, large K values tend to reroute all particles between and around globules as they are mainly seen as obstacles. Close to the percolation threshold, the dispersed phase becomes almost a co-continuous phase, and the roles of each phase can be reversed. When diffusion flux occurs mainly across the dispersed phase, the cross-section of contact and passage between globules limits the effective diffusion coefficient.

Effect of K
The complex effect of K on the asymptotic diffusion coefficient is analyzed in Fig. 8. A search algorithm was applied to identify the hydrodynamic regime while discarding the very terminal region with poor sampling. Confidence intervals were calculated from the repetitions for different emulsions (Table 1). Uncertainty was ranged between a few percents and a factor 2 according to tested conditions. The poorest estimates were obtained when low r D and K values were combined together. Their combination increased the duration of the sub-diffusive regime drastically and increased the total simulation cost dramatically: both simulation length and the number of particles needed to be augmented. In practice, the effect of very low r D and K values could not be explored with a reasonable cost beyond a practical lower limit of 10 −2 . However our methodology and code remain valid beyond this limit. It was verified that the asymptotic diffusion coefficient was 0 when K = 0 whatever the applied r D value. By contrast, there was no computational limitation for applying large K values. Infinite K values were simulated by assuming a value of 10 99 . This condition made the globules impenetrable and asymptotic diffusion coefficients D D d converged to r D .
For the same r D value, a variation of K between 10 −2 and +∞ generated a variation of D D d up to five folds. A greater effect was even observed for values lower than Fig. 6 a, c, e Typical trajectories during 10 6 ⋅ t , b, d, f corresponding residence times in the emulsion plotted in Fig. 4g  ( d = 0.3 ) for r D = 2 and a, b K = 0.1 , c, d K = 1 , e, f K = 10 . Five independent trajectories are plotted, residence times are averaged over the five trajectories trajectories and calculated on a 50 × 50 grid 10 −2 . The effect of K varied significantly according to the value of r D . When the diffusivity was higher in the discontinuous phase, the asymptotic diffusion coefficient did not vary monotonously with K (Fig. 8d) and reached a maximum close to K = 1 . The maximum was exactly located at K = 1 when r D = 1 (Fig. 8e). The maximum was obtained for K values slightly lower than unity for r D < 1 and above unity for r D > 1 (Fig. 8f-j). K → ∞ did not lead to any maximum.
It was noticeable that the effect of d was more pronounced for values of K and r D above unity. For r D values higher than 1, the maximum was obtained for the lowest concentration in globules and was always lower than the r D value. For r D values lower than 1, the trend was the opposite, the maximum was obtained for the highest d values and reached values higher than r D .

Effect of r D
Mean square displacement presented in Fig. 8 suggested that r D was the main scaling coefficient of D D d values. Figure 9 plots on a log-log scale the values of D D d versus r D to identify the value of the critical exponent which would correspond to the asymptotic model D D d ∝ r D . Besides, the plotted curve y = x makes it possible to identify the particular conditions, which verified D D d > r D .
The critical exponent was lower than 1 for almost all conditions except when r D and K values were together lower than 2. values slightly higher than 1 yielded effective diffusion coefficients higher than in the continuous phase only for intermediate K slightly lower than 1 (Fig. 9b-d). As a rule of thumb, decreased with r D . The decay varied from values close to 1 down to 0.25 in this study. The decaying rate was minimal for large K values. It was almost 0 when K → ∞.
The results simulated with a complete reflection against globules ( K → +∞ ) are compared with the well-known tortuous model D = (1 − d)D c [56]

Comparison of simulated results with the model proposed by Kalnin and Kotomin [35]
A more general attempt of validation was thought by comparing our simulations with the predictions of Eq. (29) in Kalnin and Kotomin [35]. In this model, diffusion occurs in concentric spheres/discs. The polar/spherical topology does not match the descriptions in Cartesian coordinates, but it offers an acceptable description when all globules are noncontact. In our notations, the analytical model in 2D reads: Equation (14) accepts a free parameter associated with the volume ratio of concentric regions, 0 ≤ Φ ≤ 1. It was fitted to our simulated results obtained for the same d value and morphology. Simulated and analytical estimates are compared in Fig. 11 for all tested conditions where a good convergence was achieved. The agreement was globally correct, but negative deviations were observed for inhomogeneous emulsions ( E 1 ) and positive ones for the smallest D D d values. Figure 11b provides an interpretation of the deviations by showing that Φ cannot be captured with the sole effect of d , but that it depends also non-linearly of K or r D . Confidence intervals on Φ were obtained indeed by adjusting independently Φ to the results simulated with similar K and r D values. The comparison of the width of confidence intervals with the effect of d shows that contribution of d is as high as the other variations due to K or r D . As a result, Φ should be envisioned as a nonlinear function of d , K , r D and of the morphology. The presented results extend therefore the model and the conclusions of Kalnin and Kotomin [35] on more realistic morphologies.

Diffusivities in 3D emulsions
Simulations on 3D emulsions generalize the results obtained on 2D ones. They are more expensive to calculate and require computers with multiple cores (up to 256 cores were used) to reach convergence in a couple of hours. Asymptotic diffusivities are plotted in Fig. 12 as isocontours. Contour plots exhibited a half square shape confirming the predominance effect of r D (rate of diffusion in the continuous phase) when K values were greater than unity 10 −1 . The effect of K was significant mainly for values lower than 1. Oscillations were observed for some isocontours and were related to the non-linear combination of the effects of r D and K around unity, as already shown in 2D emulsions (Fig. 8). Concentration in globules and polydispersity modified the location of the maximum effect of K scale at constant r D values.

Conclusions
Mass transfer in food is known to interact with the food structure. Effective properties are not only scarce, but they do not explicit the complex relationship with local composition and microstructure. Multiscale modeling of (14) mass transfer is with this respect much more preferable as it opens the possibility of optimizing properties and exploring virtually new trends regardless of the existence of any underlying theory and of the food feasibility. This work proposes a comprehensive framework to simulate mass transfer over discrete times using Langevin dynamics (LD) capable of mimicking the random walk of diluted solutes in polyphasic and multiconstituent food systems. The parameterization of the simulation is simplified by the use of intrinsic transport properties respectively to each phase, such as diffusion coefficients, activity coefficients or sorption isotherms. These properties, as well as detailed microstructures, are either broadly available or can be measured independently. At each time step, the continuous and probabilistic interpretation of the first "interface hit" enforce both mass conservation and local thermodynamical equilibrium across the picked jumping direction. The integration scheme can handle potentially any kind of interfaces or frontiers, including interfaces between two fluids, one fluid and a solid, two solids or interactions with macromolecules. The only condition is that the normal direction at the "hit" point is known. In addition, the analogy with the Kramers interpretation of chemical reaction paths opens the extension to diffusion-reaction problems and chemisorption. A detailed case-study is shown covering hundreds of configurations met in emulsionned food systems when the globules are diluted or concentrated near their percolation t h r e s h o l d : K = 10 −2 , ∞ × r D = 10 −2 , 10 2 ×d = [0.1, 0.4] × E 0 , E 1 . The use of a dimensionless formulation enabled to cover almost all practical cases far beyond the domain of applicability of the Maxwell-Garnett equation and its extensions. Beyond their analogy with microscopic observations of emulsions, trajectories in 2D systems offer of visual interpretation of the geometric construction of skewed trajectories in crowded environments. Transitions between free and confined diffusion under thermodynamical control are directly observable. To accelerate convergence, it is preferable to implement the methodology on multicore processors with small systems rather than simulating a single very large system. In its 2D and 3D versions, the simulated systems might be seen as oversimplified emulsions without surfactants and relative displacements of the globules respectively to the continuous phase. It is, however, straightforward to add a third phase to represent globules as core-shell objects and to implement proper interactions with the considered solute. Mass transport in micelles and liposomes can be represented in the same fashion. Adding the contribution of sedimentation/creaming, coalescence and Oswald ripening requires integrating the Stokesian dynamics of the globules. By replacing a series of static and randomly generated emulsion by a series snapshots representing the evolution of the emulsion over time, fundamental results proven on bond-percolation models [57][58][59] can be used to derive dynamic diffusivities in a straightforward manner. The time-frequency frequency effective diffusion coefficient, D dynamic shelf-life, bioavailability and food safety. Tackling the coupling of these models across scales is the next frontier for systems (reactions, structures) evolving on different time scales. We will present the coupling between transport and oxidation reactions in a companion paper.

Compliance with ethical standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.