Modeling food matrix effects on chemical reactivity: Challenges and perspectives

ABSTRACT

The same chemical reaction may be different in terms of its position of the equilibrium (i.e., thermodynamics) and its kinetics when studied in different foods. The diversity in the chemical composition of food and in its structural organization at macro-, meso-, and microscopic levels, that is, the food matrix, is responsible for this difference. In this viewpoint paper, the multiple, and interconnected ways the food matrix can affect chemical reactivity are summarized. Moreover, mechanistic and empirical approaches to explain and predict the effect of food matrix on chemical reactivity are described. Mechanistic models aim to quantify the effect of food matrix based on a detailed understanding of the chemical and physical phenomena occurring in food. Their applicability is limited at the moment to very simple food systems. Empirical modeling based on machine learning combined with data-mining techniques may represent an alternative, useful option to predict the effect of the food matrix on chemical reactivity and to identify chemical and physical properties to be further tested. In such a way the mechanistic understanding of the effect of the food matrix on chemical reactions can be improved.

KEYWORDS:

• Food matrix
• kinetics
• thermodynamics
• modeling
• activity coefficients
• fingerprinting

Introduction

At its most basic level, food quality depends on how compounds in the food interact, via chemical and physical reactions. Such reactions are subject to basic laws from chemistry and physics, irrespective from where they occur, and they can be fully characterized in terms of thermodynamics and kinetics. However, when a chemical reaction occurs in real foods, in many cases different results are obtained as compared to the same reaction occurring in simple aqueous solutions. This happens because foods do not behave like ideal solutions. Indeed, foods are multicomponent, concentrated systems, consisting of water-insoluble materials (for instance, cell membranes and lipid droplets), complicated aqueous solutions of ionic and nonionic compounds, with various phases that can also be in different physical states. Even more importantly, however, results are likely to be different when the same reaction occurs in different foods: the food matrix changes the rate of the reaction and/or its equilibrium position. The food matrix can be defined as the whole of the chemical components of food and their molecular relationships, the chemical composition of food, and the way those components are structurally organized at micro-, meso-, and macroscopic scales. Recognizing the effect of the food matrix on chemical and enzymatic reactions means to recognize that, even starting from the same initial conditions and applying the same physical conditions to the surroundings, the same reaction would not evolve in the same way in ideal solutions as in real foods or, even more importantly, in different foods. This is equivalent to say that differences among food matrices are responsible for the apparent variability that is observed in a certain chemical behavior (for instance, the stability of a certain nutrient upon thermal treatment) within a “population” of food products.

The capacity of quantitatively describing and predicting the effect of the food matrix on chemical reactivity would be of paramount importance for food scientists because it would lead to a more rational approach to food product development and handling. Unfortunately, a fully mechanistic description of the effect of the food matrix on chemical reactivity is difficult: foods are very complex systems. Here we suggest that an empirical approach based on chemical/physical fingerprinting combined with data mining could represent a complementary strategy to predict this effect as well as to generate hypotheses to be tested with a more fundamental approach. The aim of this paper is to provide tools for understanding and predicting the effect of the food matrix on chemical reactions in food and as such, for controlling the final quality attributes of food. We will first provide two examples of the variability caused by the food matrix on thermal stability of important bioactive compounds, namely glucosinolates and vitamin C which reflects itself in different rate constants for the same reaction in different foods. Then, we will describe two approaches for the description and prediction of the food matrix effect on chemical reactivity: (1) a mechanistic approach based upon fundamental laws of chemistry and physics and (2) an empirical approach based on machine learning and data mining.

Effect of the food matrix on thermal stability of bioactive compounds: Glucosinolates

Glucosinolates (GLs, β-thioglycoside N-hydroxysulphates linked with a sulphur β-d-glucopyranose) are secondary metabolites that occur in Brassica vegetables (broccoli, cauliflower, cabbage, etc.) (Fahey et al., 2001). GLs are hydrolyzed by a group of endogenous β-glucosidases termed myrosinase (MYR). MYR is stored separately from glucosinolates in plants, but will come into contact with GLs upon cell damage. The hydrolysis of GLs produces isothiocyanates, nitriles, epithionitriles, or thiocyanates depending on certain conditions, such as pH, vitamin C content and the presence of epithiospecifier proteins (Hanschen et al., 2017). Isothiocyanates have been demonstrated to be potentially health-beneficial for humans because of their anticarcinogenic properties (Steinmetz, 1996; Traka, 2009; Miller, 2012; Biswa Nath Das, 2013).

During the past 20 years, much research has been conducted on the effect of postharvest treatments on nonenzymatic GLs thermal degradation (Fahey et al., 2001). GLs thermal degradation has been mainly described by first-order kinetics. In Figure 1, a graphical depiction of the variability in rate constants reported in the scientific literature for GLs degradation is provided. Figure 1 only includes a selection of rate constants from literature. A more comprehensive overview of the kinetic parameters and of the food matrices in which they have been derived is available in the Appendix. From Figure 1 and the Appendix, it emerges that the thermal breakdown of the same GL shows different rate constants in the different Brassica vegetables. For instance, heating ground samples of red cabbage, broccoli, Brussels sprouts, pak choi, and Chinese cabbage (previously microwaved and blended in liquid nitrogen to rule out MYR activity) at 100 ˚C resulted in different degradation rate constants for three GLs. Brussels sprouts were the vegetables in which the highest degradation rate constants for all three GLs were observed (Dekker et al., 2009). The authors hypothesized that the different chemical environments of the vegetables, such as differences in pH and cellular composition, can promote or inhibit GLs thermal degradation. In another study, thermal degradation of four aliphatic and four indolic GLs during boiling, steaming, and stir-frying of three Brassica oleracea varieties, the Dutch “Boerenkool” (kale), the Italian “Broccolo Lavagnino,” and “Nero di Toscana” were compared (Giambanelli et al., 2016). The thermal degradation rate constants were estimated by second-order kinetics, and the vegetables were only chopped (to a similar size). The results showed that for all the cooking and for all the GLs, the highest GLs degradation rates occurred in the “Nero di Toscana” variety. Such matrix effect can be detected also in different genetic lines of the same variety. Hennig and others (2013) investigated the nonenzymatic GLs thermal degradation in Brassica oleracea doubled haploid populations (Hennig et al., 2013). The authors hypothesized that since the metabolic composition is determined by biochemical traits that are (partly) genetically regulated, the thermal degradation of GLs is genetically regulated as well and, therefore, genomic regions can be identified that influence the thermal degradation rates of GLs. GLs degradation was modeled using first-order kinetics and the rate constants were determined for five aliphatic and two indolic GLs. A threefold difference in GL degradation rate constants among the doubled haploid population was observed for most of the GLs, whereas degradation rate constants for glucoraphanin differed by ninefold. The nonenzymatic thermal degradation of GLs was also investigated in broccoli with different a_w (Oliviero et al., 2012). It was shown that a drier matrix increased the thermal stability of GLs up to 100°C but not at 120°C. GLs degradation has been also studied in dilute aqueous solutions at different pH condition, water activity a_w, and in the presence of Fe (II) and vitamin C (Hanschen et al., 2012). In that study it was shown that the addition of broccoli sprout powder to an aqueous model system significantly reduced thermal stability of aliphatic GLs. In line with this observation, in a more recent study it was shown that “dilution” of the broccoli matrix with ground potato, corn starch, lentil protein, or onion increases the thermal stability of GLs, in particular in broccoli-onion mixtures (Giambanelli et al., 2015). The authors hypothesized that the protective effect may have been caused by flavonoids of which onions are among the richest sources (Slimestad et al., 2007), but could not find a correlation between different flavonoid concentrations (from three onion varieties) and GLs thermal stability.

Figure 1. Selection of rate constants for thermal degradation of glucosinolates published in the scientific literature. Measures of uncertainties have not been reported for clarity of representation. RSD% was in no case > 10%. More information on the food matrices where the rate constants have been derived can be found in the appendix.

Effect of the food matrix on thermal stability of micronutrients: vitamin C

