In Silico Prediction of Food Properties: A Multiscale Perspective

Several open software packages have popularized modeling and simulation strategies at the food product scale. Food processing and key digestion steps can be described in 3D using the principles of continuum mechanics. However, compared to other branches of engineering, the necessary transport, mechanical, chemical, and thermodynamic properties have been insufficiently tabulated and documented. Natural variability, accented by food evolution during processing and deconstruction, requires considering composition and structure-dependent properties. This review presents practical approaches where the premises for modeling and simulation start at a so-called “microscopic” scale where constituents or phase properties are known. The concept of microscopic or ground scale is shown to be very flexible from atoms to cellular structures. Zooming in on spatial details tends to increase the overall cost of simulations and the integration over food regions or time scales. The independence of scales facilitates the reuse of calculations and makes multiscale modeling capable of meeting food manufacturing needs. On one hand, new image-modeling strategies without equations or meshes are emerging. On the other hand, complex notions such as compositional effects, multiphase organization, and non-equilibrium thermodynamics are naturally incorporated in models without linearization or simplifications. Multiscale method’s applicability to hierarchically predict food properties is discussed with comprehensive examples relevant to food science, engineering and packaging. Entropy-driven properties such as transport and sorption are emphasized to illustrate how microscopic details bring new degrees of freedom to explore food-specific concepts such as safety, bioavailability, shelf-life and food formulation. Routes for performing spatial and temporal homogenization with and without chemical details are developed. Creating a community sharing computational codes, force fields, and generic food structures is the next step and should be encouraged. This paper provides a framework for the transfer of results from other fields and the development of methods specific to the food domain.


INTRODUCTION
Food engineering and modeling tend to be more difficult than in other areas as foodstuffs are inherently complex objects; biological processes initially define their organizations, structures, and compositions. They evolve rapidly during harvesting, transportation, and distribution. Food processing profoundly changes all the properties of food: composition, structure, physical states of the constituents, degradation, and production of new chemical compounds (Datta, 2016). In other words, successful modeling should capture the evolution of the physicochemical properties with the degree of processing but also the effects of the initial biological variability and of supply chain fluctuations. Storage, preservation, final preparation and digestion of processed foods also profoundly alter these properties (Bornhorst et al., 2016). Modeling in food addresses usually only one part of the complexity: only the evolution of an "averaged" food product is sought and not the variability between food products. Thermodynamic, transport and mechanical properties are therefore "averaged" over a representative structure, composition, and assumed constant or dependent on macroscale state variables (e.g., temperature, water content, etc.). With the rapid progress in modeling and simulation in other engineering disciplines, including materials, polymer science, and soft matter, this chapter argues that the implicit paradigm of averaging any aspect of food transformation may be challenged in the future. The lack of extensive and freely accessible "mechanistic" food properties databases motivated the quest for alternatives to classical modeling (Nesvadba et al., 2004). When they exist, they do not consider compositional and structure effects essential to processed foods characteristics. Moreover, simulating flows, transfer, reactions, creation of structures and their mechanical properties in realistic conditions remains beyond the reach of conventional approaches without adequate parameterization on real-life food products (Roos et al., 2016).
Directly calculating physicochemical properties of foods from pure phase data or atomistic structures paves the way to the direct simulation (from first principles) of their processing and preservation regardless of the food or its process. In this new perspective, the incomplete macroscale model in food is supplemented by an explicitly given microscopic model. Such a framework differs from earlier semi-empirical models tabulated in Ref. (Gulati and Datta, 2013). This article reviews a selection of techniques enabling the emergence of macroscopic behaviors from the knowledge of structures or molecules and their potential applications in food engineering. Indeed, the popularization of high-throughput measurement techniques encourages the development of modeling techniques capable of using directly as inputs 3D images from X-ray microcomputed tomography, laser-scanning confocal microscopy in UV, visible and midinfrared ranges Raman confocal microscopy). Ideally, the approaches need to be sufficiently robust and tailored to accommodate various strategies (fully or partially in silico) to replace food missing properties and incorporate chemical and structural food characteristics. The growing development of equation-free multiscale computation (Kevrekidis and Samaey, 2009), along with the widespread use of parallel and cloud computing, makes the passthrough to the scale where the decision is taken possibly without requiring an explicit macroscopic-level step. "Food system analysis" at any step of its transformation, storage, oral processing and digestion may nourish innovation in food but also reaches external communities, including medical communities, industry and authorities. Linked decisions across the scales with arbitrary levels of entanglement (Gear et al., 2003) could be used to evaluate the impact of existing or future practices.
The manuscript borrows several concepts and techniques that we have developed to fill in the gaps in continuum numerical approaches, and more generally, to introduce chemical information, which does not fit well in evolutive equations, governed essentially by conservation equations. From a thermodynamic point of view, we address in this paper problems which are entropy-driven and therefore not directly related to heat transfer or mechanical work. The vantage point adopted in the proposed computational framework suggests that two multiscale problems coexist: spatial and temporal. Abundant literature describes the first one. It is approached from the perspective of mathematical homogenization and applied to the prediction of diffusion coefficients. This technique should not be confused with techniques for averaging over volume elements. It does not focus on the observation scale to define heterogeneities but instead defines a very large equivalent medium in which local variations become negligible. The mathematical tool is asymptotic analysis. The same concept is required to predict temporal evolutions at the food scale from scales close to the vibration of atoms or molecules. The coherent macroscopic behavior emerges from the interactions of microscopic agents (solutes, macromolecules, cells). Therefore, spatial and temporal microscopic details are not smoothed out but rather integrated up to the scale of the food, usually using a procedural computing approach. The feasibility of combining chemical-level, microscopic-level and food-level is illustrated in the interactions between the food and its packaging. The chosen methodologies have sufficiently low computational cost to operate in loops and support food production problemsolving: e.g., resource and waste minimization, valorization byproducts, safety, security and health issues (Bazilian et al., 2011). But what are the challenges and opportunities of modeling across scales?

The Scale Problem
It is tempting to think that all scales can be modeled continuously from atoms to the cycle of production, distribution, and consumption of food by a population of consumers. The grounds why this is not possible in a classical sense are illustrated in Figure 1 regarding the farthest scales from the observer (human scale): either in the direction of atoms or on the scale of a sufficiently large system (including more than one single food).
At the finest scale, there are more molecules/entities in a 1 L water bottle (55.6 moles, i.e., 3.3 × 10 25 molecules; i.e., 10 26 atoms) than stars/entities in the observable universe [from 10 22 to 10 24 according to the European Space Agency (ESA, n.d.)]. The most extensive simulation of a biophysical system at the atomic scale on an exascale supercomputer (130,000 cores) included only 10 9 atoms/entities (Jung et al., 2019). If it were water, a cube of 215 nm side would contain this number of atoms. A water cube of 42 nm edge would require a more reasonable power of 1,000 cores. These volumes are insufficient to represent any food transformation (freezing/thawing, cooking, flow in a pipe), but they are sufficient to access fundamental properties at the nanometric scale. For example, the coalescence of two drops of 58.5 nm in diameter has been recently simulated at the atomic scale (Perumanath et al., 2019). The role of thermal fluctuations in the construction of a common liquid-gas interface has been demonstrated. This example illustrates three critical properties: • Structural chemical information can only be maintained on scales commensurable with molecules and segments of polymers; beyond, structural and chemical information will have to be encoded using the principles of thermodynamics and mechanics on a larger scale.
• Statistical physics or mechanics can describe phenomena below the thermodynamic limit (TL), that is, when thermal fluctuations are significant. TL is much larger than molecules and also applies to colloidal systems (Frenkel, 2014). • The union of systems smaller than TL does not allow the reconstruction of a macroscopic system by sampling. The property of additivity or extensivity is recovered only beyond TL.
There are no absolute rules to define the critical length scale for TL; TL tends to be longer in the presence of long-range forces or for systems with large surface-to-volume ratios in monophasic systems (Sheehan and Gross, 2006;Frenkel, 2014). In emulsions, foams, and suspensions, TL is much larger and reaches few micrometers.
The continuum scale (CS) lies between TL and an upper limit called "large scale" (LS). The "single-scale" approach is associated with finite-element and finite volume techniques. At these scales, there is no need to assemble forces on discrete systems and to assess the concept of thermodynamic and mechanic equilibrium locally. Statistical formulations are cast in conservation equations and resolved generically with the machinery of differential and integral calculus. Therefore, the food properties and driving forces are effective and capture most of the microscopic details in one single formulation. The properties are allowed to vary either with position, time, and intensive variables such as temperature, pressure, or local concentrations. The effective porous medium formulation describes well the effective heat transfer properties in heterogeneous food products such as dry foods. However, it is less efficient in deriving transport and mechanical properties as it does not consider phases spatial arrangement and morphology (Datta, 2007). Similarly, the corresponding driving forces (e.g., osmotic or swelling pressures) are well described for systems with low TL values (i.e., random mixtures of molecules). However, they may be more difficult to derive for systems presenting an organization at nanoscale or memory effects (i.e., dislocation, plasticity, glassy state) (Lucarini et al., 2020). Two examples easily illustrate the loss of details of the continuum approach at the food scale. The fracture of dry foods (e.g., chips) or rehydration of dry foods (e.g., cereals) will be described at continuum scale from the classical macro-strain without considering the anisotropic "micro strain" in cell walls and the effect of drying on the rigidity of cell walls. Similar microscopic effects occur when a native starch suspension is heated in a heat exchanger (Plana-Fattori et al., 2016). Gelatinization and swelling of starch granules will modify the viscosity locally and finally at a much larger scale. Microscopic gradients persist when two components with different viscosities are thermomixed such as during chocolate conching (Greiner et al., 2014).
Complex food simulators can be assembled from low-level mechanistic components for one particular food product. The concept of "large-scale" (LS) occurs when the same model is used for a broader range of food products. Its strict definition is more elusive and would correspond to situations where the underlying assumptions of at least one constitutive model cease to be verified. This model can be involved in either the prediction of effective properties or driving forces. The reasons could be: the raw biological materials are different (different variety, maturity, different tissues), the morphologies are different (different distributions for dispersed systems), the properties of the continuous phase are different (e.g., the addition of a thickening agent, fat replacers), different pretreatments (higher shrinkage, compaction), aging effects, etc.