Vitamin C is composed of L-ascorbic acid and its oxidized form, L-dehydroascorbic acid. Vitamin C is not synthesized in humans and its deficiency causes a disease known as scurvy. The main dietary sources of vitamin C are citrus fruits, Brassica vegetables, and many other fruits and vegetables including melons, kiwi, and bell peppers. Storage and thermal treatments may drastically reduce vitamin C content in foods. The main degradation pathway for ascorbic acid is thought to be its oxidation to dehydroascorbic acid, followed by its irreversible hydrolysis to 2,3 diketogulonic acid, which possesses no vitamin C activity (Giannakourou and Taoukis, 2003). However, a nonoxidative pathway has been also described where ascorbic acid is degraded without producing dehydroascorbic acid (Peleg et al., 2016). A marked effect of the food matrix on vitamin C thermal stability is recognizable, as reported in Table 1. For instance, the degradation rate constants of vitamin C during storage were reported to be different among several citrus fruits of the same genus (Burdurlu, 2006). The authors stored grapefruit, tangerine, lemon, and orange juice at 28, 37, and 45 ˚C and modeled the vitamin C degradation rates by first-order kinetics. At 28 ˚C the reaction rate constant was higher in lemon juice than in orange juice, whereas at 45 ˚C the opposite was observed. The degradation rate constants of vitamin C in grapefruit and tangerine juices were similar and intermediate to those of lemon and orange juices. Stešková et al. (2006) reviewed the thermal stability of vitamin C in fortified food during storage, modeling the degradation rate by first-order kinetics (Stešková, 2006). The fortified foods investigated were milk, evaporated milk, several types of bread, bran and potato flakes, cereals, cocoa powder, and apple chips. At the same storage temperature (25 ˚C), the highest degradation rate constants were observed in the breads, whereas in all the other foods, reaction rate constants were much lower, except for the fortified milk that showed intermediate values. The authors hypothesized that the difference in water content between the breads and the other dry foods (such as flakes and cocoa powder) can explain the different degradation rates of vitamin C. However, the fortified milk had higher water content than bread which showed the highest degradation rate constants. The same authors also reviewed the degradation of vitamin C in fortified beverages during storage. The beverages investigated were apple, cranberry, grapefruit, pineapple, tomato, and vegetable juices, but also grape and orange drinks, dry fruit drink mix and a carbonated beverage. The carbonated beverage showed highest degradation rates and the dry fruit drink mix the lowest (Stešková, 2006). Other studies have reported that vitamin C is more stable in dry matrices or in matrices with lower a_w (Laing et al., 1978; Uddin et al., 2001) and its dependence on a_w has been often used to optimize the air drying process to retain the vitamin C content ((Mishkin et al., 1984; Mishkin, 2006; Karim., 2009; Jin et al., 2014; Jin, 2014). Vitamin C stability during frozen storage was found in the order okra > green peas = green beans > spinach (Giannakourou and Taoukis, 2003). The different thermal stability of vitamin C in different foods can be attributed to multiple factors, for example, the different tissue structure, mechanical damage during harvesting, intrinsic enzyme (ascorbate oxidase) and sulfhydryl group content, and the presence of metal ions, such as Fe³⁺ and Cu²⁺, which act as catalysts for degradation. Last but not the least, the level of oxygen in food is also crucial since it acts as reactant in vitamin C oxidation.

Table 1. Overview of kinetic parameters reported in the scientific literature for the thermal degradation of vitamin C in different food matrices.

Compound	Matrix	Reaction order	Reaction rate constant	Activation energy	Reference
Vitamin C	Meet based model system with different aw	First order	k92˚C (min^–1)	kJ mol^−1	(Laing et al., 1978)
	0.7		7.9	58.6
	0.8		8.5	71.1
	0.9		10.5	66.9
Vitamin C	Vegetables	First order	k-20˚C × 10^–3 (d^–1)	kJ mol^−1 ± SD	(Giannakourou and Taoukis, 2003)
	Green peas		2.13	98 ± 10
	Spinach		4.54	112 ± 23
	Green beans		2.23	101
	Okra		1.05	106
Vitamin C	Kiwifruits at different aw	First order	k40˚C × 10^–3 (min^−1)	kJ mol^−1	(Uddin et al., 2001)
	0.51		0.05	38
	0.62		0.15	43
	0.75		0.50	52
	0.83		0.74	62
Vitamin C	Fruit juices	Zero-order reaction	k40˚C (min^−1)	kJ mol^−1	(Ozkan et al., 2004)
	Orange juice		0.06	40
	Sour cherry nectar		0.23	28
	Grape juice		0.16	26
	Pomegranate juice		0.04	71
Vitamin C	Soft Drinks	First order	k23˚C (min^−1)	Not estimated	(Stešková, 2006)
	Dry fruit drink mix		11.7
	Apple juice		73.0
	Cranberry juice		39.8
	Grapefruit juice		39.8
	Pineapple juice		46.5
	Tomato juice		41.5
	Vegetable juice		73.0
	Grape drink		51.5
	Orange drink		41.5
	Carbonated beverages		96.3
Vitamin C	Juices	First order	k45˚C ± SD (× 10^2, min^−1)	kJ mol^−1± SD	(Burdurlu, 2006)
	Orange		0.41 ± 0.05	105.4 ± 18.0
	Lemon		0.34 ± 0.05	53.5 ± 4.2
	Grapefruit		0.40 ± 0.03	77.0 ± 4.6
	Tangerine		0.41 ± 0.05	79.0 ± 2.9

Vitamin C	Fruits	Two-fraction first order	k1 at 100˚C (min^−1)	k2 at 100˚C (min^−1)	Ea1 kJ/mol^−1 ± SD	Ea2 kJ mol^-1 ± SD	(Hiwilepo-van Hal et al., 2012)
	Marula		7.2 ± 2.1 (× 10^−3)	7.9 ± 5.6 (× 10^−4)	29 ± 16	119 ± 26
	Guava		1.2 ± 0.3 (× 10^−1)	0.53 ± 1.8 (× 10^−5)	58 ±10	190 ± 89
	Mango		1.3 ± 0.5 (× 10^−1)	no reliable estimation obtained	39 ± 14	no reliable estimation obtained

From Figure 1 and Table 1 it is clear that kinetic parameters for the same reaction may differ among different food matrices (e.g., different juices or juices at different a_w) and those differences depend on the particular chemical milieu wherein the reaction occurs. In the following sections we will provide some basic thermodynamic and kinetic details to understand, at least qualitatively, the effect of the food matrix and we will present two approaches to quantitatively describe and predict such effect.

Thermodynamics and kinetics of (bio)chemical reactions

Chemical reactions can be characterized in terms of their thermodynamic properties determining the position of the equilibrium (i.e., by equilibrium constants) and kinetic properties determining the rate at which equilibrium is attained (i.e., by rate constants and their temperature- and pressure-dependence). Let’s consider, for the sake of simplicity, a reversible elementary bimolecular reaction:(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1)

When equilibrium is attained the following expression holds:(2) $K_{th} = \frac{a_C a_D}{a_A a_B}$(2) where $K_{th}$ is the thermodynamic equilibrium constant for the reaction, $a_C,a_D,a_A$, and $a_B$ are the thermodynamic activities of products C and D and reactants A and B, respectively. Thermodynamic activities are related to molar concentrations through the following relation:(3) $a = y*c$(3) where $a$ is the thermodynamic activity, $y$ is the molar activity coefficient and $c$ is the concentration of the chemical species in mol L⁻¹. Under ideal conditions (gas state or very dilute solutions)$y$ = 1 and thermodynamic activities equal concentrations. From eq. 2(2) $K_{th} = \frac{a_C a_D}{a_A a_B}$(2) and 3(3) $a = y*c$(3) , it follows that, in nonideal conditions, the following expression holds for experimentally observed concentrations:(4) $K_{obs} = K_{th} \frac{y_{C }{y}_D}{y_{A }{y}_B}$(4)

In nonideal conditions, the observed equilibrium constant $K_{obs}$ is different from the one that would be measured under ideal conditions $K_{th}$, because activities do not coincide with concentrations ($y$ values are ≠ 1). The dependency of the equilibrium constant on temperature is expressed by the van ‘t Hoff equation:(5) $\frac{d\ln K_{th}}{d(1/T)} = - \frac{\Delta H_r^{\circ }}{R}$(5) where $T$ is the absolute temperature, $\Delta H_r^{\circ }$ is the standard enthalpy change of the reaction and $R$ is the gas constant.

Transition state theory has been developed to describe the effect of temperature on rate constants. According to the transition state theory, reactants are in quasi-equilibrium with the so-called transition state (or activated complex). The rate of the forward reaction in eq. 1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) depends on the rate at which this transition state converts into the reaction products. This is expressed as follows:(6) $A + B \stackrel{K_f{}^{‡}}{\iff } TS \stackrel{k_{id, f}}{\to } C + D$(6)

The same holds for the reverse reaction:(7) $C +D \stackrel{K_r{}^{‡}}{\iff } TS \stackrel{k_{id, r}}{\to } A + B$(7) where $TS$ is the transition state, $k_{id,f}$and $k_{id,r}$ are the rate constants for the conversion of TS to C + D and to A + B, respectively. $K_f{}^{‡}$ is the equilibrium constant for the formation of the TS from reactants:(8) $K_f{}^{‡} = \frac{a_{TS}}{a_A a_B}$(8) where $a_{TS}$, $a_A$, and $a_B$ are the thermodynamic activities of $TS$, A, and B, respectively.

$K_r{}^{‡}$ is the equilibrium constant for the formation of the TS from products:(9) $K_r{}^{‡} = \frac{a_{TS}}{a_C a_D}$(9) where $a_{TS}$, $a_A$, and $a_B$ are the thermodynamic activities of $TS$, A, and B, respectively.

Under the representation of eq. (6(6) $A + B \stackrel{K_f{}^{‡}}{\iff } TS \stackrel{k_{id, f}}{\to } C + D$(6) ) and (7(7) $C +D \stackrel{K_r{}^{‡}}{\iff } TS \stackrel{k_{id, r}}{\to } A + B$(7) ), based on the law of mass action, the rate of the forward reaction in eq. (1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) ) is as follows:(10) $Rate = \frac{\text{d}\left[\text{A}\right]}{\text{d}t} = \frac{\text{d}\left[\text{B}\right]}{\text{d}t} = - k_{id,f} a_{TS} + k_{id, r }{a}_{TS}$(10)

Which, combined with eq. (3(3) $a = y*c$(3) ), (8(8) $K_f{}^{‡} = \frac{a_{TS}}{a_A a_B}$(8) ), and (9(9) $K_r{}^{‡} = \frac{a_{TS}}{a_C a_D}$(9) ) gives the following expression:(11) $Rate = \frac{d\left[A\right]}{dt} = \frac{d\left[B\right]}{dt} = -k_{id,f} K_f{}^{‡}\left[\text{A}\right] \left[\text{B}\right] \frac{y_A y_B}{y_{TS}} + k_{id,r} K_r{}^{‡}\left[\text{C}\right] \left[\text{D}\right]\frac{y_{C }{y}_D}{y_{TS}}$(11)

[A], [B], [C], and [D] are concentrations of reactants. In nonideal systems, $y$ values are ≠ 1 and, if they are not known, they are, by default, mathematically incorporated in the observed rate constant. In such situation, the following relationship holds:(12) $k_{reac,f} = k_{id,f} \frac{y_{A }{y}_B}{y_{TS}}$(12) where $k_{reac,f}$ is the observed rate constant for the forward reaction in eq. 1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) . Under ideal conditions, $y$= 1 and the observed rate constants are equal to the ideal ones. The same holds true for the reverse reaction.

The rate constant depends on the temperature of the system and its temperature dependence may be expressed in terms of the Eyring equation:(13) $k_{id,f }= \frac{k_{B}T}{h_p}\text{ exp} \left(-\frac{ \Delta G_{‡}^{\circ }}{RT}\right)$(13)

$k_B$ is Boltzmann’s constant, $h_p$ is Planck’s constant, $T$ is the absolute temperature, $R$ is the gas constant, and $\text{Δ}{G}_{‡}^{°}$ is the standard Gibbs energy change related to the formation of the TS from reactants.

The formalism expressed in the previous formulas only holds true in solutions that are ideally mixed and in general in all those cases where the rate of encounter of reactants is not limiting. When reactants diffuse slowly in the reacting medium, their frequency of collision can be the rate limiting step of the reaction and the observed rate can be very slow, despite the fact that the chemical reaction can be very fast (very high $k_{reac,f}$). In such situations, the observed rate constant is expressed by the following formula (only the forward reaction is shown):(14) $\frac{1}{k_{obs}} = \frac{1}{k_{reac}} + \frac{1}{k_{diff}}$(14) where $k_{reac}$ is the rate constant for the conversion of the reactants into products (eq. 6(6) $A + B \stackrel{K_f{}^{‡}}{\iff } TS \stackrel{k_{id, f}}{\to } C + D$(6) ) and $k_{diff}$ is the rate of the diffusion of the reactants. If the reaction occurs in conditions such that the diffusivity of reactants is not hampered, that is, $k_{reac} \ll k_{diff}$, $k_{obs}$ is entirely dominated by $k_{reac}$. The reaction is said to be reaction controlled. When molecular diffusion of reactants is hampered such that $k_{reac} \gg k_{diff}$ then, $k_{obs}$ is mostly determined by $k_{diff}$ and the reaction is said to be diffusion controlled.

If the reaction system is homogeneous, that is, all the physical and chemical properties are the same in any infinitesimal volume, then the rate of the reaction will be the same in any infinitesimal volume of the system and is also equal to the average rate. If the system is heterogeneous, that is, with properties different in different infinitesimal volume, then the reaction rate will vary locally and the average rate will be obtained by ideally integrating the reaction rate over the entire reaction volume.

Mechanistic modeling of the food matrix effect

When a chemical reaction occurs in food, the observed rate is, most likely, different from what would be measured under ideal conditions (a very dilute solution) and different from food to food. These observed rates are therefore “apparent” rates compared to the true rates measured under ideal conditions. The position of the reaction equilibrium can be affected also. This effect is buried within observed equilibrium and rate constants. In short, the food matrix will change thermodynamic and kinetic properties of the reaction by acting on concentrations of reactants and products (see eq. 2(2) $K_{th} = \frac{a_C a_D}{a_A a_B}$(2) and 11(11) $Rate = \frac{d\left[A\right]}{dt} = \frac{d\left[B\right]}{dt} = -k_{id,f} K_f{}^{‡}\left[\text{A}\right] \left[\text{B}\right] \frac{y_A y_B}{y_{TS}} + k_{id,r} K_r{}^{‡}\left[\text{C}\right] \left[\text{D}\right]\frac{y_{C }{y}_D}{y_{TS}}$(11) ), activity coefficients (eq. 4(4) $K_{obs} = K_{th} \frac{y_{C }{y}_D}{y_{A }{y}_B}$(4) and 12(12) $k_{reac,f} = k_{id,f} \frac{y_{A }{y}_B}{y_{TS}}$(12) ), and diffusivity (eq. 14(14) $\frac{1}{k_{obs}} = \frac{1}{k_{reac}} + \frac{1}{k_{diff}}$(14) ) of reactants and products, as well as on the temperature sensed by the reactants in each part of the system (eq. 5(5) $\frac{d\ln K_{th}}{d(1/T)} = - \frac{\Delta H_r^{\circ }}{R}$(5) and 13(13) $k_{id,f }= \frac{k_{B}T}{h_p}\text{ exp} \left(-\frac{ \Delta G_{‡}^{\circ }}{RT}\right)$(13) ). This is summarized and visually depicted in Figure 2.

Figure 2. The multiple effect of a food matrix on the rate and the equilibrium of a generic reaction. To be able to convert into products C and D, reactants A and B must come into contact to form an intermediate transition state TS. This, in turn, decomposes into products. The food matrix can affect the rate of this reaction in multiple ways: (1) Parallel reactions may deplete one of the reactants or both, thus changing concentrations of reactants. A special case is the effect of pH which changes the concentration of the reactive form of the reactant. (2) In the presence of catalysts the reaction may proceed through a different path which leads to a different transition state (TS’). (3) Partitioning of (one of) the reactants between different phase (e.g., oil and water phase) changes the concentrations of reactants in each phase. (4) Short-range and long-range interactions with ionic and nonionic solutes may alter reactants and transition state chemical potential (here the change in chemical potential is represented by reactant B being of different size which has obviously nothing to do with its molecular size). (5) The presence of obstacles (surrounding molecules, pores, etc.) or changes in viscosity may augment the effective path which reactant molecules must go through before they can collide, or may slow down reactant molecules. (6) Different thermal properties may produce a thermal profile within the food during nonequilibrium conditions, which changes the energy of reactants kinetic and, thus, frequency of collisions.

Mechanistic modeling the effect of food matrix on chemical reactivity means to be able to predict its effect on reactivity based on fundamental laws of chemistry and physics. In other words, it means to disentangle the effect of food matrix from the observed rate constant by quantifying its effect on concentrations of reactants (and intermediates), activities, diffusivities, as well as on temperature of the food systems. In the following subsections, theoretical frameworks for the quantification of each single effect are described.