Continuum vs. Discrete Descriptions
The core of this article is to highlight the different strategies to derive properties for food considered not as a uniform rigid body but as a soft matter organized at different scales. The choice of considering a discrete collection of interacting entities (molecules, phases) or an equivalent homogeneous continuum depends not only on the problem and the choice of the end-user but also on whether the chemical information needs to be considered explicitly or not. At this point, it is essential to note that no measurement is required for none of these methods beyond some information on the food composition and the spatial organization of phases herein. Counterintuitively, asymptotic analysis linking microscopic descriptions and standard mathematical formulation in simulation software (Comsol Multiphysics, OpenFoam, Fluent) are more abstract and require a minimum mathematical background. Conversely, particle methods used below TL (e.g., molecular or coarse-grained approaches) or above TL (Lagrangian descriptions) are more flexible and look more physical in particular to describe surface effects (coalescence, wetting, capillarity, adsorption, adhesion). The article presents both approaches, and both are sources of unavoidable errors. Replacing measured properties by computed properties (see calculable properties in Table 1) requires rigorous constructions and proper validation. The present article illustrates both methodologies when the sought property (transport coefficient, thermodynamic property) is: • either a scalar or a symmetric tensor (anisotropic diffusion without no drift); • either defined as a unique value or as a distribution of values.

Problems in Food that Could be Resolved by Modeling at Multiple Scales
The concepts of modeling across the scales are still emerging in food, and the examples covering from molecular to very large scales are scarce in the literature. The authors developed experience in the field initially from molecules (Vitrac et al., 2006;Durand et al., 2010;Gillet et al., 2010;Vitrac and Gillet, 2010) to consumer exposure  with the intent of encouraging the industry and authorities to use high throughput predictive techniques for compliance testing and risk assessment. From our experience, the main obstacles to the generalization of large-scale modeling in food applications are: -the lack of case studies showing how to combine deterministic and non-deterministic modeling in food engineering curriculum; -the lack of a computational framework to combine tools obtained with different approaches (statistical physics, deterministic, stochastic); -the difficulty to separate the phenomena occurring at different scales; -the resistance of some practitioners and professionals who think that food cannot be reduced to algorithms.
In related areas of food engineering, we argue that "modeling across the scales" has already produced very encouraging and competitive results. challenges identified in the literature, which can be considered resolved or in the process of being resolved. The list is not exhaustive and indicates the dynamics of this theme for the coming decade. It is important to note that the most crucial point is not so much the multiscale modeling itself but the ability to reverse engineer the property being studied. The problem is no longer: "which property I need to supply to my model, but what does the model tell me about the property that could achieve the desired effect." The introduction of the chemical and biochemical information allows tackling formulation and shelf-life determination problems from a mathematical point of view without the constraint of available raw materials or technologies. Food design becomes then a virtual activity supported by reusable libraries of self-similar problems. The chosen solution could be refined without extra costs for specific markets and consumers. Predictive property modeling methods are well documented in the literature but have not been implemented in the multiphysics software used by the food science and food engineering community.

THE PROBLEM OF SPATIAL HOMOGENIZATION
Spatial homogenization is essential for modeling from 2D and 3D images. This section describes the general principles of spatial homogenization, as shown in Figure 2. The approach is applied to transform heterogeneous structures into effective media, whose properties can be encoded into commercial simulation software. Homogenization theory was elaborated in the 70s and 80s for composite materials; its application to structured foods (solid and semi-solid) is a natural extension. Classical homogenization is preferred to alternative approaches as it has received more formal proofs than other methods such as volume averaging and does not suffer complex entropy constraints for closure, such as the hybrid mixture theory (Battiato et al., 2019).

General Formulation
Homogenization (time-or spatial-) aims to calculate effective properties for processes spanning over several time scales (e.g., combing rare events along with fast evolutions) or over heterogeneous structures. Any spatial-homogenization problem (not time-dependent) along the dimensionless coordinate x (socalled slow variable) is recast as: where L() is a differential operator (linear or not) and y is the fast variable (fine-scale mapping) inducing the local fluctuations of the property. u ϵ is the field (e.g., concentration, strain, temperature) calculated at the scale of the microstructure, so-

Rheology
• Newtonian Keffer et al. (2004) Due to the limitations of the accessible time frames, transport coefficients (viscosities, diffusivities) are calculated at very high shear rates far beyond those applicable experimentally • Non-Newtonian Cummings and Evans (1992) Mechanics • Elasticity Datta (2016), Cummings and Evans (1992), Misof et al. (1997) T g values are overestimated by molecular modeling due to the high cooling rates applied • Glass transition temperature: T g Han et al. (1994) Elasticity can be directly studied at the molecular scale. Plasticity and ruptures can be studied out-of-equilibrium on small systems under strong strain rates. The observed irreversibility can differ from those observed macroscopically due to finite size effects and studied time scales shorter than relaxation times • Thermo-mechanical behaviors Shenogina et al. (2012) • Food structuring Dickinson (2012) Transport • Permeability Venable et al. (2019) Random walks can be studied for diffusivities greater than 10 −14 m 2 ·s −1 . See Section 3.2, Section 4.3 and Section 4.4 • Diffusivities Vergadou and Theodorou (2019), Dubbeldam and Snurr (2007) • Solubilities Vergadou and Theodorou (2019) Thermodynamic properties • Phase diagram Fan et al. (1992), gelation Lu et al.
Invoking Eq. 1 on the entire macroscale domain can be particularly expensive (typically O −n or more) and possibly not achievable for very small scales ϵ ≪ 1 or when the properties need to be updated regularly. The cost is usually reduced by assuming a separation of scales between the microscopic details and the macroscale domain simulated (ϵ ≪ 1). It implies that beyond a critical scale ϵ c , u ϵ converges to the macroscale solution U, which can be calculated at a much lower cost by solving the homogenized version of Eq. 1: Finding an explicit homogenized equation L(x) from microscopic details may be a difficult mathematical problem and requires an adequate closure problem. Such equation exists as soon as y ≫ ϵ c ϵ (see Figure 2B) and is associated with effective properties, whose meaning is similar in various physical problems, but whose values reflect different microstructures and underlying phenomena, such as transport on distances smaller than the mean-free-path of molecules or transports controlled by rare events (e.g., trapping, surface adsorption, interactions with macromolecules). Many software packages can evaluate U(x) using finite differences, finite volumes, finite elements, or higher-order spectral techniques based on proper property values. It is worth noticing that the volume source/sink term f(x) appearing on the right-hand side of Eqs 1, 2 is assumed to be macroscopically imposed and, therefore, not subjected to microscopic fluctuations. Otherwise, the homogenization procedure cannot converge to a unique solution.

Principles
The microscopic quantities a ϵ (x, y x ϵ ) (static) or a ϵ (x, y x ϵ , t) (dynamic) may be diffusivities (this example), permeabilities, elasticity moduli, reaction rate constants, etc. They converge to the effective properties when y ≫ 1. It is worth noticing that homogenization does not aim at averaging the microscopic quantity a ϵ (x), but the underlying conservation laws instead. By considering transport phenomena and due to the intertwining with other microscopic phenomena such as ad/ absorption, a ϵ (x) should only be envisioned as the property Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 relating the microscale flux j ϵ (x, t) a ϵ (x)∇u ϵ (x, t) with the microscale driving force ∇u ϵ (x, t).
The corresponding microscopic advection-diffusion problem reads: Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 8 By assuming that a and b vary smoothly over a periodic cell (of volume V) larger than ϵ c , Eq. 3 homogenizes as: where b is the average of b over V and A is the homogenized operator (stiffness matrix), defined as z zx i [a ij z zx j ] ∇ · [A∇], and capturing the effective transport properties along each coordinate.

Intuitive Homogenization Procedure
The simplest case arises from the homogenization of Eq. 3 in onedimension (n 1) at steady-state and in the absence of advection, which reads: By setting F(x) x 0 fdx and C ϵ a constant, the double integration of Eq. 5 leads to: By invoking the separation of scales for ϵ → 0 and that the problem is periodic along y, one gets C Thus the homogenized transport coefficient A 〈 1 a ϵ 〉 −1 , where 〈X〉 represents the volume average of X. The following result can be generalized for arbitrary boundary conditions and dimensions. The homogenized matrix corresponding to A ϵ a ϵ  (Zhikov et al., 1994)]: The presented results generalize standard approximations for effective media, such as the Maxwell Garnett and the Maxwell-Rayleigh approximations for spherical and non-spherical inclusions. For symmetric positive definite matrix, the following inequalities (known as Voight-Reiss' inequalities) hold: The first equality Though the studied transport equation does not include any thermodynamic consideration, particular cases of two-phase systems can be treated with the previous procedure. Soft and hardcore inclusions are coded as A ϵ (x) 0 and A ϵ −1 (x) 0, respectively, with nonzero values outside inclusions:

Asymptotic Expansions in the Absence of Microscopic Discontinuities
Explicit solutions exist only in rare cases such as periodic stratified media, but they cease to operate in large dimensions. Alternatively, the microscale solution u ϵ (x, y) can be envisioned as a first-order approximation of the solution of the homogenized problem U(x): Equation 10 accepts a statistical interpretation: . The exact calculation of expansions involves total derivatives. The procedure is accelerated by highlighting the composition properties of the microscale differential operator: where: The expressions Eq. 11 and Eq. 12 applied to approximation Eq. 10 (by assuming that the macroscale solution is twice differentiable up to the boundary) lead to: The separation of scales (ϵ → 0) enforces that Δ 1 δ u −Δ 2 U, which can be resolved as a periodic boundary value problem to the independent variable y. By denoting its general solution {ϕ k } k 1...m , one gets δu(x, y) ϕ k (y) zU zx k . This solution introduced in Eq. 13 and averaged over all possible solutions leads to the following homogenized coefficients: a ij 〈a ij + a ik zϕ j zy k 〉. A second-order approximation can be derived for periodic macroscale solutions U. The results are very general, but it requires that U is of class C 2 . This property fails to be verified if the homogenized medium includes barriers (e.g., cell walls), strong partitioning, or large jumps in transport properties.

Homogenizing With Microscopic Discontinuities via a Variational Formulation
In the presence of jumps in a ϵ (x), the homogenizing problem is better solved using a weak formulation of the partial differential equations. The sufficient condition is that u ϵ fulfills Dirichlet or Neumann boundary conditions at jumps. For a sufficiently smooth and periodic a ϵ ij (x, x ϵ ), the solution of the microscale problem at steady-state minimizes the standard variational minimization problem: The solution of Eq. 14 is inferred efficiently via a top-down approach sketched in Figure 2C. Standard finite-element techniques (triangular mesh) are used to solve the macroscale problem. The effective macroscale operator herein is approximated by solving the microscale problem only on a small subset of elements taken at quadrature points. This strategy is known as the heterogeneous multi-scale method introduced by W. E. and B. Engquist (Engquist and Huang, 2003) and discussed in Ref. (Weinan et al., 2005). The overall cost does not increase as ϵ decreases and generalizes naturally to time-dependent and nonlinear problems.

Limitations
Capturing the microscale statistics adequately is essential as the whole approach relies on an approximation of the separation of scales (ϵ ≪ 1), for which A evolves only smoothly with the position. Due to the independence of microscopic and macroscopic scales, there is no need to apply microscale fluxes with intensities as large as those applied at the macroscale. Indeed, macroscale fluxes describe an evolution of the macroscopic system, whereas microscale ones describe its evolution at some equilibrium states. The concept of equilibrium differs here from strict thermodynamic equilibrium as it assumes an equilibrium of microscale fluxes/ forces without enforcing stronger concepts, such as microscopic reversibility or time reversibility. The concept of equilibrium can be approached numerically by introducing a micro time scale and by refining the microscale solutions iteratively until reaching a steady state.

Formulation With Non-Separable Scales
The main interest of homogenized coefficients is that they can be used for a broad range of problems, different boundary conditions, and source terms. When microscopic and macroscopic scales are not fully separable, the concept of homogenized coefficients is, however, no longer applicable. Nevertheless, the concept of the homogenized conservation equation still holds. The term A∇U in Eq. 4 can be replaced accordingly by the macroscale flux J averaged at each iteration over the microscopic structure.
In the general case and the presence of fluctuations of microscale fluxes, J involves both an averaging over the cell and a temporal averaging.
Significant works exemplifying the techniques of homogenization in food or in structures of general interest are listed in Table 3. We regret the lack of implementation of homogenization techniques in standard software packages, as

Effective property References
Diffusion coefficients in porous media Auriault and Lewandowska (1997) Multicomponent gas diffusion in apples Ho et al. (2011) Permeabilities through spheres and polyhedron packaging Boutin and Geindreau (2010) Thermal conductivity for particulate foods Chinesta et al. (2002) Elasticity modulus from confocal microscopy images of food Kanit et al. (2006) TABLE 4 | Examples of direct description of the behavior of multiphasic food properties.

Method Examples (references)
Hybrid mixture theory (HMT) • Rheological properties of blueberry puree Nindo et al. (2007) • Stresses during drying of corn kernels Takhar et al. (2011) • Phase-change and thermomechanics during frying Bansal et al. (2014), Sandhu and Takhar (2018), during extrusion Ditudompo and Takhar (2015), during freezing Zhao and Takhar (2017) • Unsaturated fluid transport in food Takhar (2014) Dissipative particle dynamics (DPD) DPD methodology is usually applied at the supramolecular scale, where the properties of the food are less specific. A general methodology has been developed to devise interaction parameters in DPD models from arbitrary dynamics (at atomistic scale or from another coarse-grained simulation). Eriksson et al. (2009) The applications are numerous, from evaluating the permeability of the skin to chemicals Otto et al. (2018) to digestion phenomena Palkar et al. (2020) well as the lack of standard food structures to develop or test algorithms.

Hybrid Mixture Theories and Discrete (Lagrangian) Methods to Describe Multiphasic Systems
Averaging and homogenization methods (Battiato et al., 2019) have the same objective to characterize the macroscopic properties of heterogeneous media by applying filtering procedures to instantaneous equations similar to statistical methods in turbulence. A great deal of work, including those of S. Whitaker (Whitaker, 1986;Whitaker, 1996;Whitaker and Whitaker, 1999), has shown how the treatment of the classical local equations of mechanics could be reduced to averaged equations, for example, obtaining the law of Darcy from the Navier-Stokes equation. Closing the averaged equations requires defining a priori characteristic scales of time and space. In addition to the properties that are difficult to average, the constitutive homogenized equations do not exist for all problems and all situations. Several features are challenging to implement, including local thermodynamic equilibrium with possibly chemical reactions, non-Newtonian flows, flows in elastoporous media, large deformations (swelling or shrinkage), flows and fluxes affected by macroscopic strain, inertial flow in cavities, multiphasic turbulent flows, etc. Two families of approximations and simulations have been proposed to address non-heuristically previous limitations, with examples listed in Table 4. The first is the hybrid mixture theory (HMT). This macroscale theory uses constitutive equations, which are neither homogenized nor averaged over the microscale but directly formulated at the macroscale and subjected to the second law of thermodynamics. The approach has been first developed by Hassanizadeh and Gray in seminal papers (Hassanizadeh and Gray, 1979a;Hassanizadeh and Gray, 1979b;Hassanizadeh and Gray, 1980). All conservation equations are averaged using the Whitaker procedure to obtain macroscale field equations, but each variable remains defined precisely in terms of its microscale counterpart. The structure of the theory of mixtures of Bowen (Bowen, 1980;Bowen, 1984) has been implemented to consider the impenetrable/immiscible properties of porous media. Constitutive relations are required to close the formulation, but the approach can be seen as a natural extension of the classical theories of gas mixtures. The formulation is well adapted to describe immiscible and incompressible flows in deformable and rigid solids. For example, capillary pressure can be derived from the free energy of the mixture as a function of temperature, the strain of the solid, and fluid volume fractions. However, HMT proposes closure only for the macroscopic scale and not across the scales. One practical consequence is the lack of conservation of interfacial areas in fractal structures. Hierarchical averaging procedures have been proposed to estimate wetting properties from the microscopic geometry of pores Miller and Gray, 2005;Gray and Miller, 2006;Miller and Gray, 2008;Gray and Miller, 2009). Though these theories are complex, they highlight the risk of consistency loss in too coarsened models. Keeping thermodynamic constraints across the scales is essential.
The second family relies on discrete methods to represent inertial flows and complex interactions between solids and fluids. They offer a competitive alternative to the complex description with continuous methods of free surfaces and situations where surface tension dominates (e.g., foaming, coalescence). Discrete mechanics is based on connected objects, points, segments, surfaces, volumes (beads). The positioning of these objects in space is unnecessary, the referential is only local, and the objects are located in relation to each other (Lagrangian description). The corresponding methodologies are essentially meshless, and no spatial discretization is required. Because the variations are numerous, only two methods are briefly presented: discrete element method (DEM) and smoothed particle hydrodynamics (SPH). Both methods rely on collision operators. DEM is the oldest methodology (Cundall and Strack, 1979) and describes dry cohesionless granular flows while accounting for the exact shape of particles. Conversely, SPH is a particle-based continuum method well suited for flows and deformation problems. Continuous media are represented by overlapping particles representing a volume of fluid or solid. The methodology has been initially developed for astrophysical problems (Gingold and Monaghan, 1977;Cleary and Monaghan, 1990). The methodology outperforms any continuous methods to describe free surface flows.
Inherently in the original formulation, the concept of a "blurred" particle suffers several challenges at boundaries and in computational efficiency. It is, therefore, common to represent also solid interfaces with fixed particles. As DEM and SPH simulation methodologies do not rely on statistical ensembles, their original flavors are unsuitable for entropy-driven phenomena. Fluctuations are recovered in Dissipative Particle Dynamics (DPD), which targets all soft-matter and flow phenomena occurring at nanoscale, typically 10 nm-10 μm and 1 ns-1 ms (Español and Warren, 2017). The DPD model, introduced initially by Hoogerbrugge and Koelman (Hoogerbrugge and Koelman, 1992) and subsequently equipped by Español and Warren with a proper statistical mechanics model (Español and Warren, 1995), gained rapid popularity due to its algorithmic simplicity and versatility. Quoting Español and Warren, "just by varying at will the conservative forces between the dissipative particles, one can readily model complex fluids like polymers, colloids, amphiphiles and surfactants, membranes, vesicles, and phase separating fluids. Due to its simple formulation in terms of symmetry principles (Galilean, translational, and rotational invariances), it is a very useful tool to explore generic or universal features (scaling laws, for example) of systems that do not depend on molecular specificity but only on these general principles". Hence, the interactions in DPD schemes conserve both linear and angular momentum locally. Finally and to be exhaustive, several SDPD "smooth DPD" strategies (Uma et al., 2011;Shang et al., 2012; have been explored to get an overlap with SPH methodology.

THE PROBLEM OF TIME HOMOGENIZATION
Quoting Gorban (Gorban et al., 2006), the first coarse-graining algorithm was proposed one century ago by Paul and Tanya Ehrenfest in their paper for the scientific Encyclopedia (Ehrenfest et al., 2014). The operation proposed by the authors transforms the continuous probability density in phase space into a piece-wise function averaged over periodic cells. At the time, averaging over periodic cells was revolutionary as it enabled the system's entropy to increase even after a long time and overcame the previous limitations offered by the Liouville equation for closed systems. Nowadays, coarse graining over grids is everywhere in statistical physics both at equilibrium and non-equilibrium. Complex transport phenomena (flows, diffusion) can be studied with various dissipation coefficients (viscosity, diffusivities) and decoupled with proper filtering and subgrid modeling.
Strategies that span the time scales in condensed phases follow a different tack than those designed to span spatial scales. Generally, there is a wide gap of about ten orders or more of magnitude in time between processes like atomic vibrations, which ultimately limit simulations at the atomistic scale, and complex diffusion processes that control many dynamic properties such as shelf-life, safety, release digestion, oil uptake, etc. Any robust multiscale modeling strategy should be equipped with both spatial and time homogenization. The key ingredients for time acceleration and ergodicity are: • to get a proper probabilistic description of first passage times, • to define the parent stochastic transport mechanism (diffusion, vibrations), • to set the adequate threshold (energy, filling time, working time), • to make the time scale emerge from the considered statistical ensemble or the method of sampling, • to get a proper methodology to sample times continuously.
This section introduces the concepts of transition state theory as a common framework and illustrates several strategies to identify first passage or first hitting times.

Transition State Theory
Previous methods can zoom spatially continuously from the macroscopic scale down to atoms. Time scale is, however, not treated similarly by modelers as its definition may be twofold. It can date the events since the beginning of the simulation time or our clock time. Still, it can also represent only a reciprocal frequency without specifying the origin of simulation time. This last definition prevails in simulations involving thermodynamic/statistical ensembles, where time emerges directly from free energy variation between two macroscopic states. The dichotomy of time interpretation above a gap ranging from microseconds to milliseconds is shown in Figure 1A. It can be interpreted easily as follows. Macroscopic times define the order of events, while microscopic times are random, and only the durations between events are statistically defined. This memoryless process is called Poissonian and can be sampled causally from a single series of events (e.g., the trajectory from a single particle) or a collection of events (e.g., collection of particles). The difficulty arises from sampling the frequency of rare events (e.g., trapping-release), whose effects control the overall rate or kinetics.
In cohesive systems such as entangled polymers (e.g., packaging walls), structured food systems (e.g., emulsions) or in the presence of chemical reactions (chemisorption, cage effects, etc.), several time scales coexist due to the broad correlations between consecutive phenomena. After coarsegraining, such phenomena can be understood as a macroscopic trap-escape sequence (i.e., collection of macroscopic states), whose frequencies are very low compared to events controlled by the local thermal vibration of atoms or molecules, whose corresponding states are called micro-states. The separation between macro-and microstates on the surface of free energy of the studied system (e.g., mixture) is equivalent to the separation of two different time scales.
Transition-state theory (TST) approximates the rates of a transition between two neighboring macrosites/macrostates, denoted 1 and 2, separated by an asymmetric energy barrier ΔE 1→ ‡ and ΔE 2→ ‡ as: where Q ‡ and {Q} i 1,2 are the partition functions of the transition state and the considered macrosite i, respectively. k B and h are the Boltzmann and the Planck constants, respectively; T is the absolute temperature. As partition functions (canonical) are interpreted as the sum of the probabilities in the form of exp(− ϵ j k B T ) with ϵ j the energy associated to the state j, they are likely to be also a function of T. TST has been initially devised to calculate reaction rates along a reaction coordinate (Eyring, 1935), and it has been subsequently generalized to the diffusion in glassy polymers (Gray-Weale et al., 1997).
The concept of transition state between macrosites 1 and 2 is essential, and it assumes that the transition paths 1 → 2 and 2 → 1 follow the same route and cross a saddle point with negative curvature, denoted ‡, where the origin of the motion is indiscernible. The concept of saddle point in classical mechanics assumes that there is a separable coordinate along with the transition occurs (the curvature of the saddle point is positive everywhere except along one coordinate), or more precisely where the flux of probability is maximum in both directions. It corresponds to a favorable collective displacement of atoms or a favorable geometry configuration, suggesting a "tunnel" or a "Red Sea mechanism" connecting two states. By denoting {p i } i 1,2 the equilibrium probability of state i, the principle of microreversibility reads at any time: Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 Once the connectivity of the sites and the frequencies of the corresponding transitions established, the equilibrium distribution can be determined by integrating the macroscopic mass balance equations until to reach steady state ( dp eq i dt 0): The asymptotic residence time at any microsite is given by the harmonic average of transition rates: At this stage, it is essential to notice that the TST and the methodology treat time as a continuous variable. Space is replaced by a Lagrangian/Hamiltonian representation via a surface of potential or free (preferred) energy. A continuous representation in time and space can be recovered by averaging the residence times over all possible sites: A form of Eq. 21 is applied free-volume theories used to model trace, solute, and mutual diffusion in condensed phases. The probability of a translation (transition) is assumed to be proportional to the probability of a critical amount of free volume accumulating in the immediate vicinity of the considered molecule. The exact position of the atoms and the concept of adjacency do not need to be explicated. Similarly, the dual-mode sorption model assumes that the translation may occur in two different micro-environments (microsites). It might be tempting to consider that Eq. 16 describes a temperature activation of the Arrhenian type. This is not the case because partition functions also contain an enthalpy component. The prefactor describes the effect of thermal agitation as in Stokes-Einstein's law by a simple temperature effect. The other parameters carry the energy barriers and the effects of thermal expansion. Finally, it deserves to note that multiple transition state theory assumed that the transition rates or, more precisely, that the free energy of interactions is distributed.

Examples
The following examples illustrate mechanisms involving multiple time scales and the methodology to simulate them by dropping unnecessary details.

Justification of the Separability of Microscopic and Macroscopic Time Scales to Predict Transport Coefficients
Diffusivities of gases and organic molecules in condensed phases (liquids or solids) can be derived from the secondorder stationary displacements of the center-of-mass (CM) of one single diffusant (or averaged over the trajectory of several diffusants) by molecular dynamics simulation. The displacements controlling the translation of CM due to thermal noise (collision with atoms of surrounding molecules) between times τ and τ + t, r CM (t, τ) are defined from the positions vectors (3 × 1 vector) of CM, x CM , relatively to the center of mass of the entire system (diffusant + host), x O : When the displacements of CM are truly independent, the quadratic displacement increases linearly with the lag operator t with 〈r CM (τ, t)〉 all t 0: where X T is the transpose of X . The diffusion matrix D(t) is positive-defined, whose eigenvectors code for the main direction of translation on time scale t. It accepts a unique Cholesky decomposition: where Λ(t) is a diagonal matrix containing the eigenvalues and where the columns of B(t) are the associated eigenvectors. Displacements fluctuate on the time scale t at a rate drCM(τ,t) dt equal to B(t)Λ(t) 1 2 X, where X is a random vector, each of whose components follows a reduced centered normal distribution.
Solutes larger than voids exhibit a preferential direction as shown for methoxybenzene (anisole) in Figure 3, so that the random trajectory looks like an extruded tube on large time scales ( Figure 3A). The solute seems to jump from one sorption site to the next, separated by distances ranging from 3 to 6 nm ( Figure 3B). During the different jump events, the solute progressively loses memory of its previous orientation ( Figure 3D). The macroscopic diffusion coefficient eventually emerges as the trace limit of D(t) when t → ∞ (note that t 1 f is the lag operator and not the clock time). By choosing the columns of B(t) as reference frame, one gets the following expression of the rate of increase of the mean-square displacement of CM: which leads to the practical formula of the macroscopic diffusion D(f 0): The exponent d ln(msd CM (t)) plays a central role in identifying the rate of convergence to the hydrodynamic regime, where all displacements are truly uncorrelated (exponent close to unity). Exponents lower than 1 reveal strongly correlated random walks and apparent diffusivities decreasing with simulation time. Such decreases are shown in Figures 3E-I Since the diffusion of a solute larger than voids in a rubber polymer is strongly related to the relaxation of polymer segments, which could think that the convergence of D(t) should obey some linear damping (e.g., Hooke's law) and consequently that diffusion coefficients could be easily extrapolated, when the associated time constant t 0 is known. This reasoning would lead S CM (f) to decrease in 1 f 2 and an exponential convergence to the hydrodynamic regime: (28) Figure 3I shows conversely that S CM (f) evolves from a highfrequency white noise (tumbling motions of the diffusant) towards a low-frequency white noise (macroscopic diffusion) by crossing an intermediate regime linear in 1 f (pink noise), separated by two frequencies denoted f 1 and f 2 . This noise corresponds to a random damping process with t 0 distributed between t 2 1/f 2 and t 1 1/f 1 , such reasoning leads to: These results have been used to build a solid theory of diffusion in polymers coupling the free-volume theory and the mode-coupled theory (see Section 4.4). They justify the universal cage-escape mechanism to describe random walks in condensed phases.

Diffusivities on a Surface of Potential: The Relationship Between Diffusion and Sorption at the Molecular Scale
The frequency f 1 ≤ 1/t 0 ≤ f 2 and, more particularly, its lower bound control the value of the diffusion coefficient. Ref. (Vitrac and Hayert, 2007) studied the effect of the distribution of sorption energies G i on D systematically. The frequencies f 1 and f 2 are connected with the jumping rates via TST, , with G i the local sorption energy at site i and G ‡ i → j the free energy barrier associated to the translation from i to j. The main results are illustrated in Figure 4.
The strong relationship between sorption free energy and diffusion can be easily understood on a periodic surface of free energy, comprising only two states separated by a fixed distance l ( Figure 4A). Although it may seem surprising to simulate a random walk on a two-state line, the distance visited is a random variable because the duration of jumps/transitions is also a random variable. The diffusion coefficient depends on two characteristics, the height of the smallest energy barrier G Δ and the partition coefficient K contrast between the two sites. Increasing K contrast while keeping the same average sorption leads to a dramatic decrease in D values. This effect was verified on random surface energies in higher dimensions. Beyond a critical dispersion value, preferred sites behave as attractors creating negative correlations in the long term. Simulated assumptions are close to what happens when the size of the diffusant is increased (Neyertz and Brown, 2004). On non-random surfaces (i.e., with channels), the reduction of D values is amplified by additional tortuosity effects. Because the results can significantly deviate depending on whether the microreversibility hypothesis is verified, it is preferable to restrict the TST approach to the molecular and supramolecular scales. TST-type modeling lends itself very well to the integration of trajectories by Kinetic Monte Carlo (KMC) methods. Classical molecular dynamics simulation and KMC methods can be combined to reproduce the effect of rare events over the long term (Neyertz and Brown, 2010). The methodology allows to obtain statistically independent but without guaranteeing the final D value if limiting events have been insufficiently sampled.

Diffusion in Emulsified Foods
The previous examples justify the concept of timehomogenization to derive macroscopic transport coefficients. Rare events control the entire dynamics and its activation by temperature. One can intuitively think that such phenomena do not exceed a few hundred nanoseconds and are irrelevant at the phase level of the food. This is erroneous because diffusion is a stochastic process in time and in space at all time and length scales. We used to model the diffusion in the hydrodynamic regime in discrete times as a self-similar skewed trajectory, where the direction of the translation and its span are random normal variates. Close to interfaces and, in particular, curved interfaces, the process ceases to be isotropic. By analogy with optics, the constant K contrast and the brusque change of diffusion coefficients at interfaces act as refractive indexes and modify the distribution of molecules on both sides of the interfaces. In dense emulsions with globules close to the percolation threshold, diffusants scatter and reflect towards the medium with the highest residence time. The translation from the medium with the highest diffusion coefficient controls the effective diffusion coefficient. These effects are illustrated in Figure 5 in 1, 2, and 3 dimensions with the following conventions based on the results of Ref. (Vitrac and Hayert, 2020). K is the partition coefficient between the continuous and the dispersed phase. r D D c D d is the normalized diffusion coefficient. For large K or low r D values, the molecules initially contained in the globule create a halo (i.e., close to a diffraction pattern created by an infinitesimal object) around the globule, with a little back diffusion to the reverse direction. At this stage, it is important to notice that K defines only the ratio of diffusant densities across the interface zΩ (x 0) and not the probability to cross the interface. For the configuration depicted in Figure 5A with a diffusant located initially at x x 0 , the ratio of probability to remain the globule after a dimensionless time Fo D d t x 2 0 is given by: The time limit of Eq. 30 is counter-intuitive and reflects the ability of the diffuser to make countless round trips across the interface but at different speeds. The local thermodynamic equilibrium condition imposes the density ratio (formally K) and continuity of the fluxes. Additionally, it is noteworthy that globule interfaces do not represent in this representation a transition state, as the density at the interface is not continuous. However, the global escape rate from each globule could be used as input parameters to simulate diffusion in dispersed media. Compared to Eq. 8, this technique accepts larger fluctuations in surface area, globule size and density, large drops in diffusion coefficients, and finally, partition coefficients very different from unity.

Oil Percolation in French-Fries
Cage-escape mechanisms can be found in various applications, including the percolation of oil in French-fry or potato chips-type products once the product is removed from the bath. Trapped air bubbles in the potato tissue reduce the accessible cross-section to oil and consequently slow down dramatically the capillary flow through crust defects. These effects are illustrated in Figure 6 and are based on simulations and experimental observations reported in  and in (Achir et al., 2010;Vauvre et al., 2014;Patsioura et al., 2015), respectively. We could think that the effective diffusion coefficient of oil should scale as: where here c is the surface tension, θ the contact angle, η o the viscosity of the oil, 〈r〉 the average radius of cells assuming roughly a cylindrical shape. The micro-model showed that the successive filling times of two cavities of radius r and length Δz and connected by a defect of radius r 0 and length l 0 increase considerably with each pass. If the oil-air interface is displaced without creating a bubble ( Figure 6A), the recurrence relationship between the filling times Δt i,fill is given by: with 0 ≤ ς < 1 a factor that controls the gas phase displacement mechanism [see Ref. ]. Due to the hexagonal shape of the cells and the angular details of the defects, the oil flows faster in the corners ( Figure 6B) and fills the defects/restrictions without displacing the gas phase. The destabilization of the liquid films in the defects is responsible for the sudden phenomenon of "snap-off" (collapse of two initially disjointed films), which will lead to the creation and trapping of bubbles in the defects.
In the presence of an air bubble trapped within the cavity, filling times diverge and Eq. 32 becomes: with ξ the restriction coefficient of the passage related to the presence of the bubble and υ the fraction of the volume of the cell (i.e., formally the cumulative work of the volume forces acting on the bubble) preceding that must flow before the trapped bubble is released. The comparison of Eqs 32, 33 enables to approximate the lag time as Δt i+1,lag ≈ 1 ξ 4 Δt i+1,fill for low values of ξ et υ. Experimental filling kinetics of connected cells in parenchyma cells confirm this evolution and the distribution of waiting times as shown in Ref. .
This example and the previous one illustrate the importance of time homogenization and the sampling of rare critical events. Indeed, snap-off and partitioning are elusive and equilibrium properties, respectively. Because its occurrence depends on history, its effects cannot be correctly described by a homogenization procedure, assuming a particular steady-state.

Molecular Modeling as an Enabler of Multiscale Modeling
Most of the food industrial and technological processes require including chemical details in models to reach predictive descriptions beyond qualitative ones via dimensionless models. It could be tempting to consider that at the bottom of the hierarchy of the modeling methods across the scales, we have the Hamilton's equations, their approximate potential energy functions, and the methodology to resolve them numerically via molecular dynamics (MD). Millions to billions of atoms can be simulated nowadays up to hundreds of nanoseconds but not and beyond and at the cost of a time step imposed by the vibration frequency of the hydrogen atom. Brute-force MD is part of the solution but not the solution. Other sampling methodologies, such as thermodynamic integration, particle insertion, non-trivial statistic ensembles, coarse-graining, must be combined to reduce the overall cost of any information added at the atomistic or molecular scale. The blueprint for turning the dream into reality was hypothesized in the eighties by Herschel Rabitz (Rabitz, 1989). The outline of the modeling hierarchy to reproduce bulk observations from atomistic simulations is shown in Figure 7.
As a general principle, all simulations should validate that the density and pair radial distributions (or structure factor) of pure components are well verified before studying the system's dynamic response alone or in interactions. This section reviews the experience explored by the authors on molecular modeling and its applicability to food engineering cases. The field is evolving very fast and is the most promising for the future. As already argued by authors (Nguyen et al., 2017a), it starts to be less expensive and faster to launch thermodynamic calculations at the molecular scale than doing the experiments. As discussed in Chapter IV of Ref. (Tadmor and Miller, 2011), molecular dynamics of large systems can be combined or compared with many conventional techniques, such as linear elasticity, nonlinear thermoelasticity, hydrodynamics models of complex fluids, stochastic models with multiple scales, gas-kinetic scheme, quasi-continuum method.

Principles of Molecular Modeling
Molecular modeling (MM) is an umbrella term covering various techniques, where forcefields bring "life" to simulations, as illustrated in Figure 8. These techniques can be rooted in the seminal works of Martin Karplus, Michael Levitt, and Arieh Warshel, whom the 2013 Nobel Prize has awarded in Chemistry for "the development of multiscale models for complex chemical systems." A mechanical work leading to displacements of atoms or entities (blobs, group of atoms) can be derived from this potential independently of the path followed by the entire system. Only initial and final steps at the end of each time step count. The considered potentials (between two or more particles) may have several interpretations, quantum or classical, related to hardcore or softcore interactions. This section is not exhaustive but gives an overview of the most relevant techniques for food applications.

Forcefields
Molecular modeling uses forcefields to calculate the pair (distance), angular/planar, quadrupole, and octupole interactions between atoms, as sketched in Figure 8A. Most often, forcefields rely on classical Newtonian mechanics. Forcefields representing chemical bonds are either fitted to experimental data (type-I forcefields such as bead-spring model) or fitted to quantum calculations (type-II forcefield with additional corrections to mimic the vibration spectra of molecules). Non-covalent potentials are associated with van-der-Walls and Coulombic interactions defined by partial charges usually applied at the center of atoms. Some software packages, such as LAMMPS (Sandia National Laboratories) and Materials Studio (Dassault Systèmes BIOVIA), offer flexible formulations for soft and solid-state materials and coarse-grain systems. Others, such as Gromacs and Charmm, offers forcefields optimized/specialized for biological molecules, with an emphasis on proteins and lipids. Integrated environments (Avogadro, MedeA LAMMPS from Materials Design, Materials Studio) facilitate the many operations (setup, post-treatments) needed to build and interpret molecular simulations. Some generic visualization tools such as VMD and Ovito offers versatile visualizations for multiscale simulations. For

Methods to Explore the Potential Surface
This section introduces briefly the different principles of molecular mechanics. Additional details on forcefields, integration schemes and inversion procedures can be found in the reference textbook of Frenkel and Smit (Frenkel and Smit, 2002). These details are essential to understand the parameterization of open-source (Gromacs, LAMMPS, Charmm, CP2K) or proprietary/ commercial software packages (NAMD, Tinker, Materials Studio). At this stage, it is worth noticing that classical oscillators do not behave as quantum ones, whose more probable state coincides with a minimum of potential energy. The likeliest states associated with classical oscillators in dynamics correspond to kinetic energy minima (speed minima) and, therefore, maxima potential energy (see Figure 8C). According to the fluctuation- Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 21 dissipation theorem, the states corresponding to each maximum are unlikely, but their long-term average in molecular dynamics converges towards the likeliest state at equilibrium. This result will be obtained for microscopically reversible phenomena integrated in time with a time-reversible scheme. However, the convergence to realistic states may be very slow due to explicit numerical schemes, which are stable only for time steps much smaller than the period of the vibration of the lightest atom, typically 1 picosecond in the presence of hydrogens. When trajectories are not required but only configurations or distributions, Monte-Carlo sampling offer atoms better results. It can explore the surface energy in a non-physical manner and overcome any energy barrier separating two feasible states. When the sampling involves a sequence of trial moves (e.g., translation, rotation, volume change), each move must also be picked randomly to enable the reverse course and preserve the detailed balance. The macroscopic quantities are in both cases extracted from ensemble averages.

Coarse-Graining and Boltzmann Inversion
Due to the implicit time integration, acceleration can be gained only by reducing the number of freedoms coupled with longer time steps made possible by merging several atoms (or even molecules) into one heavier particle and therefore subjected to lower frequency vibrations. The iterative coarse-graining methodology is shown in Figure 9. The final goal is to preserve the convergence towards the final property by devising a proper mean force potential ϕ CG (R) from detailed pair energies ϕ ij (r ij ). Even if the potential is not defined uniquely, methodologies providing canonical averages will minimize errors on generated configurations, energies, and entropy comparatively to the full detailed ones. Three methodologies have been usually preferred for coarse-graining. (R) and g super atoms (R) are the radial distribution of the CG model at step k and of the real one, respectively. a is a relaxation constant. In practice, the inversion requires averaging the potential update ϕ • Boltzmann inversion (Reith et al., 2003) is a generalization of the previous method to polymers: the Helmholtz free energy F(R) −k B T ln g all atoms (R) is an initial guess for the potential g initial CG (R); the methodology is applied separately to intra-and intermolecular potentials (e.g., angular) with an order determined by the intensity of each potential (ϕ stretch Though the same methodologies could be applied to adjust forcefields to experimental data, they are presented in a logic of bottom-up modeling and simulation built on molecular mechanics. The knowledge gained at a small scale is used to reconstruct the macroscopic behavior. At small scales, forcefields are inherently repulsive, and the packing of particles is very compact. Excluded volumes vanish at a higher level of coarsegraining, and some overlap must be accepted to preserve the bulk density. At this stage, it is noteworthy that described procedures apply only to conservative forces. Dissipative forces required for DPD simulations must be adjusted to dynamic properties such as viscosity or self-diffusion. Finally, and compared with atomic force fields, coarse-grained forcefields depend on the system's state (solid, liquid, gas), possibly on temperature and local composition. They are, therefore, not recommended for the calculation of state diagrams without independent validation.

Thermodynamic Equilibrium Between Phases
Studying the effects of temperature and pressure on phase transition is an easy task with the generalization of efficient thermostats and barostats. Enforcing equilibrium between phases in molecular simulation is comparatively much more difficult. Small systems do not behave as large ones. Water did not freeze at the right temperature as it was confined. Water and oil will not create a nice interface within a small cubic simulation box. Droplets are not spherical, and their shape fluctuates with the number and the arrangement of the molecules inside each droplet. Besides, surface tension in droplets smaller than Tolman's length deviates from its macroscopic value. The finite size of tested systems complicates the construction of large droplets at high dilution. Setting a dimension longer than the two others forces the creation of one flat interface (or two). In this case, the transport of small solutes can be tracked between the two phases at reasonably low concentrations. However, the statistics are too poor to get reliable partition coefficients, and the strategy cannot be generalized to polymers.
Chemical equilibration is the chief difficulty in reaching thermodynamic equilibrium. It requires an exchange of molecules between the two phases. Liquid-vapor and liquidliquid equilibria can be efficiently studied in the Gibbs ensemble (Panagiotopoulos, 2002) for moderately dense multicomponent fluids. Figure 10 illustrates the principles. The phase coexistence properties of multicomponent systems are studied from a single Monte-Carlo simulation involving two distinct physical regions with different densities and compositions. The regions (boxes) can be dissociated without any common interface. Three types of perturbations are performed on each box: 1) a random displacement of molecules that ensures equilibrium within each region, 2) an equal and opposite change in the volume of the two regions that results in equality of pressures, and 3) random transfers of molecules that equalize the chemical potentials of each component in the two regions.
The great advantage of the Gibbs method over the alternative techniques for studying phase coexistence is that, in the Gibbs method, densities and compositions of the coexisting phases are resolved simultaneously. The coexistence line between phases does not require the determination of chemical potentials as a function of temperature and pressure. The list of operations must be picked from a repertoire of trial moves in a random order to preserve the microscopic reversibility. However, the Gibbs ensemble method remains impractical for very cohesive phases such as polymers and very ordered systems (crystals, liquidcrystals). Free-energy calculation methods must be preferred (see Figure 11).