Reactants concentration

Reactants can be consumed (or formed) through parallel reactions with (or from) other food components which would change their actual concentrations in the system. For instance, the condensation between a reducing sugar and a particular amino group will proceed at a slower rate when the sugar is simultaneously consumed in a parallel reaction with another amino group, or due to sugar isomerization. When parallel reactions involving reactants have been identified, their effect on the rate of the reaction of interest can be accounted for by reformulating the expression of the rate law. For instance, in the simpler case where reactant A in eq. 1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) is also involved in a parallel irreversible reaction:(15) $\text{A }+\text{ E }\to \text{ F }+\text{ G}$(15)

The rate law of eq. 11(11) $Rate = \frac{d\left[A\right]}{dt} = \frac{d\left[B\right]}{dt} = -k_{id,f} K_f{}^{‡}\left[\text{A}\right] \left[\text{B}\right] \frac{y_A y_B}{y_{TS}} + k_{id,r} K_r{}^{‡}\left[\text{C}\right] \left[\text{D}\right]\frac{y_{C }{y}_D}{y_{TS}}$(11) would be reformulated as follows:(16) $Rate = \frac{d\left[A\right]}{dt} = \frac{d\left[B\right]}{dt} = -k_{reac,f}\left[\text{A}\right] \left[\text{B}\right] + k_{reac,r} \left[\text{C}\right] \left[\text{D}\right] - k_{reac,p}\left[\text{A}\right] \left[\text{E}\right]$(16) where $k_{reac,p}$ is the rate constant for the parallel reaction. Of course, much more complex reaction networks may occur in food (for example also including parallel reactions where reactants are generated). In such situation the deployment of more sophisticated modeling approaches like multiresponse modeling are necessary (Van Boekel, 2008).

Reactants concentration can also change because of changes in the moisture content of food. If moisture loss is homogeneous throughout food, change in reactant concentration is inversely proportional to the change in moisture. However, in transient, nonequilibrium situations where the moisture content varies in time and throughout the food, mathematical models are needed to describe the diffusion of water and thus the actual moisture distribution throughout the food. This is important, for example, during high temperature processing of solid food, like cooking and frying. Changes in moisture content also result in changes in food structure and rheology and thus in the diffusivity of reactants (see next sections).

Catalysis

The presence of a catalyst changes the kinetics of a chemical reaction by providing an alternative route for the reaction to occur. The kinetics of the reaction in such a situation can be dealt with in the same way as described above (see eq. 16(16) $Rate = \frac{d\left[A\right]}{dt} = \frac{d\left[B\right]}{dt} = -k_{reac,f}\left[\text{A}\right] \left[\text{B}\right] + k_{reac,r} \left[\text{C}\right] \left[\text{D}\right] - k_{reac,p}\left[\text{A}\right] \left[\text{E}\right]$(16) ), that is, by reformulating the rate law as a function of two parallel reactions, one catalyzed and the other noncatalyzed, each characterized by its own reaction rate. Suppose a simple irreversible unimolecular conversion:(17) $\text{A }\to \text{ B}$(17)

The rate of this reaction is expressed as follows:(18) $Rate = \frac{d\left[A\right]}{dt} = - k_{reac}\left[\text{A}\right]$(18)

The alternative catalyzed route is:(19) $\text{A }+\text{ C }\to \text{ B }+\text{ C}$(19) where C is the catalyst. The rate of the catalyzed reaction is:(20) $Rate = \frac{d\left[A\right]}{dt} = - k_{reac,c}\left[\text{A}\right] \left[\text{C}\right]$(20)

The overall rate is therefore:(21) $Rate = \frac{d\left[A\right]}{dt} = - k_{reac}\left[\text{A}\right] - k_{reac,c}\left[\text{A}\right] \left[\text{C}\right] = \left(-k_{reac}- k_{reac,c}\left[\text{C}\right]\right) \left[\text{A}\right] = k_{obs } \left[A\right]$(21)

Thus, the rate constant of the noncatalyzed reaction and the rate constant of the catalyzed reaction both contribute to the observed rate constant. In food systems, acids and bases and metal ions represent the most common catalyst. Acids and bases can catalyze a reaction through the undissociated acid or the undissociated base or through generation of H⁺, OH⁻ ions, that is, through changes in pH. Each of the catalyst (undissociated acid, undissociated base, H⁺, OH⁻) provides an alternative route to the uncatalyzed reaction. Each route is characterized by its own reaction rate and each rate is proportional to the catalyst concentration. The overall reaction rate can be expressed as a function of multiple parallel reactions, as described in eq. 21(21) $Rate = \frac{d\left[A\right]}{dt} = - k_{reac}\left[\text{A}\right] - k_{reac,c}\left[\text{A}\right] \left[\text{C}\right] = \left(-k_{reac}- k_{reac,c}\left[\text{C}\right]\right) \left[\text{A}\right] = k_{obs } \left[A\right]$(21) . On top of that, pH can change the ratio between charged and neutral species of a reactant or product, that is, their pH-dependent species distribution (of course for ionisable reactants and products). This ratio can be calculated from the well-known Henderson-Hasselbalch equation. For instance, the relative distribution of the dissociated and undissociated forms of an acidic group is linked to pH by the following expression:(22) $pH = p{K}_a + {\log }_{10} \frac{a_{A^{-}}}{a_{HA}}$(22) where $K_a$ is the acid dissociation constant of the acidic group, and $a_{A^{-}}$ and $a_{HA}$ are the thermodynamic activities of the conjugated base and acid, respectively. By applying the above equation the actual activities of the reactive species ($a_{A^{-}}$ and $a_{HA}$) can be calculated at any pH and thus, the effect on reaction rate can be predicted.

Partitioning of reactants among different phases

In multiphase foods, reactants can partition between different phases, which changes their concentrations in each phase. A simple example is the distribution of a reactant between an aqueous phase and a micellar phase, the so-called micellar effect (van Boekel, 2009). When amphipathic compounds like mono- and diglycerides or phospholipids are present in solutions above the critical micelle concentration, a micellar phase can form which is in equilibrium with amphipathic molecules in solution. The distribution of reactants between the aqueous and the micellar phase can result in two distinct “reactors” wherein the same reaction occurs at two different rates because of differences in reactant concentrations. An example is the catalysis of the reaction between sorbic acid and thiols by surfactants (Wedzicha and Zeb, 1990) or bovine serum albumin (Wedzicha and Zeb, 1991). Another relevant example is the partitioning of reactants in emulsions. In such systems, reactants and intermediates can partition between the aqueous phase, the lipid phase or the interphase depending on their polarity and surface activity. The condensation between a reducing sugar and an amino group occurs with a different rate in oil in water nano-emulsions compared to aqueous solutions because of the pH-dependent partitioning of the amino acid between the lipid and the aqueous phases (Troise et al., 2016).

Partition of solutes in heterogeneous, multiphase systems can be quantitatively described by partition coefficients. A partition coefficient, $P$, is defined as:(23) $P = \frac{a_1}{a_2}$(23) where $a_1$and $a_2$ are the thermodynamic activities of the solute in phase 1 and phase 2, respectively, when the equilibrium between the two phases has been attained. Under ideal conditions, thermodynamic activities can be replaced by concentrations in eq. 23(23) $P = \frac{a_1}{a_2}$(23) . A partition coefficient depends on the solute, on the phases and also on temperature. The interfacial region between two phases represents an additional phase (even though it is not a true phase). This partition, together with the volume fraction of each phase in the system determines the molar concentration of each compound in each phase. The overall reaction rate is the sum of the contribution to the reaction from each phase. An example is provided for emulsions by (Romsted and Bravo-Diaz, 2013) but it can be formally extended to other systems made up of different phases or pseudophases (e.g., a micellar system). If we consider the bimolecular reaction in eq. 1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) , as long as the rate of partitioning between phases is much faster than the reaction itself (i.e., equilibrium has been attained for the reactants between the phases), the contributions to the overall reaction rate of the three separate reaction pathways in the continuous oil, interfacial, and aqueous regions are proportional to the concentrations of A and B in each region. This leads to the general equation:(24) $k_{obs} = k_2 \left[A_T\right] \left[B_T\right] = k_w \left[A_W\right] \left[B_W\right]{\text{ Φ}}_W + k_{O }\left[A_O\right] \left[B_O\right]{\text{ Φ}}_O + k_I \left[A_I\right] \left[B_I\right] {\text{Φ}}_I$(24) where $k_2$is the second-order kinetics rate constant; $k_w, k_O$, and $k_I$ are the rate constants for the reaction in each phase; [A] and [B] the reactant concentrations and φ are the volume fractions with φ_W + φ_O + φ_I = 1. Subscripts W, O, and I represent the water, oil, and interface phases, respectively. Eq. 24(24) $k_{obs} = k_2 \left[A_T\right] \left[B_T\right] = k_w \left[A_W\right] \left[B_W\right]{\text{ Φ}}_W + k_{O }\left[A_O\right] \left[B_O\right]{\text{ Φ}}_O + k_I \left[A_I\right] \left[B_I\right] {\text{Φ}}_I$(24) reduces to only one of the three terms displayed on the right-hand side if the reaction occurs entirely in one phase, that is, the reactants are soluble only in one phase. If the partition coefficients happen to be known, the amount of reactants in each phase can be calculated and eq. 24(24) $k_{obs} = k_2 \left[A_T\right] \left[B_T\right] = k_w \left[A_W\right] \left[B_W\right]{\text{ Φ}}_W + k_{O }\left[A_O\right] \left[B_O\right]{\text{ Φ}}_O + k_I \left[A_I\right] \left[B_I\right] {\text{Φ}}_I$(24) can be used to predict the overall rate. However, it must be stressed that rate constants may be different in different phases (i.e., $k_w, k_O$, and $k_I$in eq. 24(24) $k_{obs} = k_2 \left[A_T\right] \left[B_T\right] = k_w \left[A_W\right] \left[B_W\right]{\text{ Φ}}_W + k_{O }\left[A_O\right] \left[B_O\right]{\text{ Φ}}_O + k_I \left[A_I\right] \left[B_I\right] {\text{Φ}}_I$(24) may be different) because reactants are exposed to different milieus or because of differential properties of the solvent, for example, its dielectric constant. Partition coefficients must be determined experimentally, even though partition between oil and water can be also predicted by empirical and mechanistic models (Moln r et al., 2004; Molnár et al., 2004; Cheng et al., 2007).

Chemical potential

Since foods are nonideal systems, the chemical milieu around the reactants can change their chemical potential. This is a purely nonspecific effect exerted by the bystander molecules on those involved in the target reaction (in the sense that the bystander molecules do not chemically react directly with the molecules being affected). Neutral cosolutes will interact with the reactants through short-range forces, like dipole-dipole or dispersion forces. Ionic cosolutes will interact with neutral or charged reactants through long-range interactions like Coulomb forces. These interactions can be quite strong and may occur frequently in food which contains minerals but also dissociated acids and bases and also proteins and carbohydrates which carry charges (which depend on pH). The chemical potential of reactants can also be altered by a mere steric effect, that is, the physical occupancy of part of the available space by neighbor molecules. A typical example of this is the molecular crowding effect. This results in thermodynamic activity and diffusivity of reactants strongly deviating from those in ideal (aqueous) model systems (Ellis, 2001; Zhou et al., 2008). Since molecular crowding is ultimately a steric effect, this will be sensed more by high-molecular-weight compounds, like proteins and polysaccharides. However, the shape of the crowder also plays an essential role. Under crowding conditions, reaction equilibria are moved towards states that excludes the least volume to all the other macromolecules, for example, aggregation of macromolecules and folding of proteins are favored. Since protein unfolding is hindered under crowding conditions, heat stability of technologically relevant food enzymes is therefore expected to be higher in a crowded environment (Sasahara et al., 2003; Oliviero et al., 2014). As a rule of thumb, molecular crowding is expected to speed up transition-state-limited reactions which lead to a smaller free volume and to slow down diffusion-limited reactions because of the hindrance of reactants diffusion. The molecular crowding effect has been studied in the biochemistry domain for a while but its impact on kinetics and thermodynamics of food related chemical reactions has been very limitedly explored. Few notable exceptions are available in the scientific literature: the increase in egg white protein resistance to UV light in a crowded environment (Manzocco and Nicoli, 2012, 2015) and the increase in the rate of glycation of proteins in crowded systems (Perusko et al., 2015) (Weng et al., 2016).

In ideal solutions, that is, solutions where the average interactions between solute-solute, solute-solvent, solvent-solvent molecules are the same as in a pure liquid, the thermodynamic equilibria and the rate laws are expressed in term of concentrations of reactants. In nonideal solutions, the rate laws are expressed in terms of thermodynamic activity of the reactants and will vary accordingly. Once the activities of the probe molecules have been measured or estimated, the changes in equilibria and the reaction rates can be quantitatively accounted for (see eq. 4(4) $K_{obs} = K_{th} \frac{y_{C }{y}_D}{y_{A }{y}_B}$(4) and 12(12) $k_{reac,f} = k_{id,f} \frac{y_{A }{y}_B}{y_{TS}}$(12) ). Despite the fact that the relationship between activity and concentration is very simple (see eq. 3(3) $a = y*c$(3) ), the prediction of the activity coefficients of reacting species from fundamental laws of physics is quite difficult, especially in very complex systems like foods. Several mechanistic, empirical, and semiempirical models have been proposed to predict the value of the activity coefficients of probe molecules in nonideal systems.

Deviation from ideality may arise from short-range and long-range interactions between molecules. Long-range interactions refer to electrical potential experienced by charged molecules. Short-range interactions include, for example, dipole–dipole interactions and volume exclusion interactions. In electrolyte solutions nonideality arises from electrostatic interactions between solutes-solvents and solute-solute and those interactions are governed by Coulomb’s law, but also from short-range interactions. The first model to be proposed to estimate the activity coefficients of electrolytes was based on the Debye–Hückel theory (Hückel, 1924). In this theory the solution is considered as made up of charged ions of finite size surrounded by a cloud of opposite, point charges at a distance of minimum approach dispersed in a structureless medium of a determined relative permittivity. The model can only successfully predict the mean activity coefficient for dilute solutions when ion pairs are absent. To overcome this limitation, the mean spherical approximation (MSA) theory can be used (Blum, 1975). Although it has been developed for the description of thermodynamic properties of electrolyte solutions, those of nonionic cosolutes can be described as well. In the MSA approach, all the solute species, charged or uncharged, have finite size. The expression for the activity coefficients comprises a contribution from the electrostatic part and one for the volume exclusion effect via the hard sphere contribution, which reveals how space can be occupied by the cosolutes, which are approximated to have spherical shape. A very attractive alternative to calculate activity coefficients for electrolyte systems is represented by the electrolyte perturbed-chain statistical associating fluid theory (ePC-SAFT), which derives from the combination of SAFT-based equations with electrolytes theories. ePC-SAFT proved more accurate in predicting activity coefficients of biochemical compounds over Debye–Hückel theory (Hoffmann et al., 2013; Hoffmann et al., 2014). PC-SAFT accounts for molecular shape and size of solutes and their interactions with solvent by conceiving molecules as chains composed of spherical segments which interact through forces depending on their distance and diameter.