Calculations of Activity Coefficients Within the Flory Approximation
Activity coefficients and excess chemical potentials are essential properties in food where local composition effects modify strongly how food constituents are transferred inside the food or loss during processing or storage. The same properties control how chemical contaminants are transferred to the food from the packaging or the storage environment. This section reports how the mean-field theory of Flory-Huggins formulated initially on lattices can be turned into a tailored and robust predictive method by using atomistic calculations as inputs. The benefit is twofold. The theory enables first-order predictions without explicitly representing the condensed phase (interacting liquids, entangled polymers). Atomistic calculations with detailed forcefields can replace obsolete group contribution methods and calculations on lattice symmetric imposing symmetric coordination numbers for A surrounding B, and B surrounding A.

Prior Definitions
G. N. Lewis proposed the concept of the thermodynamic potential of substance i in the phase α, μ i,α (T, P), while verifying that RT P the molar volume of the gas i. After integration at a constant temperature, the change in chemical potential from pressure P ref to P is RT ln P P ref for a pure ideal gas. Lewis generalized to phase mixtures (ideal or not) by replacing pressure by a function f i , called fugacity and assessing the capacity of a substance literally to flee: . Reference values at the temperature T μ ref i,α and f ref i,α are physically related, but the choice of the reference chemical potential or the reference fugacity is arbitrary. Fugacities of substances absorbed on solids can be determined directly by molecular modeling using the grandcanonical ensemble (μVT). Its utility lies in its resemblance with experimental measurements (e.g., sorption microbalance): the system interchanges particles with a reservoir so that the chemical potential can be by changing the number of particles changes accordingly. As for the Gibbs ensemble, the method is prohibitive for condensed phases. Flory-Huggins' theory offers a scale-invariant to estimate activity coefficients, c v i,α , relative to volume fractions ϕ i,α in any condensed phase, as summarized in Table 5. The presented relationships expressing c v i,α in function of ϕ i,α are equivalent to conventional sorption isotherms for binary and ternary mixtures. Indeed, activities and c v i,α are related as: FIGURE 10 | Strategies to represent implicitly thermodynamic equilibrium between phases in molecular simulations.
Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 23 By considering that partial molar volumes are similar in both phases, solute partition coefficients between the phases α and β are given by: Defining the activity coefficients respectively to volume fractions instead of molar fractions as in the regular solution theory offers several advantages. Volume fractions enable to consider situations where the solute i is larger than surrounding molecules and also the converse situation where i is infinitely smaller than polymer segments. At this stage, it is important to note that the reference state of solute i is an amorphous state. If its reference state is solid at temperature T, a fugacity correction equals to the ratio of fugacities between solid and liquid states must be added [see Ref. (Nguyen et al., 2017b)]. Flory-Huggins sorption isotherms are equivalent to Henry ones at infinite dilution and describe well swelling effects (Oswin type isotherms) at high concentrations. Sorption is assumed to be reversible (no hysteresis). Generalization to sorptiondesorption cycles with hysteresis when T g is crossed during sorption can be found in Ref. (Kadam et al., 2014). For liquids, the Kirkwood-Buff (KB) solution theory (Kirkwood and Buff, 1951) at a microscopic scale can provide exact macroscopic predictions even for interacting liquids and vice-versa, as shown by Arieh Ben-Naim (Ben-Naim, 1977). The pro and cons of KB and FH theories are discussed in (Nguyen et al., 2017b). Corrections derived from KB integrals and pair radial distributions are mainly required for compressible mixtures and large solutes, as shown in . The term n compressible k around i captures the positional entropy reduction due to the specific arrangement between a large solute i and the solvent molecules.