In solutions of neutral solutes, molecules experience only short-range interactions with neighboring particles. The UNIFAC model is one of the models that can be used to predict activity coefficients in such situations (Fredenslund et al., 1975). It is a modification of the UNIQUAC model. In both the models, the activity coefficients are split into two terms, a combinatorial term γ^C and a residual term γ^R. The first term quantifies the deviation from ideality due to differences in molecular shape. The latter is an enthalpy term that is based on the interactions between solutes. The difference between UNIFAC and UNIQUAC is that UNIFAC uses the functional groups present in interacting molecules to calculate the activity coefficients. In that respect, UNIFAC is a group contribution model, that is, a model that considers a whole solution as a mixture of functional groups rather than distinct molecules. Another example of group contribution model is the Savage-Wood Additivity of Groups (SWAG) model (Savage and Wood). This model is based on the assumption that each functional group X of solute i interacts with each group Y of cosolute C with a group-characteristic contribution that is independent of all other functional groups or relative positions and stereochemistry in solute i and cosolute C. The assumptions lead to an expression for the pairwise Gibbs energy interaction parameters in terms of specific group interactions. Other empirical models (Margules equation, van Laar equation, Wilson model, PC-SAFT, etc.) have been proposed for nonelectrolytes. A recent interesting development is represented by the conductor-like screening model for realistic solvation (COSMO-RS) model (Klamt and Schuurmann, 1993; Klamt, 1995; Klamt et al., 1998) and its subsequent reimplementations (e.g., COSMO-SAC) (Lin and Sandler, 2002). COSMO-RS model is not a group contribution model and does not require functional group parameters. This model uses quantum chemistry to calculate the chemical potential of molecular species in real solvents. First, COSMO equations are used to calculate the screening charge density on the surface of each molecule. From the pairwise interaction energies of surface patches, the chemical potential of a chemical species in a real solvent or mixture is predicted. The chemical potential of the solution or the mixture is approximately the weighted sum of the screen charge densities of its components. In analogy with the UNIFAC and UNIQUAC models, the COSMO-RS is a Gibbs excess free energy model and splits the chemical potential in a combinatorial and a residual contribution. However, compared to UNIFAC and UNIQUAC models (and other data interpolation models), the advantage of COSMO-RS models is that it does not need any adjustable, empirical group- or system-specific parameter, and chemical potentials are automatically predicted based on general or compound-specific properties, that is, COSMO-RS are more “mechanistic” compared to UNIFAC and UNIQUAC. In this respect, COSMO-MS has better predictive ability of thermodynamic properties of “new” compounds but is relatively more inaccurate in prediction of thermodynamic properties of compounds for which large amount of experimental data are available (this allows empirical parameters of (semi-)empirical models such as UNIFAC and UNIQUAC to be fitted to experimental data). Another important limitation of the COSMO-RS model is that it is, at the moment, inapplicable to electrolyte solutions (Klamt et al., 2010). The deviation from ideality due to purely volume exclusion effects can also be predicted, for instance via scaled-particle theory based on hard sphere approximation (Zhou et al., 2008). With the scaled-particle theory, within certain practical conditions, experimentally measured activity coefficients of proteins in crowded solutions can be estimated surprisingly accurately using simple structural models, in which rigid globular proteins are represented by hard particles having a size and shape similar to that of the protein at low resolution (such as spheres or spherocylinders), and large random coil polymers are represented by a random matrix of long rigid rods.

Outside chemical engineering, applications of thermodynamic models to predict chemical potentials of biological reactions are rare and can be found mostly in the pharmaceutical and the agrochemical areas (Klamt et al., 2002; Klamt et al., 2003; Klamt, 2005). A summary of the possible approaches can be found in (Held and Sadowski, 2016). Their application for the prediction of the activity coefficients of probe molecules in food is even more rare. Examples are the predictions of activity coefficient of aroma compounds in model juices or real juices, mainly based on the UNIFAC model (Carelli et al., 1991; Sancho et al., 1997) and those of amino acids in sugar solutions (Lakshmi and Nandi, 1976). More recently, the activity coefficient of calcium ions in milk-based systems enriched with disaccharides has been modeled using the MSA theory and used to describe the ion partitioning in milk and milk with added electrolytes (Gao et al., 2010a; Gao et al., 2010b). Other examples are the reaction of carmoisine, a food dye, with hydrogen sulphite ions (Wedzicha and Rumbelow, 1981) and the effect of ionic strength on dissociation of hydrogen sulphite, a food preservative (Wedzicha and Goddard, 1991). Thermodynamic activities can be measured also. To this end, equilibrium methods can be used, that is, methods based on partition of the compound between two phases such as osmotic measurements. This is because partition coefficients are based on activities rather than concentrations (see eq. 23(23) $P = \frac{a_1}{a_2}$(23) ). The activity of water, a_w in solutions, for instance, is measured from its vapor pressure above the solution which is based on the water partition between the liquid and the vapor phase. Measurement of an electric potential by means of an ion-selective electrode is also an option for charged molecules (pH measurement is one such example).

A specific case of thermodynamic activity that has a paramount importance in food is water activity, a_w. a_w can affect chemical reactivity in a direct way when water is a reactant (e.g., in a hydrolysis reaction) or in an indirect way, by changing concentration of reactants (but only if the change in a_w is associated with change in water content) or chemical potential of reactants. Concentration of solutes in food at reduced a_w also results in an increase in viscosity, volume occupancy or tortuosity and therefore, to a slower molecular diffusion.

Diffusivity

Differences in physical, structural and rheological properties of food affect the mobility of reactants within the food matrix which has a large effect on bimolecular (and higher order) reactions. Generally, molecular diffusivity can be described by two different approaches: Fick’s law of diffusion and the Maxwell–Stefan approach. Fick’s second law of diffusion can be written as follows:(25) ${\left(\frac{\partial c}{\partial t}\right)}_x = D_{ }\left(\frac{{\partial }^{2}c}{\partial x^2}\right)$(25)

This equation gives the change in concentration of the molecular probe in time at different points (x) as a function of the concentration gradient. Based on Fick’s law for diffusion, the diffusive behavior of compounds is described in terms of their diffusion coefficient $D_{ }$. The diffusion coefficient can be predicted from more fundamental laws in some simple cases. For particles or large molecules in a viscous fluid, the Stokes–Einstein equation can be applied:(26) $D = \frac{kT}{6\pi \mu R}$(26) Here, k is the Boltzmann constant, μ is the solvent viscosity, $T$ is the absolute temperature, and $R$ is the radius of the diffusing particle. In a solution that is not made up of polymers, it is the viscosity of the medium that determines the molecular diffusivity; however, it must be stressed that in several food systems, molecular diffusion cannot be described as diffusion of inert spherical particles in a viscous medium. In several food systems this looks like a very strong approximation, because reactants occur in crowded environments that limit molecular mobility by confining them within discrete molecular boundaries, like pores, and complex structural properties may be present. At molecular scale, the effect on molecule diffusion depends on the number, size, and shape of the surrounding molecules. When molecules are moving in a crowded environment, the steric hindrance exerted by the impenetrable macromolecules would create effective pores which constrain the free movement of the probe. This effect is known as molecular confinement.

A typical example is the diffusivity of a molecule in a plant cell or a dried fruit. In such a situation, molecular diffusivity is severely hindered, which reflects itself in diffusion coefficients smaller than the theoretical ones and which are more properly referred to as effective diffusion coefficients, D*. Effective diffusion coefficients are not just a function of solvent viscosity and density, but some structural features have to be considered, like the volume fraction that is available to the solute, the tortuosity (the ratio between the actual path travelled by a fluid element between two points divided by the straight line path between the same two points) and the pore size compared to the size of the diffusing molecule.

In the framework of the Maxwell–Stefan (MS) approach to diffusion the rate of diffusion depends on the balance between the driving forces to diffusion and the resistance exerted by the surrounding molecules and expressed in friction coefficients between each pair of diffusing molecules. This can be expressed as follows:(27) $- \frac{RT}{a_i} \frac{d{a}_i}{dz} = {\varphi }_{i,j} X_J \left(v_i - v_j\right)$(27)

The left-hand side of eq. 25(25) ${\left(\frac{\partial c}{\partial t}\right)}_x = D_{ }\left(\frac{{\partial }^{2}c}{\partial x^2}\right)$(25) expresses the “force” that causes the motion of component i compared to component j, whereas the right-hand side expresses the friction exerted by component j towards component i when the two compounds move with a relative velocity equal to $\left(v_i - v_j\right). X_J$ is the molar fraction of component j and ${\varphi }_{i,j}$is the friction coefficient between components i and j. Compared to the Fickian approach, the MS approach can account for all the forces responsible for molecular diffusion (electrical, gravitational) apart from chemical forces (concentration gradients) by adding relevant expressions of the driving forces in the left hand side of eq. 25(25) ${\left(\frac{\partial c}{\partial t}\right)}_x = D_{ }\left(\frac{{\partial }^{2}c}{\partial x^2}\right)$(25) . Moreover, with this approach the thermodynamic nonideality is incorporated in the driving forces for the diffusing probe molecule that is a function of its activity ($a_i$) rather than its concentration. If the activity coefficients can be calculated (see previous section) or measured, then the effect of food matrix on diffusion can be explicitly taken into account, whereas it is concealed in the diffusion coefficients of the Fickian approach. In other words, the MS diffusion coefficients are relatively less sensitive to changes in the food matrix (through the effect of differences in the chemical potential of the diffusing molecule) compared with the Fickian diffusion coefficient. The MS approach to diffusion has been applied in the food domain only to a limited extent. In one study, a MS multicomponent approach was developed to model salt and water diffusion in cheese, considering cheese as a three-component system consisting of NaCl (component 1), water (component 2) and a matrix of protein and fat (component 3) (Payne and Morison, 1999). In another study the diffusion of salt and water through sausage during drying was modeled (Costa-Corredor et al., 2010).