Calculation of χ T i,k at Atomistic Scale
The calculation procedure of excess chemical potentials is summarized in Figure 11, with technical details given in Refs. (Gillet et al., 2009;Gillet et al., 2010;Vitrac and Gillet, 2010;Nguyen et al., 2017a;Nguyen et al., 2017b;Zhu et al., 2019a). The Flory-Huggins isotherm comprises a positional entropy term and an enthalpic term proportional to the extent of {χ i,k } k P,F . The first term is determined by the number of molecules k, which are displaced by the introduction of i assuming that the mixture is incompressible (first approximation) or compressible (from radial distribution). The second term is essentially enthalpic, but it also contains the non-positional entropic contribution capturing the fluctuations of the contacts between molecules.
Both properties can be studied by molecular modeling (packaging of molecules, pair interactions energies). The advantage of the approach is the capacity to extrapolate the results to arbitrary binary, ternary and quaternary composition. Aqueous and hydroalcoholic solutions have been specifically discussed in Ref. . Flory-Huggins coefficients {χ i,k } k P,F are defined as the dimensionless mixing energy (enthalpy) in excess relative to pure compounds; it is defined for polymers (k P) and liquids (k F): where 〈X〉 T represents an ensemble-average of X, n P is the length used in the approximation of the polymer and n F 1.
In agreement with the original Flory approximation, enthalpies 〈h A+B 〉 T are estimated by summing energies of contact ϵ AB when B is contacting the seed molecule A. 〈h B+A 〉 T represents the same energy when B is used as a seed molecule. In practice, ϵ AB is calculated by choosing an orientation randomly for the contact molecule and by  (Zhu et al., 2019a)]. The formula applies only to homogeneous phases. For heterogeneous phases (e.g., emulsions), the effective partition coefficient/activity coefficient must be spatially homogenized.
Sorption of solute i in a liquid k = F or in a mixture of liquids F 1 and F 2

Binary mixtures
For binary mixtures, i + P (polymer), the activity coefficient is given by: i,F being the Flory-Huggins interaction coefficient between i and F. r i,k represents the number of molecules of F displaced by the insertion of the substance i in F. For most of the migrants, r −1 i,F is expected to be larger than unity in water and ethanol and lower than unity in oil. r i,F can be approximated as VF Vi with n compressible k around i accounting for the partial compressibility of the molecules of i and F. Rigorously, partial molar volumes should be used instead of molar volumes

Ternary mixtures
The activity coefficient in mixture of two miscible liquids F 1 and F 2 can be estimated with: with χ i,F1 ,F2 being a ternary Flory-Huggins coefficient whose contribution can be neglected in the absence of a specific ternary complex in the solution i,P being the Flory-Huggins interaction coefficient between i and F Tertiary mixture The activity coefficient in a wet polymer associated with a volume fraction of water ϕ w , depends on the three pairs of Flory-Huggins coefficients, χ i,P , χ w,P (water in dry P), χ i,w (solute in water), as The corresponding Flory-Huggins interaction coefficient in a copolymer AB, χ i,(AB) , reads: Eq. 51 averages of all possible contacts between i − A and i − B, where the term χ AB ϕ A ϕ B represents the additional cohesion energy brought about by the interactions between A and B. It can be generalized to more complex copolymers while the contacts can be assumed perfectly random (i.e., no macro-or microphase separation) Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 translating it along a random line until at least one point of contact is established between the van der Waals surfaces of the contact and seed molecules. The process is repeated for all conformers and stereoisomers considered. Finally, h A+B is estimated as the product of contact energies and the number of neighbors z AB (number of B molecules surrounding A): Equation 37 assumes that ϵ AB and z AB are statically independent (zero covariance). For polymers, the property of independence is achieved with a large enough n P value so that the surface of contact of the polymer is independent of the length of the considered polymer. The main advantage of the whole approach is that there is no need to represent entanglements in the polymer and free-volume. The shape of the backbone of an infinitely long chain with shorter oligomers prevents head and tail atoms from coming in contact with any van der Waals surface. Cooperative hydrogen bonding is captured by using a value of n cooperative greater than one. The latent heat of vaporization of water can be correctly approximated by 〈h water+water 〉 T using a value of n cooperative different to unity. This value depends on the type of forcefield used to simulate water. As an example, the rigid water model governed by the TIP4P forcefield gives an acceptable value with n cooperative 1 whereas n cooperative 4 is required with the same forcefield but using three-point charges (forcefield TIP3P). The number 4 reinforces in this case that any water molecule is on average involved in 4 hydrogen bonds of similar strength.
Contact energies are calculated irrespective of any temperature consideration. The effect of temperature is recovered by weighting the distribution of contact energies with the Boltzmann factor B(ϵ): At the expense of calculating two integrals, Eq. 38 can be used to estimate χ (n k ,T) i,k at several temperatures. Conformers need to be generated so that they represent their conformations (radial distribution, constraints of torsion) in the corresponding condensed phase. In practice, they are sampled from molecular dynamics simulations of the equivalent condensed phase.