Temperature

Finally, the rate of the reaction depends on the temperature that is sensed by the reacting molecules but during food processing, heat transfer is very often in the unsteady state, in which temperatures are changing and materials are warming or cooling. In such situations a temperature profile is generated in time and within food which depends on the thermal properties of food. For this reason, food with different thermal properties that is exposed to the same external temperature will experience a different thermal gradient before thermal equilibrium is reached. The change in the surface temperature of food during heating/cooling is determined by its resistance to heat transfer from the external medium, for example, air, water, steam or oil, through natural or forced convection. This transfer can be calculated with the following general formula for heat transfer:(28) $Q = Ah \left(T_b - T_s\right)$(28) where $Q$ is the quantity of heat flowing, $h$ is the surface heat-transfer coefficient, $A$ is the surface area, $T_b$ is the bulk fluid temperature, and $T_s$ is the food surface temperature. The surface heat transfer coefficient is largely independent from the food matrix but mostly from the heating/cooling fluid.

The heat is transferred from surface to interior (or the other way around) by conduction. This can be calculated by using the Fourier equation, which is the equivalent of the Fick equation for heat:(29) ${\left(\frac{\partial T}{\partial t}\right)}_x = \alpha { }_{ }\left(\frac{{\partial }^{2}T}{\partial x^2}\right)$(29) where $\alpha$ is the food thermal diffusivity, which in turn is function of food thermal conductivity (λ), specific heat (C_p) and density (ρ):(30) $\alpha = \frac{\lambda }{\rho c_p}$(30)

Specific heat (C_p) is defined as the amount of heat required to increase the temperature of a unit mass of the substance by a unit degree. As it happens, C_p is an additive property, that is, C_p of a composite material is almost equal to the weighted sum of C_p for its components, that is, water, proteins, fats, carbohydrates, ice, and ashes (Choi, 1986). C_p of each component is, in turn, a function of temperature and can be calculated at different temperature from empirical equations (Choi, 1986). The thermal conductivity of a material is defined as a measure of its ability to conduct heat. The dependence of thermal conductivity of food components from temperature can be predicted from empirical equations (Choi, 1986). Unlike C_p, λ for a composite food cannot be predicted simply from its relative composition in the major components, but structural features like the geometrical distribution of the phases must be considered. Examples are the parallel models (the components are assumed to be placed parallel to heat flow), series models (components are perpendicular to the heat flow) and the Krischer model (mixed distribution). Other models such as Maxwell–Eucken and Kopelman models have been developed for two-component food systems consisting of a continuous and a dispersed phase (Sahin and Sumnu, 2006). However, heat can be transferred in food also through other mechanisms like convection, that is, mixing of fluid masses at different temperature or like distillation, that is, evaporation and condensation of vapour within a dispersed gaseous phase (e.g., through bubbles in bread). An additional and obvious complication in food is that, during thermal treatments, thermal properties may change because of the chemical and structural changes occurring.

Food matrix and enzymatic reactions

The effects displayed in Figure 2 still hold true when an enzymatic reaction is considered. The reaction rate will be affected by the concentration of substrates for the enzyme, their chemical potential and diffusivity, which, in the case of enzymatic reactions, also include the frequency of enzyme-substrate encounters. On top of that, the effect of the food matrix on the enzyme itself should be considered. Food matrix can affect enzyme catalytic activity via modification of enzyme conformation. This may happen via a number of chemical inhibitors or activators directly binding to the enzyme or indirectly through, for example, pH. In addition, the volume exclusion effect is expected to affect enzymatic reactions through modification of enzyme conformation, the quaternary structure of enzymes as well as the formation of the enzyme-substrate complex (Minton and Wilf, 1981). An example is provided once again by the GLs-MYR system described above. For instance, vitamin C is known to increase the enzymatic activity of MYR (Burow et al., 2006), and thus the conversion of GLs into breakdown products whereas pH and the concentration in Fe²⁺ affect the activity of the ESP which converts the intermediate product of MYR activity into epithionitriles (Hanschen et al., 2017). As a result the relative ratio between GLs enzymatic breakdown products strongly depend on food matrix properties.

Empirical modeling of the food matrix based on fingerprinting

The approach described so far, aims at predicting the effect of the food matrix on chemical reactivity based on the knowledge of the underlying mechanisms affecting reactivity in biological systems. A fundamentally different approach for modeling food matrix effects on chemical reactivity is represented by empirical black-box modeling based on machine learning and data mining techniques. In this approach, the parameter of interest is not predicted through mechanistic models, but through statistical correlations between the parameter of interest and a set of food matrix characteristics/properties. This hypothesis-free data-driven approach is briefly described in Figure 3. The first step is the quantification of the quality attribute of interest. This can be, for instance, the thermal stability of a target nutrient, measured by the rate constant for its degradation at a certain temperature. The second step is the extensive quantitative description of the chemical composition and the physical structure of several food matrices. We can name this comprehensive description of the food matrix its fingerprint. Chemical fingerprints can be obtained through a wide array of analytical techniques, such as liquid or gas chromatography (combined or not with mass spectrometry), spectroscopy, and electrophoresis. For example, chemical fingerprints can be constituted by the concentrations of a set of compounds, in the form of normalized intensities in a mass spectrum or areas of peaks in a chromatogram. It is important to stress that the identification of each single chemical descriptor is not formally necessary as long as each single variable is quantified with accuracy and precision. The physical structure of the system can also be described by combining outputs from microscopy, spectroscopy or imaging techniques. Spectroscopic data are inherently multivariate because the absorbance (or transmittance) of the sample is recorded at different wavelengths. When imaging techniques are used, these have to be translated into quantitative datasets, which is not straightforward. The authors believe that such translation of microscopy and other imaging techniques, along with the conversion of images into numerical data, would lead ultimately to empirical models that are able to predict complex quality attributes from the comprehensive numerical representation of the food matrix composition and structure. This scenario can be deemed realistic due to the fast improvement of calculation computer power as well as suitable data processing techniques. Hyperspectral imaging is a powerful technology for remotely inferring the material properties of foods. Hyperspectral images may be thought of as made up of a series of images of the same object taken at different wavelength or as a collection of pixels each represented by a spectrum. Taken as a whole, a hyperspectral image is simultaneously a chemical and structural fingerprint of the test material (Dale et al., 2013). Classical vibrational spectroscopy techniques, such as fluorescence spectroscopy and infrared spectroscopy also give information on the physical structure of test samples. Those techniques are typically used to describe the chemical composition of a test sample, but structural information can also be isolated with proper data processing. Other instrumental approaches to obtain quantitative information on the food matrix are, for instance, magnetic resonance imaging (MRI) and X-ray tomography combined with image analysis. Quantitative measures of physical/structural properties of the sample, such as porosity and tortuosity of a solid material, can be derived from these techniques (Schoeman et al., 2016).

Figure 3. Empirical modeling of the food matrix based on fingerprinting. This holistic hypothesis-free data-driven approach would consist in obtaining the values of a thermodynamic or kinetic parameters for several food matrices (1(1) $\text{A }+\text{ B }\rightleftarrows \text{C }+\text{ D}$(1) ) along with chemical and physical/structural fingerprint of each of them (2(2) $K_{th} = \frac{a_C a_D}{a_A a_B}$(2) ) and translate this fingerprint into a multivariate quantitative dataset (MQD, 3). Chemical fingerprinting can be obtained through many analytical techniques like, chromatography. The physical structure of the system can also be quantitatively described by combining outputs from microscopy, spectroscopy and imaging techniques combined, for example, with image analysis. Once MQDs are obtained for several matrices multivariate statistical techniques are applied to calibrate the MQD against the measured thermodynamic or kinetic property (4(4) $K_{obs} = K_{th} \frac{y_{C }{y}_D}{y_{A }{y}_B}$(4) ). This empirical model allows (1) to predict the thermodynamic of kinetic behavior of future samples once their fingerprinting is obtained (2) to identify a set of chemical and physical markers that can then be used to set off new hypothesis-driven experiments in order to develop proper mechanistic models.