Examples
The whole approach is versatile and can cover a large number of situations, as exemplified in Figure 12. The study detailed in Ref. (Lu et al., 2008) shows that partition coefficients between a random copolymer material and water are determined for a broad range of organic solutes and isomers. Conventional group methods, such as solubility parameters, drop important molecular details such as structural isomerism, stereoisomerism, tacticity, pattern order, and local chain composition. The examples from Ref. (Nguyen et al., 2017b) show how the calculations can be done inexpensively from homopolymers and then subsequently extended to copolymers ( Figure 12A) via Eq. 51 or calculated directly on copolymers ( Figure 12B). The results are similar For genuinely random copolymers, and Eq. 51 averages all microscopic contacts, including specific hydrogen bonding, satisfactorily. Nevertheless, mixtures of homopolymers must be considered differently by averaging the energy of mixing and not contact energies, as discussed in the supplementary information of Ref. (Nguyen et al., 2017b).
By combining atomistic calculations with Eq. 35, partition coefficients were calculated continuously as a function of the polymer's acetylation rate. Acetate groups were assumed randomly distributed along the chain. Confidence intervals represent the uncertainty in the convergence of averages with minimum covariances ( Figure 12C). Extensive sampling of pair configurations (from several hundred million to several hundred billion) facilitates identifying specific interactions between chemicals. Boltzmann-weighted ensemble averaging (Canonical ensemble) makes it possible to extrapolate sampled results to different temperatures. The only condition is that the conformations of tested conformers remain representative of the state of each component in the mixture. For production, results can be tabulated (sampled pair contact energies, coordination numbers, specific volumes, etc.) and stored for future use. Figure 12D shows the prototype of an in-house software used to calculate partition coefficients between arbitrary plastic materials and food simulants. The calculation procedure is predictive for any organic solute without net charge and prediction errors in the same range as experimental errors. The applicability to food constituents has been tested successfully for starch-water mixtures at rubber state, as met during deepfrying (Nguyen et al., 2017a). The concept of free volume is fundamental to liquids and amorphous polymers. Free volumes are defined as the free space between molecules or chain segments; they control important properties of amorphous phases such as their glass transition temperature, their thermal expansion, and their transport properties (White and Lipson, 2016). The law of Stokes-Einstein predicts diffusion coefficients in fluid phases, which are proportional to 1) temperature, 2) the reciprocal viscosity of the liquid, and 3) the reciprocal diameter of the solute. The hydrodynamics theory behind the Stokes-Einstein relationship is, however, not sufficient to illuminate the activation of diffusion by temperature. Indeed, the fluidity (the reciprocal viscosity) or equivalently the friction coefficient is also activated by temperature. Besides, the theory falls short of explaining the extreme mass dependence of diffusivities as observed in amorphous polymers (gels, thermoplastics, thermosets). Although diffusion at constant volume requires a little activation, k B T, the cooperation of many degrees of freedom in the host environment is needed to induce the translation of gases or organic solutes larger than surrounding rigid segments. Thermal expansion of van-der-Walls liquids and rubber polymers (ca. 0.1% every 10°C in polyethylene) is the dominant cause of the acceleration of the diffusion in these phases and not directly temperature. Figure 13 justifies that molecular modeling alone cannot address today all aspects of the mechanism of diffusion for molecules larger than voids in cohesive systems (glasses, entangled polymers). There are two limitations. Firstly, the relaxation modes of entangled polymers or liquids at glassy state cannot be simulated at atomistic scale, as decorrelation times are not accessible to brute force calculations before the next decade. Secondly, large molecules do not translate as a single rigid unit but proceed by small reorientations of internal subunits. A macro-state jump of the center-of-mass succeeds on distances larger than the gyration radius of the diffusant involves, therefore, a large collection of micro-state reorientations. We thought falsely that coarse-graining applied to both the diffusant and the surroundings could alleviate previous time scale limitations. But, because the coarse-graining process removes important atomistic details governing the rate of translation of real diffusants, it created artifacts instead. Logterm molecular dynamics up to 0.06 ms did not succeed in reproducing the strong mass dependence observed at solidstate, as discussed in Ref. (Durand et al., 2010). Heterogeneous dynamics were replaced by smoother ones consistently with Rouse modes observed at molten state when the entire chain of the polymer can translate instead of small segments at solid-state. For example, the displacements of hydrogen atoms along the polyethylene chain control the diffusion rate of helium. In contrast, the shaft motions of connected rigid monomers (-CH 2 CH 2 -) govern the translation of methane (CH 4 ). Larger molecules require activating the relaxation of several polyethylene segments, causing hops to span over a broad range of time scales, as previously discussed in Section 3.2.2.

Theory of Vrentas and Duda and its Extensions
Free-volume theory and its extensions, shown in Figure 14, offer an elegant framework to address the problem of many FIGURE 13 | Time scale limitations met by atomistic simulations: (A) comparison between macroscopic fluctuations and those accessible today to molecular modeling; (B) projection of capacity of calculations for the next decade in polymers; (C) scaling of diffusion coefficients of gases and organic molecules with their van-der-Waals volume V vdw in rubber (natural rubber) and glassy polymers (polyvinyl chloride: PVC) at 298 K (Berens and Hopfenberg, 1982).
Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 28 correlations between the translation of the center-of-mass of the solute and the displacements of the surrounding molecules or segments theoretically. Using theory instead of brute force calculations removes restrictions on solute size and the distance to the glass transition temperature T g . The mean-field approximation accepts a rapid parameterization from simple geometry calculations on conformers and solute self-diffusion coefficients determined experimentally or by molecular modeling. By assuming that the translation of a solute is only possible if the free volume fluctuation is larger than a critical volumeV p , D is obtained by superimposing all the translation modes for V ≥V p , one gets FIGURE 14 | Principles of free-volume theory. (A) Interpretation of free volumes in liquids/amorphous polymers at equilibrium (rubber state) and subjected to various degrees of subcooling (glassy state); (B) Connection between the original theory of hole-free volume of diffusion (Cohen and Turnbull, 1959;Vrentas and Duda, 1977b;Vrentas and Vrentas, 1994a;Vrentas and Vrentas, 1994b) and its recent extensions (Fang et al., 2013;Fang and Vitrac, 2017;Zhu et al., 2019b) The idea behind Eq. 39 initially postulated by Cohen and Turbull (Cohen and Turnbull, 1959) is that a diffusivity may exist for each type of free volume. We could imagine that all molecules do not translate at the same time to the same fluctuations of free volumes. In the original theory, volumes are expressed in excess of an irreducible volume at a critical temperature T cr : where α(T) is the coefficient of thermal expansion of the fluid. This hypothesis is not realistic in the case of polymers, because the redistribution of interstitial free volumes is associated with relaxation of the polymer chain and thus with very high activation energies (Watanabe, 1999). Vrentas and Duda (Vrentas and Duda, 1977a;Vrentas and Duda, 1977b) removed the objection by arguing that only one hypothetical fraction of free volumes contributed to diffusion; they called it hole-free volumes (hFV). Such theoretical volumes would be redistributed within the polymer without a specific energy barrier. The authors consolidated their theory and arguments during the following 2 decades (Vrentas and Vrentas, 1998). The final theory connects self-diffusion coefficients with mutual and trace diffusion coefficients in liquids and polymers at glassy and rubber states. The canonical formulation for a rigid solute in a rubber polymer at infinite dilution reads [see Ref. (Vrentas and Vrentas, 1998) for details]: where ξV * 2 is the critical hFV to induce a translation; V p FH 2 is the critical hFV available for the solute (1) and the polymer (2). ξ is the ratio between the critical volumes of the jumping units of the solute and the polymer.  (Fang et al., 2013) and (Zhu et al., 2019b)]; (B) an example of predictions for water in polyethylene terephthalate between 5°C and 190°C. The only degree of freedom is the ratio of thermal expansion between glassy and rubber states, r, and characterizing the degree of subcooling of the processed polymer (see Figure 14A). It needs to be fitted or measured independently on the finished product.
Frontiers in Chemical Engineering | www.frontiersin.org January 2022 | Volume 3 | Article 786879 Because all the parameters can be inferred from quantities independent of the diffusion coefficients or simulated by molecular modeling, it is a predictive and consistent theory from first principles and with Boltzmann statistics. The theory has been later extended (Fang et al., 2013;Fang and Vitrac, 2017;Zhu et al., 2019b) to cover flexible molecules and not only rigid solutes. Earlier mechanisms were kept, a rigid segment can translate or reorient when a void of matching the size of the rigid pattern (blob) opens during redistribution (Fang et al., 2013). The dynamics of the center of mass corresponding to several connected segments is particularly complex as all blobs cannot translate at once. This mechanism was studied by molecular dynamics and mapped onto quantum oscillators to correct the probability of each microstate. The new description predicted successfully that organic solutes (aromatic, aliphatic) in solid polymers would translate with all rigid segments blocked except one. With the help of mode coupled theory, all hFVrelated effects were recast in the behavior of linear probes and more specifically in the scaling exponent, α lin relating D with the number of connected blobs, as D ∝ M −α lin (Zhu et al., 2019b): The revised hFV theory is equipped to capture diffusivities in rubber and glassy polymers of any solute consisting of one single blob or of many blobs of the same size or not, symmetric or not. Figure 15 provides the practical formula and scaling curves valid for 12 polymers (representing the following polymer families: polyolefins, polyesters, polyvinyls, and polyamides). An example of predictions is shown for water in polyethylene terephthalate (PET). The results are remarkable as PET is a polymer, which represents one-third of all food packaging materials; it is the only polymer that can be recycled for food contact so far, and no freevolume theory has been previously published for this polymer. Calculations of the shelf-life of beverages based on diffusion coefficients and sorption isotherms are presented in Ref. (Zhu et al., 2019c).

GUIDELINES FROM A MICROSCOPIC PERSPECTIVE
Diffusive or dispersive phenomena have an entropic origin; they are therefore everywhere and characterize the propensity of matter to separate under the effect of thermal agitation. In foods, they are affected by food structure, free volumes, and molecular interactions. The concepts of entropy (thermal, configurational, conformational) are introduced frugally in continuum mechanics via the selection of constitutive laws next to the conservation laws. These constitutive laws are first or second-order approximations of more complex microscopic descriptions, whose complexity is generally reduced to an apparent diffusivity. In the context of food and food packaging, effective diffusivities can range from 10 −3 to 10 −19 m 2 ·s −1 . This section reviews the strategies for estimating them by adopting a microscopic and mechanistic perspective. The first guidance uses a classification of diffusive phenomena to choose between spatial and temporal homogenization and to suggest a technique. Most concepts are flexible and can be adapted to the phases, domains, species, scales, dimensions, and coupling involved. The second guidance illustrates how apparent radial forces as real as osmotic forces can emerge from microscopic descriptions to describe large classes of entropy-driven problems. This introductive example offers an invitation to discover a highly active research field aiming to understand and control microscopic phenomena.

Mechanism-Based Strategies to Use and Derive Effective Diffusivities
Analyzing transport properties in foods at all scales at once is either intractable or useless due to the inherent variability of the biological matter. By contrast, replacing microscopic and molecular composite effects with effective diffusivities for relevant geometry domains (skin, crust, core, packaging layer, etc.) allows the use of generic multiphysics simulation software. They are not developed for food products and processes and use continuum mechanics to resolve coupled transfer equations. By repeating the simulations for different property values corresponding to different compositional and structural effects, macroscopic properties of processed foods (e.g., barrier or shelf-life performances) can be related to microscopic or molecular ones and possibly adjusted to reach the desired target. When the properties cannot be tabulated in advance, the effective properties need to be updated regularly. Parallel computing and client-server architecture (usually available in commercial packages) facilitate the development of hierarchical modeling between scales. A whole framework so-called FMECAengine 3D has been developed by us with combine three main loops: one on microscopic properties so-called [E]valuation, a second one [D]ecision devoted to 3D at the macroscopic scale, and a final one [R]esolve relying on multicriteria optimization with microscopic properties as free parameters and product characteristic as goals [see (Zhu et al., 2019c)].
The choice of the appropriate local and global levels, where the equivalent diffusivity is calculated and valid, respectively, remains a challenge. Intuitively, one might think that the smallest level would require the least computational effort (e.g., fastest mixing or decorrelation), but this does not consider the cost of refined details. For most problems, temporal and spatial problems are very intertwined. The chief question is which variable needs absolutely to be considered continuous. As a rule of thumb, it is preferable to consider continuous in-time representations when energy barriers control long waiting time distributions.
The following example will show that considering time and space continuous requires introducing proper probability measures. Indeed, what is the probability of two particles moving randomly to meet in diluted regime over a given duration? When a collection of particles are considered, the many interactions problem (collisions are very likely) removes the need for detailed space-time statistics as soon as the density and packaging are well reproduced. Overview in Section 3.2.1, Section 4.2, Figure 3, Figure 8 •      As general guidance, we propose a balance by separating problems above and below the thermodynamic limit (see Figure 1B). Table 6 lists practical strategies to extract effective diffusivities from poor selective mechanisms such as viscous flow, bulk, and Knudsen, surface diffusion to very selective ones such as molecular sieving, capillary condensation, diffusion solubility. Above the thermodynamic limit, diffusion phenomena are self-similar, but scaling remains particularly difficult to predict in the presence of multiple domains, components, or coupling with other driving forces. Below the thermodynamic limit, the many interactions with components, interfaces, and phases may delay the convergence dramatically to the thermodynamic limit and lower diffusivities.

An Analytical Example of How to Treat Microscopic Mechanisms of Diffusion With Radial Forces
Spatial heterogeneities and viscoelastic relaxation are examples preventing the time and length scales from being separable. They appear in various diffusive phenomena ranging from dissolution to bubble growth. At the microscopic scale, relating the diffusive behavior to the composition gradient as driving forces becomes rapidly impractical. Concentration fields frequently fluctuate due to morphological details (e.g., obstacles, disorder, correlation, network dead-ends), thermodynamic effects (supersaturation, sorption, crystallization . . . ) and external coupling (osmotic effects, swelling, shear, flow, Marangoni convection). All dispersive phenomena share a common feature: the position and velocity of individual particles cannot be known for sure. The behavior of a single particle can be described in terms of probability density functions in the laboratory-fixed coordinates (via the Fokker-Planck equation) or in a Lagrangian reference frame (via the Langevin equation). However, some concepts are more adapted than others to describe certain phenomena such as order-disorder transitions (mixing, melting, solubilization, elastic behavior). The concept of radial force is one of them and can be deduced from the representation of entropy in statistical mechanics (Neumann, 1980;Roos, 2014).
In the food context, radial forces can be naturally derived and from the differentiation of the chemical potential of one single particle (atom, blob, molecule, fluid volume, colloidal particle), denoted μ Δr 1,r , by assuming that it is located within a spherical shell of radius r and thickness Δr. For an eventually non-ideal mixture [see the discussion in Ref. (Krishna and Wesselingh, 1997)] and Eq. 39 in Ref. (Zhu et al., 2019a), the magnitude of the radial force reads for each particle: where ρ Δr 1,r 1/(4πr 2 Δr) is the density of one single particle in the shell and Γ 1 1 + z ln c v 1 z ln ϕ 1 | T is the thermodynamic factor associated with scaling the activity coefficient c v 1 with the volume fraction ϕ 1 .
From Eq. 43, all displacements of species due to chemical potentials (without external forces) can be associated with a radial and centrifugal force with an average magnitude proportional to the reciprocal distance between the source and recipient regions. Conventional descriptions of random walks are recovered by equilibrating f r with the viscous drag ζ 1 d〈r〉 dt associated with the friction coefficient ζ 1 k B T D1 , where D 1 is the local diffusivity. Integrating the force balance in spherical coordinates leads to the following differential equation: Equation 44 is locally valid for any particle following a random walk in a medium governed by D 1 and Γ 1 . Its integration can be carried out either using Monte-Carlo sampling or with analytical solutions of the probability density function p r , which gives the probability that a particle colliding with its environment reaches a distance between r and r + dr during t. For an isotropic, uniform, and infinite medium, the distribution has a Gaussian form: where σ is the mode of the distribution (likeliest value of r), 〈r〉 2 π √ σ, 〈r 2 〉 3 2 σ 2 , 〈 1 r 〉 2 π √ 1 σ . Eq. 44 becomes: with the unique solution σ 2 4ΓD 1 t or equivalently 〈r 2 〉 3 2 σ 2 6ΓD 1 t, which is the same as the Einstein equation for a free random walk for Γ 1.
The approach is very general and can be applied to any particle and pair of particles (gel junction, particles representing an interface or an ordered phase) with the addition of proper forces (Coulombic, surface tension, rubber elasticity). In the general case, the distribution is unknown, but it is the most probable distribution maximizing entropy (or the amount of missing information). It can be derived by invoking the central limit theorem or from more straightforward rules. For any finite region divided into cells, the probability of reaching each cell on the long term is proportional to the cell volume divided by the local activity coefficient with a normalization, which depends on the distance to the neighboring regions as shown in Ref. (Vitrac and Hayert, 2020).

CONCLUSION AND PERSPECTIVES
This article tries to propose some generic recipes for those who want to get into the "mathematical cooking" of food products. The core of the methods consists in presenting how chemical structural information, cell structure and possibly phase organization can be encoded in a coherent mathematical scheme that can be used with classical simulation software at the scale of continuous media. The main tool is the asymptotic analysis. Depending on the importance of the thermodynamic contribution, it must be applied to space or time dimensions. The available computing power necessarily limits the possibilities of prediction and extrapolation but the bricks of a rigorous construction are illustrated. They already offer a significant progress in a context where few properties have been tabulated on raw materials and foods during processing. Indeed modeling and simulation in food is a generally overused term because it means both a descriptive mathematical model and a mechanistic model that can be used to extrapolate, generalize, or even imagine new food products. But the properties are missing at fine scale and the goals can be multifaceted (nutritional, sensorial, cooking behavior, shelf life, health, environmental, etc.). Cooking, drying, freezing and packing food on a computer does not seem very appealing, but the ingredients are there this time. The substances, the macromolecules, the processing conditions can be included in the multi-scale modeling scheme. Not everything can be solved yet, but great steps have been taken in recent years: interactions with models from other disciplines, continuous increase of computing power, many open-source calculation codes, new high-throughput 3D reconstruction tools, more food examples, a better acceptance of model predictions by authorities. The framework is there; now, we must imagine what we can do with the available tools and imagine the new tools we would need for food. We miss quantum methods for foods; they would be beneficial for understanding the complex biochemical reactions that occur during food processing: production of allergens, toxic compounds, depletion of compounds of nutritional interest, digestion, and toxicology of food constituents. This article was written during the Covid-19 pandemic; molecular modeling has been proposed and is in the process of understanding how this new virus works. These new models give us hope that food construction and deconstruction rules can be reviewed very soon.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.