Once the quality attribute of interest and the fingerprints have been obtained for several food matrices, a calibration model is developed that correlates the property of interest, that is, the dependent variable, to the multivariate quantitative description of the food matrix, that is, the independent variables. The derivation of the calibration model from a multivariate independent variable requires the use of special modeling tools. Examples of such modeling techniques include: basic multivariate regression algorithms, such as principal component regression (PCR) and partial least square regression (PLSR); and also statistical learning algorithms inspired by biological neural networks such as artificial neural network (ANN); support vector machine (SVM); back propagation (BP) methods; self-organizing maps (SOM); or genetic algorithms (GA). These pattern recognition techniques are also said to be supervised because the models are trained against a desired output, that is the prediction of the property of interest in a new sample. The learning algorithms inspired by biological neural networks have the advantage over traditional MLR, PLSR, and PCR to be able to deal better with highly nonlinear systems.

Once the calibration model has been developed, this can be used to predict the property of interest in a new food, once its fingerprint is obtained. The advantage of using this modeling approach to analyze a complex problem is accompanied by the drawback that it only allows the prediction of the behavior of a certain system, but does not provide any knowledge on the food system. Application of machine learning combined with data mining is not a completely new approach in food science. This type of approach has been, for instance, used for authenticity assessment and fraud detection (Capuano and van Ruth, 2012; Cozzolino, 2014; Nunes, 2014), food spoilage detection (Ellis et al., 2004; Nicolaou and Goodacre, 2008). More recently, an interesting use of such approach has been proposed to discriminate between different processing methods, storage conditions or food quality classes (Vervoort et al., 2012; Kebede et al., 2013; Vervoort et al., 2013; Grauwet et al., 2014). Here we propose that this approach can also be used to predict the effect of the food matrix on chemical reactivity, for example, which is the expected rate constant or equilibrium constant for a certain reaction in a certain food. In addition to regression (where a quantitative output is predicted), classification can be also applied to allocate food matrices in broad qualitative categories, for instance more or less “protective” towards a certain nutrient. Along with prediction or classification, this approach can be used to select discriminant features that can be, in turn, tested in hypothesis-driven mechanistic models. This can be done by selecting single or patterns of independent variables that are most correlated to the dependent variable. The only example of such a use of this empirical approach reported in the scientific literature so far can be found in (Hennig et al., 2013) where an untargeted metabolomics approach, followed by random forest regression, was applied to identify metabolites associated to thermal GL degradation in a segregating Brassica oleracea population.

Conclusions and future perspectives

Foods are very complex systems from a chemical and structural point of view. There are, therefore, several food-related factors that can affect the evolution of a chemical reaction. We have named the totality of these effects the food matrix effect. The food matrix effect is responsible for the biological variation we observe in the thermodynamics and kinetics of the same reaction occurring in different foods. The striking differences that may be observed between reactions occurring in a food, compared to an ideal aqueous system, should also alert food scientists about the proper choice of model systems when the findings must be extended to real food with an acceptable level of confidence. In this report we have presented two different approaches to describe the effect of the food matrix on reactivity. The first is a mechanistic approach based on fundamental laws of chemistry and physics. This approach should permit to predict the behavior of a food system based on the detailed knowledge of its composition and structure. Two considerations are evident with respect to this approach. The first is that the definition of mechanistic is somewhat arbitrary, since at a certain point empirical or semiempirical models must be employed. For instance, the effect of a cosolute on rate constants can be incorporated into kinetics via a fundamental law of chemistry, that is the relationship between concentration and chemical activity, but the prediction of an activity coefficient is sometimes based on (semi-)empirical models, though much progress has been made in this respect. The second consideration is that, if it is true that each single effect can be quantitatively described based on mechanistic models, the simultaneous occurrence of several effects makes the theoretical prediction of the behavior of a real food system still very far to reach. For instance, the change in one ingredient in a food formulation can simultaneously change the food pH, the hydration properties, the chemical potential, the reaction pathway, and also the diffusion behavior of reactants and intermediates. And those changes may occur in different directions so that the behavior of the system may be perceived as casual or nonlinear. In this respect, a holistic or “system” approach, which would take into account the complexity of the food system and the potentially nonlinear relationship between certain chemical or physical properties and the quality attribute of interest, would be useful. System dynamics modeling is just one of those “system” approaches. The second approach to predict the effect of food matrix is a purely empirical, correlative approach, based on an extensive characterization of the chemical and structural properties of the food matrix and a correlation of the obtained matrix fingerprint to chemical reactivity. It is clear that this approach would require reliable analytical techniques capable of producing large amount of precise and accurate data, as well as intensive computational power, which is not expected to be a limiting factor in the near future. This approach may be used for accurate predictions of the effect of food matrix on reactivity and, by the identification of key characteristics of interest of the food matrix, can be the starting point for further mechanistic investigations. In such a way, this inductive, top-down, hypothesis-free approach will proceed hand in hand with the deductive, hypothesis-based, bottom-up approach in generating knowledge about the effect of food matrix on chemical reactivity.

Acknowledgment

The authors declare no conflict of interest.

Appendix.

Overview of kinetic parameters reported in the scientific literature for the thermal degradation of glucosinolates in different food matrices

Table

Compound	Matrix	Reaction order	Reaction rate constant	Activation Energy	Reference
	Brassica vegetables	First order	k100˚C × 10^–2 (min^–1)	Not estimated	(Dekker et al., 2009)
gluconapin	red cabbage		0.2
	broccoli		0.1
	Brussels Sprouts		2.1
	pak choi		0.3
	Chinese cabbage		1.9
glucobrassicin	red cabbage		0.8
	broccoli		1.5
	Brussels sprouts		4.9
	pak choi		2.5
	Chinese cabbage		2.7
4-methoxyglucobrassicin	red cabbage		1.7
	broccoli		5.0
	Brussels sprouts		6.8
	pak choi		3.3
	Chinese cabbage		2.9
	Broccoli with different water contents (%)	First order	k100˚C ± SD (× 10^–2, min^–1)	kJ mol^−1 ± SD	(Oliviero et al., 2012)
glucoraphanin	13		2.7± 0.3	159 ± 7
	34		6.7 ± 0.4	96 ± 5
	56		4.1 ± 0.3	99 ± 3
	82		1.5 ± 0.1	113 ± 2
glucobrassicin	13		1.7 ± 0.1	170 ± 7
	34		3.4 ± 0.4	82 ± 5
	56		2.6 ± 0.3	86 ± 2
	82		1.3 ± 0.1	110 ± 3
	Brassica vegetables	Second order	ksteaming temp ± SD (min^–1)	Not estimated	(Giambanelli et al., 2016)
Total indolic glucosinolate	broccoli lavagnino		0.026 ± 0.004
	nero di Toscana		0.039 ± 0.008
Total aliphatic glucosinolate	broccoli lavagnino		0.093 ± 0.041
	nero di Toscana		0.358 ± 0.137
	Brassica vegetables		kfrying temp ± SD (min^–1)
Total indolic glucosinolate	broccoli lavagnino		0.146 ± 0.017
	nero di Toscana		0.138 ± 0.019
Total aliphatic glucosinolate	broccoli lavagnino		0.183 ± 0.044
	nero di Toscana		1.076 ± 0.221	Not estimated	(Hennig et al., 2014)
	Brassica oleracea, doubled haploid population	First order	k100˚C (min^−1)
glucoraphanin	AG 6105 (2009)		0.011
	AG 1017 (2009)		0.019
	AG 6026 (2009)		0.030
glucobrassicin	AG 6105 (2009)		0.019
	AG 1017 (2009)		0.038
	AG 6026 (2009)		0.045
	Brassica vegetables from different seasons and year	First order	k100˚C ± SD (× 10^–2, min^–1)	Not estimated	(Hennig et al., 2012)
glucoraphanin	broccoli (autumn 2010)		1.4 ± 0.1
	Chinese kale (autumn 2010)		1.5 ± 0.1
	broccoli (spring 2011)		2.6 ± 0.0^a
	Chinese kale (spring 2011)		2.4 ± 0.1
glucobrassicin	broccoli (autumn 2010)		3.0 ± 0.1
	Chinese kale (autumn 2010)		3.1 ± 0.1
	broccoli (spring 2011)		3.7 ± 0.1
	Chinese kale (spring 2011)		3.8 ± 0.1
	Broccoli binary systems	Second order	k100˚C ± SD (min^−1)	Not estimated	(Giambanelli et al., 2015)
glucoraphanin	control (only broccoli)		0.038 ± 0.003
	broccoli/onion ratios (1/9)		0.015 ± 0.002
	broccoli/onion ratios (3/7)		0.020 ± 0.002
	broccoli/onion ratios (1/1)		0.029 ± 0.003

^a 0.0 = 0.036.