# Modelling the Maillard reaction during the cooking of a model cheese

Emmanuel Bertrand, Xuan Mi Meyer, Elizabeth Machado-Maturana,
Jean-Louis Berdagué, Alain Kondjoyan

► To cite this version:

Emmanuel Bertrand, Xuan Mi Meyer, Elizabeth Machado-Maturana, Jean-Louis Berdagué, Alain Kondjoyan. Modelling the Maillard reaction during the cooking of a model cheese. Food Chemistry, 2015, 184, pp.229-237. 10.1016/j.foodchem.2015.03.097. hal-01852722

HAL Id: hal-01852722

https://amu.hal.science/hal-01852722

Submitted on 6 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Open Archive Toulouse Archive Ouverte (OATAO)

OATAO is an open access repository that collects the work of some Toulouse researchers and makes it freely available over the web where possible.

This is an author's version published in: http://oatao.univ-toulouse.fr/20582

Official URL: https://doi.org/10.1016/j.foodchem.2015.03.097

To cite this version:

Bertrand, Emmanuel and Meyer, Xuân-Mi and Machado-Maturana, Elizabeth and Berdagué, Jean-Louis and Kondjoyan, Alain Modelling the Maillard reaction during the cooking of a model cheese. (2015) Food Chemistry, 184. 229-237. ISSN 0308-8146
# Modelling the Maillard reaction during the cooking of a model cheese

Emmanuel Bertrand<sup>a,*,1</sup>, Xuân-Mi Meyer<sup>b,c</sup>, Elizabeth Machado-Maturana<sup>a</sup>, Jean-Louis Berdagué<sup>a</sup>, Alain Kondjoyan<sup>a</sup>

<sup>a</sup>Institut National de la Recherche Agronomique (INRA), UR 370 QuaPA, 63122 Saint-Genès-Champanelle, France

<sup>b</sup>Université de Toulouse, INPT, UPS, Laboratoire de Génie Chimique, 4 Allée Emile Monso, F-31030 Toulouse, France

<sup>c</sup>CNRS, Laboratoire de Génie Chimique, F-31030 Toulouse, France

## ABSTRACT

During processing and storage of industrial processed cheese, odorous compounds are formed. Some of them are potentially unwanted for the flavour of the product. To reduce the appearance of these compounds, a methodological approach was employed. It consists of: (i) the identification of the key compounds or precursors responsible for the off-flavour observed, (ii) the monitoring of these markers during the heat treatments applied to the cheese medium, (iii) the establishment of an observable reaction scheme adapted from a literature survey to the compounds identified in the heated cheese medium (iv) the multi-responses stoichiokinetic modelling of these reaction markers. Systematic two-dimensional gas chromatography time-of-flight mass spectrometry was used for the semi-quantitation of trace compounds. Precursors were quantitated by high-performance liquid chromatography. The experimental data obtained were fitted to the model with 14 elementary linked reactions forming a multi-response observable reaction scheme.

*Chemical compounds studied in this article:*

2-Methylbutanal (PUBChem CID: 7284)

Formic acid (PUBChem CID: 284)

Furfural (PUBChem CID: 7362)

5-Methylfurfural (PUBChem CID: 12097)

Furfuryl alcohol (PUBChem CID: 7361)

Isomaltol (PUBChem CID: 18898)

Lactulosyllysine (PUBChem CID: 3082392)

Maltol (PUBChem CID: 8369)

Melanoïdins

Pyrazine (PUBChem CID: 9261)

**Keywords:**

Processed-cheese

Maillard reaction

Volatile compounds

Multi-responses modelling

## 1. Introduction

Processed cheese derives from the transformation of dairy (Gouda, Cheddar, Emmental, butter and various milk powder) and non dairy (emulsifiers, texturizing agents, aroma) ingredients into one homogenous generally spreadable product with a long shelf life (Caric, 2000; Kapoor & Metzger, 2008).

The mechanical and thermal settings necessary to get a microbiologically safe product with both colour and texture desirable for consumers are well known and could readily be modelled. For instance, the thermal settings necessary to get a microbiologically safe processed cheese could readily be calculated from the

parameters of Bigelow and Weibull (van Boekel, 2002). Colour defects promoted by the application of inadequate thermal settings in relationship with the composition of the cheese medium and linked to the Maillard reaction have been extensively studied (Bley, Johnson, & Olson, 1985a, 1985b). The rearrangement of caseins by emulsifying salts giving rise to a creamed texture has been the topic of much research (Lee, Buwalda, Euston, Foegeding, & McKenna, 2003; Panouille, Durand, Nicolai, Larquet, & Boisset, 2005).

In fact consumers request a microbiologically safe product with optimal colour, texture and taste. However it is very unlikely that these four responses reach their optimal properties for the same formulations and processing parameters. Therefore the best compromise has to be found. Methods for an accurate quantitation of colour, texture and microbiological safety exist. However this is not the case for taste and flavour and their accurate quantitation remain a major analytical challenge. For this reason, there is still

* Corresponding author. Tel.: +33 (0) 473 40 55 21; fax: +33 (0) 473 40 78 29.
*E-mail address: emmanuel.bertrand@univ-bpclermont.fr (E. Bertrand).

¹ Permanent address: Institut Pascal, UMR CNRS 6602, Axe GePEB - BP 10448, F-63000 Clermont-Ferrand, France.
a need to integrate the development of flavour into the multimodal strategies for the optimisation of processed cheese quality.

Control and optimisation of flavour properties has been recently described as "the ultimate challenge for the food and flavour industry" (Parker, 2013). Up to now the multi-response stoichiokinetic models of the Maillard reaction have been mostly applied to model systems (Brands & van Boekel, 2001) but rarely to real food products (Parker et al., 2012). In contrast to model systems, that are usually buffered and exposed at a single temperature, the pH of the food matrix is susceptible to decrease as weak acids such as formic or acetic acids are formed. In addition, the temperature of the cheese medium varies in order to stick to the thermal parameters that are usually applied during the elaboration of industrial processed cheese. Moreover most of the multi-response stoichiokinetic model studies do not focus on the formation of the odorous volatile compounds. (Brands & van Boekel, 2002, Martins & Van Boekel, 2005). Indeed, the formation of these volatiles and odorous compounds is a crucial step in the Maillard reaction and many elementary reactions are combined, as has been shown by studies using stable isotopes labelling (Yaylayan, 1997, 2003). Finally, the quantitation of volatile compounds in complex multiphasic food matrices, whose desorption properties potentially vary during their mechanical and thermal elaboration remains until now, a major analytical challenge (Samavati, 2013).

In a previous study (Bertrand et al., 2011), we determined from a qualitative point of view that lipid oxidation, caramelisation and Maillard reaction are responsible for most of the changes occurring in the volatile fraction during the thermal treatments applied to processed cheeses. Some of the molecules originating from these reactions are already known as responsible for "off flavour" defects. In particular, we identified two molecules, maltol and furaneol, produced during the Maillard reaction as the main contributors to "overcooked" defects.

The aim of this work is to move from a qualitative approach toward a quantitative model that could be integrated into multi-criteria optimisation strategies for the prediction of processed cheese quality. Therefore, this study is focused on the Maillard reaction, as it was found to be the main source for off-flavour identified during cooking (Bertrand et al., 2011). We were led (i) to extract an observable reaction scheme from the data contained in the volatile fraction of the processed cheeses and the literature available, (ii) and to model the evolution of the key compounds using a multi-responses stoichiokinetic model. Such a model consists of an intricate network of reactions. A detailed guideline concerning the establishment and the resolution of a multi-response stoichiokinetic model can be found for instance in the book of Van Boekel (2009a, chap. 8, 2009b, chap. 14).

The present work must not be seen as a fundamental work (generally conducted in a simplified binary mixture with perfectly controlled pH and temperature) aiming at the removal of theoretical and analytical locks for a better understanding of the Maillard reaction, but as an attempt to get the best of the current knowledge available in order to improve the flavour quality of industrially processed cheeses. This is consequently the first step of a sequential design aiming at a better understanding of the processed cheese system.

# 2. Material and methods

## 2.1. Composition, formulation and cooking of the cheese

### 2.1.1. Composition and formulation

Micellar casein native and milk permeate were purchased from Ingredia (St-Pol-sur-Ternoise, France). Anhydrous milk fat was from Campina (Amersfoort, The Netherlands). Sodium chloride

and citric acid were purchased from Sigma-Aldrich (St. Louis, MO) and were of analytical grade. Deionised water and a mixture of sodium polyphosphates (Pitkowski, Nicolai, & Durand, 2008) were also used. The final composition per 100 g of cheese was approximately 60 g of water, 20 g of fat, 12 g of protein, 6 g of lactose and some other minor constituents. More details about the manufacturing process are provided in our previous study (Bertrand et al., 2011).

### 2.1.2. Cooking system

A cooking system was designed for heating a cheese sample of about 10 g to a final temperature of 80-150 °C as quickly as possible and to maintain this temperature for a given time. It is possible to reach 150 °C in about 3 min and 30 s. The specifications of the system, its operational parameters and performances are also described in Bertrand et al. (2011). The temperature was measured by using a type-K thermocouple placed at the core of the cheese medium in a set of preliminary experiments. It was removed during the experiments in order to prevent any contamination of the volatile fraction by the probe.

Cooking conditions used for this study are shown in Fig. 1. Each point corresponds to a triplicate run including formulation, cooking and analysis steps. As the device is conceived to withstand the pressure for cooking at temperature above 100 °C, it is not possible to regularly take a sample during a single backing and only one sample, corresponding to the final stage of a single experiment could be taken out. Because of this, the number of samples taken is necessarily reduced. In order to get the most information possible, the times intervals were reduced when working at higher temperatures.

## 2.2. Assay of non-volatile compounds

The method used to quantify the sugar content (lactose and galactose) was adapted from Rocklin and Pohl (1983). Two grams of processed cheese were dissolved in 100 mL of deionised water; 2 mL of acetic acid (10%v/v) were added and the pH of the solution was adjusted to 4.6 with about 2 mL of 1 M sodium acetate. The sugar content was quantified by ion exchange HPLC and detected by amperometry. The quantitation limits for the two products were in the range of 10 mg for 100 g of processed cheese.

The free amino acid contents were measured chromatographically according to French standards (AFNOR XP V 18-113, January 1998 and AFNOR XPV 18-114, January 1998 for tryptophan).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Overview of the experimental design. Numbers represent the time spent at the selected temperature (warm up period excluded). The thermal treatments were carried out according to the methodology described in Bertrand et al. (2011). Processed cheese samples were analysed by solid-phase microextraction associated with comprehensive gas chromatography-time of flight mass spectrometry.</figcaption></figure>
Depending on the free amino-acid the quantitation limit was between 1 and 10 mg for 100 g of processed cheese.

Furosine is a non-volatile compound that is an indirect marker of lactulosyllysine, an Amadori compound formed from lactose and lysine. The furosine content was quantified by ion-pair reversed-phase HPLC. The separation was isocratically conducted at 25 °C using a Phenomenex (Torrance, CA) Jupiter column (250 × 4.6 mm, particle diameter 5 µm) with an 80/20/0.1% (v/v/v) water, acetonitrile, TFA eluent and detection at 280 nm using the reference method (ISO:DIS 18329, 2001). The Amadori content was then calculated with a 3.1 conversion factor according to Brands and van Boekel (2001). The quantitation limit was in the range of 0.1 mg of furosine in 100 g of cheese.

The colour was measured in the CIE Lab colour space on the upper side of the cheese medium using a Konica-Minolta spectrophotometer CM-2500d. Five measurements were performed in SCE mode with a viewing angle of 10° and a D65 illuminant and the average measurement was kept. The colour index (ΔC) was calculated using the following formula:

$$ \Delta C = \sqrt{(L_0 - L)^2 + (a_0 - a)^2 + (b_0 - b)^2}, $$
with L, the luminance (0 (black) < L < 100 (white)), a the red component (−120 (green) < a < 120 (red)) and b the yellow component (−120 (yellow) < b < 120 (blue)), while L₀, a₀ and b₀ are the values for the raw cheese medium.

### 2.3. Solid-phase microextraction coupled to comprehensive gas chromatography/time-of-flight-mass-spectrometry

A 75 µm solid-phase microextraction (SPME) Carboxen/PDMS fibre (Supelco, Bellefonte, PA) was used for the static headspace extraction of volatiles from the processed cheese sample. The fibre was conditioned before analysis according to the manufacturer's recommendations (1 h at 175 °C). For each sample, 1.5 g of processed cheese were dropped in a sealed 20 mL vial and pre-incubated for 30 min at 60 °C. The SPME needle was then exposed to the cheese headspace for 30 min at 60 °C and the volatiles were thermally desorbed from the fibre in the liner of the GC injection port for 1 min at 280 °C. Splitless injection was carried out with a Combipal autosampler (CTC Analytics AG, Zwingen, Switzerland) on a 6890 N chromatograph (Agilent, Massy, France) integrated in a GC × GC/TOF-MS LECO Pegasus instrument (LECO Corporation, St. Joseph, MI). The first chromatographic separation was done using an SPB-5 capillary column (5% diphenyl-, 95% dimethylpolysiloxane, length 30 m × internal diameter 0.32 mm × film thickness 1 µm; Supelco, St-Germain-en-Laye, France and the second separation using a DB-17 capillary column (50% dimethyl-, 50% diphenylpolysiloxane, length 2.50 m × internal diameter 0.178 mm × film thickness 0.30 µm; J&W Agilent, Santa Clara, CA). To achieve reliable identification, the data obtained were cross-matched against spectral databases (NIST/EPA/NIH 2005 v2.0d, Gaithersburg, MD; Wiley 275K, 1996; Masslib 1999; MSP Köfel, Koeniz, Switzerland) and comparisons with databases of linear retention indices (Kondjoyan & Berdagué, 1996). Data were expressed in arbitrary units of area.

### 2.4. Multi-response modelling

The multi-response stoichiokinetic model was developed using Matlab R2009a® software with the associated optimisation and statistics toolboxes. The system of ordinary differential equation (ODE) obtained was solved using time discretisation by the ode15s solver of Matlab. As the system was placed under transient conditions (during the temperature increase), the resolution of the system required the use of very different time scales. Therefore a solver for stiff systems such as the ode15s was found to be more

suitable than the ode45 (that turned out to be unsuccessful) (Shampine & Reichelt, 1997). The ode15s solver is based on the implementation of the numerical differentiation formulas of order 1–5. The fit of experimental data by the model was made according to the method of maximum likelihood. Given the number of parameters to be settled, we choose to conduct this recognition procedure under positivity constraint. In fact, the values of the rate constants and activation energy are necessarily positive. Therefore we used the fmincon function from Matlab which allows to bound the vector of parameters to adjust. However, this function also restricts the number of algorithms available for the minimisation procedure. We used the active-set algorithm (Rakowsha et al., 1991) in order to minimise the opposite of the likelihood function. After the fit of the model parameters to the experimental data, the confidence intervals were calculated from the Fisher information matrix. Briefly, the covariance matrix of the estimated parameters is obtained from the inverse of the Fisher information matrix. It leads to the determination of the standard error made for each of the parameters and allows the determination of the corresponding confidence intervals using the table of Student critical values.

The fitting procedure was made according to two steps. In the first one, rate constants were adjusted at a single temperature (80, 100, 120, 135 and 150 °C, respectively). If the rate constant follows Arrhenius's law, the linear regression of ln(k) as a function of 1/T will give rise to a straight line with slope −Eₐ/R and intercept ln(k₀), where k₀ is the pre-exponential factor [s⁻¹] and Eₐ the activation energy [kJ mol⁻¹]. Given 4 degrees of freedom, this will be the case, if the absolute value of the Pearson's correlation coefficient (r) associated with the regression line is greater than 0.88 (p < 0.05). In the second step, the values obtained from these linear regressions were taken as the initial values for the parameters and their uncertainties were used to set the associated boundaries of the vector of parameters to be adjusted by the fitting algorithm.

## 3. Results and discussion

### 3.1. Extraction of an apparent reaction scheme of the Maillard reaction occurring in processed cheese

Among 346 volatiles compounds identified by GC × GC/TOF-MS, 81 are significantly affected by the thermal treatment applied to the processed cheese medium (Bertrand et al., 2011). These structures came primarily from the degradation of lipids, involving hydrolysis and oxidation reactions, from the degradation of amino acids, particularly involving the Maillard reaction, and from the degradation of carbohydrates, including caramelization reactions, without any intervention of an amino acid. This study also suggested that the products of the Maillard reaction affect the smell of the cooked cheese medium more than those produced by the oxidation of fatty acids. Furaneol and maltol can be considered as the main cause of “overcooked” defects encountered in processed cheese. In this context, we chose to focus our modelling efforts on the Maillard reaction.

About 40 of these compounds could be associated either with the Maillard reaction or caramelisation. On the one hand, some compounds such as furfural or 2-methylpyrazine present very important variations associated with the thermal treatments being applied to the processed cheese medium. However, they have not been identified as flavour carriers through olfactometric analyses. This category of constituents would be good reaction markers. They will provide information on the extent of the reaction and could be easily included in a reaction scheme. However, they will be of little interest for the direct determination of the flavour of the product. On the other hand, some compounds such as furaneol could only be detected and quantitated in a very small number of
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Observable reaction scheme of the Maillard reaction applied to cheese media. Numbers (R01, R02, up to R14) represent the reaction that is considered. Green colour stands for the components that are quantitated in the cheese matrix. Orange colour stands for the one that are semi-quantitated by the analysis of the volatile constituents or by colour measurements. The compounds that are not quantitated are represented with the red colour. The different components are abbreviated as follow: Y1 lac: lactose; Y2 deo: deoxyribose; Y3 af: formic acid; Y4 fm: furfuryl alcohol; Y5 ama: lactulosyllysine (measured by its acid hydrolysis product: furosine); Y6: gal: galactose; Y7: dg1: 1-deoxyglucosone; Y8: mal: maltol; Y9: isom: isomaltol; Y10: dg3: 3-deoxyglucosone; Y11: hmf: 5-(hydroxymethyl)furfural; Y12: mf: 5-methylfurfural; Y13: fur: furfural; Y14: glyo: glyoxal; Y15: 23bd + 3h2b: 2,3-butanedione in redox equilibrium with 3-hydroxy-2-butanone; Y16: mb: 2-methylbutanal; Y17:aminoc: aminoketone; Y18: pyr: pyrazine; Y19: mela: melanoidins; lys: lysine; isol: isoleucine; aceta: acetaldehyde. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \frac{d[\text{aminoac}]}{dt} = k_{12} \times [\text{glyo}] - 2 \times k_{13} \times [\text{aminoac}] \quad (\text{R12, R13}) $$

$$ \frac{d[\text{pyr}]}{dt} = k_{13} \times [\text{aminoac}] \quad (\text{R13}) $$

$$ \frac{d[\text{mela}]}{dt} = k_{14} \times ([\text{mal}] + [\text{isom}] + [\text{mf}] + [\text{fur}] + [23bd] + [\text{mb}]) \quad (\text{R14}) $$

## 3.2. Challenges associated with the quantitation

In order to model the proposed reaction scheme, the individual components must be expressed in moles. However, the volatile compounds could only be experimentally semi-quantitated (for a given compound, data are expressed in arbitrary units of area of the chromatographic peak). This is explained by the analytical challenges posed for the quantitation of melanoidins and for the quantitation of small and highly reactive molecules in a multiphasic medium. Moreover, the texture properties of this medium are considerably varying during the treatments that are applied. For instance, Saint-Eve, Juteau, Atlan, Martin, and Souchon (2006) demonstrated the influence of gel structure associated with changes in the microstructure network of caseinate on the flavour release of aroma compounds in flavoured stirred yoghurt. They found that the release decreased for most of the aroma compounds when the yoghurt exhibited a higher viscosity. In their study, the increase of the gas/matrix partition coefficient ranged from 10% to 300% (for ethyl butanoate and diacetyl, respectively). Another study from Deleris, Atlan, Souchon, Marin, and Trelea (2008), underlines the effect of fat content on the apparent diffusivity of hydrophobic contents in

yoghurt. They found a 15-fold decrease for the apparent diffusion coefficient of linalool and even 50-fold for ethyl-hexanoate. Furthermore, it seems reasonable to suspect an effect of protein-lipid interactions on the release of aroma compounds.

Therefore, for each volatile compound (i) an apparent global quantitation coefficient ($K_{p/m}^l$) was added to convert the semi quantitated data (chromatographic peak area (p) of the compound (i)) to quantitated data (content of compound (i) in the matrix (m)). This led to add one additional parameter to the model for each volatile compound. In order to reduce the number of these additional parameters, we hypothesised that all volatile compounds have a single overall apparent coefficient $K_{p/m}$ that remains constant regardless of the texture of the matrix. It assumes that all volatiles compounds have the same partition coefficient and that they react in a similar way during the fragmentation and analysis in the mass spectrometer. However $K_{p/m}$ seems not susceptible to vary for many orders of magnitude. Atlan, Trelea, Saint-Eve, Souchon, and Latrille (2006) found that air-water partition coefficients for 12 volatile organic compounds (at 25 °C) are ranging from $10^{-4}$ to $10^{-2}$ (2 order of magnitudes, approximately). This assumption allows us to initiate the modelling work. During this work, it was not possible to adjust the experimental data correctly for the furfuryl alcohol and pyrazine compounds with only one value of apparent quantitation coefficient $K_{p/m}$.

The poor fit of these two compounds may be related to an incomplete or erroneous reaction scheme. Since the experimental data are not all quantitated, it is not possible to decide. Furthermore, these two compounds are the most polar present in the reaction scheme. They may have a partition coefficient that
differs from the other volatiles. The feedback made us add two distinct apparent quantitation coefficients $K_{fme}$ and $K_{pyr}$ for furfuryl alcohol and pyrazine, respectively, while $K$ is the apparent quantitation coefficient associated with the other volatile compounds.

The dialysis experiments conducted by Hofmann (1998) have established that in "sugar–casein" systems, the colour development is almost exclusively determined by the formation of bonds between protein oligomers and low-molecular-weight products originating from the Maillard reaction. Brands, Wedzicha, and van Boekel (2002) studied the quantification of melanoïdins formed during the cooking of a glucose–casein and fructose–casein system. They showed that the molar extinction coefficient of melanoïdins formed is not significantly different in the two systems. The low-molecular-weight products resulting from the Maillard reaction from these two sugars are indeed very similar in nature. These are the products that will be fixed to the protein backbone to form melanoïdins. In this case, melanoïdins progressively enrich with dicarbonyl compounds and the carbon to nitrogen ratio of the melanoïdins will gradually decrease (Brands et al., 2002). By comparing the results obtained in these two studies, a first chromatic extinction coefficient (Cc) was estimated. It roughly connects the chromatic index to the number of dicarbonyl compounds (nc) that are bound to the caseins. $\Delta C = Cc \times nc$. This relationship is only valid in the region of the linear response of the colour index, i.e., for a colour index value lower than 30. In this case, the value of Cc is estimated in the order of $1.3 \times 10^9 \text{ mol}^{-1}$. Because the media in these two studies and our matrix are not strictly identical, it should be noted that this is again only a rough estimate and that further work of quantitation will be needed to take into account the specificities of the cheese medium. After the parameter identification procedure, a value for Cc of $10^8 \text{ mol}^{-1}$ gave a better fit to the data. The 13-fold decrease observed between the value of Cc estimated from the literature and the one obtained after identification might be due to lactose, that is the main sugar in the present study and also a disaccharide in contrast to the previous studies where glucose was used.

Table 1 shows the values obtained for the four apparent quantitation coefficients that were further used in our study. These values produce "pseudo-quantitated" data. It has been verified that the contents of the different constituents are compatible with the theoretical yields generally observed for the Maillard reaction, such as the one described by Cerny (2008).

### 3.3. Parameters initiation

Since the proposed reaction scheme is apparent and that all the reactions described are not necessarily elementary, the rate constants identified must also be considered as apparent. In particular, they will depend on the choice made for the four quantitation parameters. Under these conditions, we initially have no information on the relationship of the different rate constants with the temperature and we cannot presuppose realistic initial values to allow the initiation of the numerical simulations.

In order to determine the temperature dependence for each rate constant, we performed the parameter identification at constant temperatures from the data obtained at 80, 100, 120, 135 and 150 °C in 5 independent simulations. To do this, we neglected the effect of the temperature rise (the first 3 min and 30 s of the thermal treatments). This is justified for the kinetics conducted at 80 and 100 °C for which the duration of temperature increase (respectively 35 and 50 s) is insignificant compared to the duration of the total heating time (up to two hours). In contrast, in the case of the kinetics conducted at 120, 135 and 150 °C, this may lead to an underestimation of the velocity constants. Since the aim of these five simulations is only to determine the order of magnitude

<caption>**Table 1**</caption>

Adjusted values for the quantitation constants ($K$) used for the conversion of arbitrary surface units into moles of volatiles compounds (fm: furfuryl alcohol, pyr: pyrazine) and for the conversion of the chromatic index into moles of dicarbonyl compounds fixed on the casein backbone of melanoidins.

<table><thead><tr><td>Constant</td><td>Value</td></tr></thead><tbody><tr><td>K</td><td>10<sup>11</sup></td></tr><tr><td>K<sub>fm</sub></td><td>25 × K</td></tr><tr><td>K<sub>pyr</sub></td><td>1.5 × K</td></tr><tr><td>K<sub>c</sub></td><td>10<sup>8</sup></td></tr></tbody></table>

of the values used to initialize the parameters for the fitting procedure at variable temperature, we considered that this underestimation was acceptable. It will be corrected during the final fitting procedure conducted at variable temperature.

Table 2 shows the values obtained for the 14 rate constants after the fitting procedure conducted at each temperature. The absolute value of the Pearson's correlation coefficient ($r$) associated with the linear regression of the 14 rate constants as a function of $1/T$ show that these rate constants, except the one related to the formation of melanoïdins ($k_{14}$), follow Arrhenius's law.

The rate constant $k_{14}$ does not follow Arrhenius's law, probably due to some analytical difficulties. Indeed, the measurements of chromatic index are saturated in the case of very high temperature treatments. In this regard the quantification of melanoidins by spectrocolorimetry methods, as was assayed by Brands et al. (2002) would have been more appropriate, since it makes it possible to avoid the saturation of the detector through dissolving the samples before reading. If the treatments made at 150 °C are removed, it occurs that $k_{14}$ also follows Arrhenius's law. In order to keep the model as simple as possible, we therefore considered that all the rate constants follow Arrhenius's law.

It is noteworthy that the linearisation of the Arrhenius's law via the logarithm may cause a loss of the normal distribution of the gradient estimator around its mean value. In this case, the boundary estimations made on account of the Student table will not be accurate. However, this method provides an order of magnitude for the different parameters. To take into account this loss of normality, it was decided to add an extra margin of 50% for the definition of the boundaries used for the parameters initialization, for the fitting procedure made at variable temperature.

### 3.4. Parameter estimation at variable temperature

The adjustment was carried out by introducing Arrhenius's law for each rate constant. The adjustment of the model at all temperatures simultaneously has been extremely difficult, due to the relative instability of the model. This is probably due to the fact that the parameters are not truly independent, since the pre-exponential factor and the activation energy of a given rate constant are closely correlated by the Arrhenius's law itself and also to the stiff properties of the system. Under these conditions, it is not excluded that the proposed solution might be derived from a local minimum. In order to minimise the existing dependency between the pre-exponential factor and the activation energy, Van Boekel (2009a, chap. 8, 2009b, chap. 14) suggests the use of a reparameterised Arrhenius's law. In this case it becomes difficult to compare the values obtained for the identified parameters with those given in the literature, as these data are pretty limited and we found reparameterized values for the chemical reactions presented in our particular study.

Fig. 3 shows the fit of the experimental data for some compounds after the fitting procedure. A close examination shows that the lactose is not properly adjusted at 80 and 100 °C. This lack of fit may be associated with a locally inaccurate or incomplete reaction
<caption>**Table 2**</caption>

Values for the rate constants as determined at 80, 100, 120, 135 and 150 °C, and Pearson coefficient obtained for the linear regression model: ln(k) = -E<sub>A</sub>/R × t + ln(k<sub>0</sub>).

<table><thead><tr><th>Constant</th><th>80 °C</th><th>100 °C</th><th>120 °C</th><th>135 °C</th><th>150 °C</th><th>r</th></tr></thead><tbody><tr><td>k<sub>01</sub></td><td>1.95 × 10<sup>-5</sup></td><td>2.15 × 10<sup>-5</sup></td><td>2.00 × 10<sup>-4</sup></td><td>2.34 × 10<sup>-4</sup></td><td>5.06 × 10<sup>-4</sup></td><td>-0.933</td></tr><tr><td>k<sub>02</sub></td><td>3.03 × 10<sup>-8</sup></td><td>8.00 × 10<sup>-8</sup></td><td>2.00 × 10<sup>-6</sup></td><td>4.55 × 10<sup>-6</sup></td><td>1.37 × 10<sup>-5</sup></td><td>-0.975</td></tr><tr><td>k<sub>03</sub></td><td>3.74 × 10<sup>-5</sup></td><td>2.00 × 10<sup>-4</sup></td><td>2.82 × 10<sup>-4</sup></td><td>2.82 × 10<sup>-4</sup></td><td>6.26 × 10<sup>-4</sup></td><td>-0.915</td></tr><tr><td>k<sub>04</sub></td><td>7.27 × 10<sup>-6</sup></td><td>1.30 × 10<sup>-5</sup></td><td>1.00 × 10<sup>-3</sup></td><td>1.26 × 10<sup>-3</sup></td><td>2.61 × 10<sup>-3</sup></td><td>-0.929</td></tr><tr><td>k<sub>05</sub></td><td>3.43 × 10<sup>-4</sup></td><td>9.80 × 10<sup>-4</sup></td><td>3.00 × 10<sup>-3</sup></td><td>3.61 × 10<sup>-3</sup></td><td>7.97 × 10<sup>-3</sup></td><td>-0.989</td></tr><tr><td>k<sub>06</sub></td><td>1.64 × 10<sup>-8</sup></td><td>4.00 × 10<sup>-7</sup></td><td>5.00 × 10<sup>-6</sup></td><td>2.77 × 10<sup>-5</sup></td><td>1.23 × 10<sup>-3</sup></td><td>-0.981</td></tr><tr><td>k<sub>07</sub></td><td>8.46 × 10<sup>-5</sup></td><td>8.40 × 10<sup>-3</sup></td><td>1.79 × 10<sup>-3</sup></td><td>1.79 × 10<sup>-3</sup></td><td>7.00 × 10<sup>-3</sup></td><td>-0.941</td></tr><tr><td>k<sub>08</sub></td><td>3.43 × 10<sup>-3</sup></td><td>5.11 × 10<sup>-3</sup></td><td>1.00 × 10<sup>-2</sup></td><td>1.00 × 10<sup>-2</sup></td><td>1.00 × 10<sup>-2</sup></td><td>-0.919</td></tr><tr><td>k<sub>09</sub></td><td>1.29 × 10<sup>-9</sup></td><td>2.41 × 10<sup>-8</sup></td><td>3.96 × 10<sup>-6</sup></td><td>2.41 × 10<sup>-5</sup></td><td>5.44 × 10<sup>-5</sup></td><td>-0.981</td></tr><tr><td>k<sub>10</sub></td><td>4.22 × 10<sup>-7</sup></td><td>1.45 × 10<sup>-5</sup></td><td>1.62 × 10<sup>-4</sup></td><td>2.12 × 10<sup>-4</sup></td><td>7.18 × 10<sup>-4</sup></td><td>-0.966</td></tr><tr><td>k<sub>11</sub></td><td>3.76 × 10<sup>-7</sup></td><td>5.50 × 10<sup>-6</sup></td><td>1.17 × 10<sup>-4</sup></td><td>2.33 × 10<sup>-4</sup></td><td>2.72 × 10<sup>-4</sup></td><td>-0.959</td></tr><tr><td>k<sub>12</sub></td><td>3.38 × 10<sup>-5</sup></td><td>1.72 × 10<sup>-4</sup></td><td>3.49 × 10<sup>-4</sup></td><td>9.16 × 10<sup>-4</sup></td><td>2.18 × 10<sup>-3</sup></td><td>-0.992</td></tr><tr><td>k<sub>13</sub></td><td>1.55 × 10<sup>-7</sup></td><td>3.68 × 10<sup>-5</sup></td><td>1.00 × 10<sup>-3</sup></td><td>2.19 × 10<sup>-3</sup></td><td>1.00 × 10<sup>-2</sup></td><td>-0.967</td></tr><tr><td>k<sub>14</sub></td><td>2.78 × 10<sup>-5</sup></td><td>4.80 × 10<sup>-7</sup></td><td>2.23 × 10<sup>-4</sup></td><td>1.62 × 10<sup>-3</sup></td><td>3.48 × 10<sup>-3</sup></td><td>-0.676</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Fit of the experimental data for some of the constituents using a multi-response stoichiometric model for the parameter values given in Table 3.</figcaption></figure>

scheme. In order to validate the reaction pathways, that are really followed during the thermal treatments applied to the cheese medium it would be necessary to carry out experiments to follow the fate of carbon-13 labelled precursors or experiments with matrices doped with selected precursors.

Table 3 presents the adjusted values for each of the 28 parameters and the uncertainty associated. The calculation of the uncertainties associated with the parameters lead to very different values depending on the constituents. It must be remembered here that these results are obtained on a real food product undergoing

major thermo-chemical changes (temperature ranging from 20 to 150 °C and pH from 7.3 to 4) and not on a "study system" submitted to a single temperature at a given and buffered pH.

To reduce the existing uncertainty for the different parameters, the number of experimental points should be increased accordingly. From an experimental point of view, recent work iteratively selected the most relevant kinetic (in a dynamic state) to perform in order to minimise the number of experiments that are necessary to conduct before the obtention of a satisfactory estimation of the model parameters (Goujot, Meyer, & Courtois, 2012). Such works
<caption>**Table 3**</caption>

Values obtained for the pre-exponential factor ($k_0$) and the activation energy ($E_a$) after the fitting procedure conducted on the experimental data obtained at every temperature for the 14 rate constants.

<table><thead><tr><th>Constant</th><th>k<sub>0</sub></th><th>[k<sub>0</sub><sub>int</sub> k<sub>0</sub><sub>exp</sub>]</th><th>E<sub>a</sub></th><th>[E<sub>a</sub><sub>int</sub> E<sub>a</sub><sub>exp</sub>]</th></tr></thead><tbody><tr><td>1</td><td>3.15 × 10<sup>12</sup></td><td>[3.06 × 10<sup>12</sup>, 3.24 × 10<sup>12</sup>]</td><td>125.0</td><td>[124.9, 125.1]</td></tr><tr><td>2</td><td>2.20 × 10<sup>10</sup></td><td>[2.19 × 10<sup>10</sup>, 2.21 × 10<sup>10</sup>]</td><td>120.0</td><td>[119.9, 120.1]</td></tr><tr><td>3</td><td>2.50 × 10<sup>5</sup></td><td>[2.49 × 10<sup>5</sup>, 2.51 × 10<sup>5</sup>]</td><td>70.0</td><td>[69.9, 70.1]</td></tr><tr><td>4</td><td>1.00 × 10<sup>14</sup></td><td>[2.88 × 10<sup>13</sup>, 1.71 × 10<sup>14</sup>]</td><td>130.0</td><td>[104.4, 155.6]</td></tr><tr><td>5</td><td>9.00 × 10<sup>13</sup></td><td>[8.70 × 10<sup>13</sup>, 9.30 × 10<sup>13</sup>]</td><td>130.0</td><td>[102.9, 157.6]</td></tr><tr><td>6</td><td>1.10 × 10<sup>12</sup></td><td>[3.28 × 10<sup>11</sup>, 1.87 × 10<sup>12</sup>]</td><td>132.0</td><td>[106.4, 157.6]</td></tr><tr><td>7</td><td>3.40 × 10<sup>9</sup></td><td>[0, 3.60 × 10<sup>10</sup>]</td><td>96.0</td><td>[66.0, 126.0]</td></tr><tr><td>8</td><td>1.10 × 10<sup>12</sup></td><td>[0, 3.16 × 10<sup>12</sup>]</td><td>120.0</td><td>[0, 298.5]</td></tr><tr><td>9</td><td>2.00 × 10<sup>6</sup></td><td>[0, 6.31 × 10<sup>6</sup>]</td><td>80.0</td><td>[0, 265.6]</td></tr><tr><td>10</td><td>2.00 × 10<sup>2</sup></td><td>[0, 9.56 × 10<sup>2</sup>]</td><td>40.0</td><td>[32.2, 47.8]</td></tr><tr><td>11</td><td>5.00 × 10<sup>5</sup></td><td>[0, 3.53 × 10<sup>6</sup>]</td><td>80.0</td><td>[21.9, 138.1]</td></tr><tr><td>12</td><td>5.00 × 10<sup>5</sup></td><td>[0, 3.89 × 10<sup>6</sup>]</td><td>60.0</td><td>[13.3, 106.7]</td></tr><tr><td>13</td><td>1.20 × 10<sup>8</sup></td><td>[0, 9.39 × 10<sup>8</sup>]</td><td>80.0</td><td>[32.8, 127.2]</td></tr><tr><td>14</td><td>1.00 × 10<sup>2</sup></td><td>[0, 1.37 × 10<sup>4</sup>]</td><td>60.0</td><td>[0, 122.5]</td></tr></tbody></table>

make sense in the case of stoichiokinetic studies that require a large number of time-consuming experimental analyses. It is observed that the uncertainties are much greater for the rate constants $k_6$ to $k_{14}$ corresponding to the reactions for which the constituents are only semi-quantitated. In addition, some unmeasured key intermediates, such as 2-deoxyribose, 1- and 3-deoxyglucosone or 5-hydroxymethylfurfural, play the role of a "buffer" in the reaction scheme. Indeed, the rate constants associated with these components (8, 9 and 10) are in this case only adjusted as compared to the consumption of the previous and formation of the following components. In this sense, the quantitation of these components would provide more robustness to the reaction scheme and reduce the uncertainty of the various parameters identified.

Jousse et al. (2002) obtain an activation energy of 128.8 kJ·mol⁻¹ for the formation of the Amadori compound. The value we estimated, 125 kJ·mol⁻¹ is of the same order of magnitude. Berg & Van Boekel (1994) show that in milk subjected to heat treatment at temperature higher than 100 °C, the lactose is mostly degraded via the Lobry de Bruyn Alberda Van Ekenstein pathway rather than by the Maillard reaction. In our study, however, it seems that the two reactions are more or less balanced. A possible hypothesis is that the mixture of sodium polyphosphates used in order to sequester the calcium ions and thus allow the product to obtain the desired texture is also a catalyst of the Maillard reaction pathway. That would equilibrate the balance between these two reaction pathways. Bell (1997) highlights the role of phosphate as a catalyst of the Maillard reaction. One last interesting feature is the evolution of the Amadori compound. In fact, this compound is adjusted to 80, 100 and 120 °C. In contrast, our model minimises the decrease occurring at 135 and 150 °C. The explanation may be related to the pH decrease. In fact, the Maillard reaction is auto-inhibited by its pH as a decrease in the pH implicates that the balance between the free amino groups in their reactive form R-NH₂ and unreactive form R-NH₃⁺ is moved towards the non-reactive form. This decrease is mainly due to the formation of weak acids such as formic or acetic acid by the Maillard reaction and caramelization pathways. Therefore, less lactulosyllysine compounds might be formed at lower pH and the balance between its formation and consumption is in favour of the decrease of this intermediate. In our study, the initial pH of the cheese medium model is 5.85. It decreases to 5.2 and 5.1 after two hours of cooking at 80 and 100 °C, respectively. To 5 and 4.7 after 20 min at 120 and 135 °C and even to 4.3 after 10 min at 150 °C. All pH values were measured at 20 °C after previous cooling of the cheese medium and are therefore independent on the temperature of the product. To take into account the inhibition of the reaction by the pH, an

attempt was made to model the pH decrease by the formation of weak acids. In the case of heated milk, Berg & Van Boekel (1994) were able to show that most of the decrease in pH was related to the production of formic acid mainly through the Lobry de Bruyn Alberda Van Ekenstein pathway. For our study however the weak acids formed are not sufficient to explain the magnitude of the pH drop observed. Van Boekel (2009a, chap. 8, 2009b, chap. 14) suggests many other elements that can induce pH variations during the processing of a food matrix, such as the water content, the temperature, the changes in ionic equilibria for example. Further work will be needed to explain the pH shift observed during the heating of the cheese medium and to include its consequences to the content of reactive amino species in the multi-responses stoichiokinetic model. This will surely be one of the key parameters to get a better fit of the data using the observable reaction scheme proposed.

## 4. Conclusion and prospects

It had been possible to model with relative success the variations of the selected markers originating from the Maillard reaction with a multi-response stoichiokinetic model. This is, to our knowledge, one of the first times that the work is done for a real food medium submitted to both temperature and pH variations.

The results obtained provide access to the activation energy and pre-exponential factors associated with the reactions of formation of volatile compounds. However, given the limited number of available experimental data, the parameters estimated present a rather high uncertainty. These results highlight the importance of choosing the best experimental conditions to achieve a good identification of the different parameters. The implementation of recent experimental planning tools should be a good strategy to select the most relevant experiments in order to accurately identify the model parameters. (Goujot et al., 2012). These new tools can iteratively select the most appropriate kinetics to conduct in order to keep the number of trials (and therefore the analyses workload) as small as possible but require the preliminary knowledge of the mechanisms at work. Therefore this study could be transposed to the stoichiokinetic model proposed as a second step of the sequential approach to improve the accuracy of the parameters.

The model presented is not yet directly transferable toward industrial equipment as most of them use direct heating. However it indicates that the thermal settings already applied in industrial conditions are already well adapted to reduce the occurrence of "overcooked" defects as they tend to minimise the time spent at temperatures above 120 °C. This particular temperature has been found as the temperature above which odorous compounds originated from the Maillard reaction, including maltol and furaneol, are formed from the degradation of the lactulosyllysine in milk-related products. This result is in agreement with the one obtained from Berg & Van Boekel (1994).

Additional studies will be needed to complete this model by introducing the effect of parameters of interest such as pH, water activity, or the initial concentration of some specific constituents (that can already been brought from the already processed raw materials such as milk powders or cheeses). In order to adapt the developed model to the reactions that are occurring during storage, it will be necessary to confirm the reaction mechanisms at lower temperatures (between 4 and 60 °C) during long time periods (up to 6 months). A combined model of the thermal history of the product from its manufacture to its consumption could then be developed and applied to the stoichiokinetic multi-response modelling procedure in order to estimate the formation rate of the odorous compounds during the whole life of the processed cheese product.
Acknowledgments

This work had the financial support of the Agence Nationale de la Recherche (French National Research Agency) under the Programme National de Recherche en Alimentation et Nutrition Humaine project ANR-06-PNRA-023. The author is grateful to Pr André Lebert for comments and valuable help carefully proofreading the manuscript and to the anonymous reviewers for their nice suggestions.

References

Atlan, S., Trelea, I. C., Saint-Eve, A., Souchon, I., & Latrille, E. (2006). Processing gas chromatographic data and confidence interval calculation for partition coefficients determined by the phase ratio variation method. *Journal of Chromatography A, 1110*, 146-155.

Bell, L. N. (1997). Maillard reaction as influenced by buffer type and concentration. *Food Chemistry, 59*(1), 143-147.

Berg, H. E., & van Boekel, M. A. J. S. (1994). Degradation of lactose during heating of milk. 1. Reaction pathways. *Netherlands Milk and Dairy Journal, 48*, 157-175.

Bertrand, E., Machado-Maturana, E., Chevarin, C., Portanguen, S., Mercier, F., Tournayre, P., et al. (2011). Heat-induced volatiles and odour-active compounds in a model cheese. *International Dairy Journal, 21*(10), 806-814.

Bley, M., Johnson, M. E., & Olson, N. F. (1985a). Predictive test for the tendency of Cheddar cheese to brown after processing. *Journal of Dairy Science, 68*, 2517-2520.

Bley, M., Johnson, M. E., & Olson, N. F. (1985b). Factors affecting non enzymatic browning of process cheese. *Journal of Dairy Science, 68*, 555-561.

Brands, C. M. J., & van Boekel, M. A. J. S. (2001). Reactions of monosaccharides during heating of sugar-casein systems: building of a reaction network model. *Journal of Agricultural and Food Chemistry, 49*(10), 4667-4675.

Brands, C. M. J., & van Boekel, M. A. J. S. (2002). Kinetic modeling of reactions in heated monosaccharide-casein systems. *Journal of Agricultural and Food Chemistry, 50*(23), 6725-6739.

Brands, C. M. J., & van Boekel, M. A. J. S. (2003). Kinetic modelling of reactions in heated disaccharide-casein systems. *Food Chemistry, 83*(1), 13-26.

Brands, C. M. J., Wedzicha, B. L., & van Boekel, M. A. J. S. (2002). Quantification of melanoidin concentration in sugar-casein systems. *Journal of Agricultural and Food Chemistry, 50*(5), 1178-1183.

Caric, B. (2000). *Processed-cheese in encyclopedia of food science and technology* (Second ed.). John Wiley & Sons Inc.

Cerny, C. (2008). The aroma side of the Maillard reaction. *Annals of the New York Academy of Sciences, 1126*(1), 66-71.

Deleris, I., Atlan, S., Souchon, I., Marin, M., & Trelea, L. C. (2008). An experimental device to determine the apparent diffusivities of aroma compounds. *Journal of Food Engineering, 85*(2), 232-242.

Goujot, D., Meyer, X. M., & Courtois, F. (2012). Identification of a rice drying model with an improved sequential optimal design of experiments. *Journal of Process Control, 22*, 95-107.

Hofmann, T. (1998). Studies on melanoidin-type colorants generated from the Maillard reaction of protein-bound lysine and furan-2-carboxaldehyde, chemical characterisation of a red coloured domain. *European Food Research and Technology, 206*(4), 251-258.

Jousse, F., Agterof, W., Jongen, T., Koolschijn, M., Visser, A., & Vreeker, R. (2002). Flavor release from cooking oil during heating. *Journal of Food Science, 67*(8), 2987-2996.

Kapoor, R., & Metzger, L. E. (2008). Process Cheese: Scientific and Technological Aspects; A Review. *Comprehensive Reviews in Food Science and Food Safety, 7*(2), 194-214.

Kondjoyan, N., & Berdagué, J. L. (1996). A compilation of relative retention indices for the analysis of aromatic compounds. *Saint Genès Champanelle, France: édition du Laboratoire Flaveur*.

Lee, S. K., Buwalda, R. J., Euston, S. R., Foegeding, E. A., & McKenna, A. B. (2003). Changes in the rheology and microstructure of processed cheese during cooking. *Lebensmittel-Wissenschaft Und-Technologie-Food Science and Technology, 36*(3), 339-345.

Martins, S. I. F. S., & van Boekel, M. A. J. S. (2005). A kinetic model for the glucose/glycine Maillard reaction pathways. *Food Chemistry, 90*(1-2), 257-269.

Panouille, M., Durand, D., Nicolai, T., Larquet, E., & Boisset, N. (2005). Aggregation and gelation of micellar casein particles. *Journal of Colloid and Interface Science, 287*(1), 85-93.

Parker, J. K. (2013). The kinetics of thermal generation of flavour. *Journal of the science of food and agriculture, 93*, 197-208.

Parker, J. K., Balagiannis, D. P., Higley, J., Smith, G., Wedzicha, B. L., & Mottram, D. S. (2012). Kinetic model for the formation of acrylamide during the finish-frying of commercial French fries. *Journal of Agricultural and Food Chemistry, 60*(36), 9321-9331.

Pischetsrieder, Monika, & Henle, Thomas (2010). Glycation products in infant formulas: chemical, analytical and physiological aspects. *Amino Acids, 42*(4), 1111-1118.

Pitkowski, A., Nicolai, T., & Durand, D. (2008). Scattering and turbidity study of the dissociation of casein by calcium chelation. *Biomacromolecules, 9*, 369-375.

Rakowsha, J., Haftka, R. T., & Watson, L. T. (1991). An active set algorithm for tracing parametrized optima. *Structural Optimization, 3*(1), 29-44.

Rocklin, R., & Pohl, C. A. (1983). Determination of carbohydrates by anion exchange chromatography with pulsed amperometric detection. *Journal of Liquid Chromatography, 6*(9), 1577-1590.

Saint-Eve, A., Juteau, A., Atlan, S., Martin, N., & Souchon, I. (2006). Complex viscosity induced by protein composition variation influences the aroma release of flavored stirred yogurt. *Journal of Agricultural and Food Chemistry, 54*(11), 3997-4004.

Samavati, V. (2013). Multivariate-parameter optimization of aroma compound release from carbohydrate-oil-protein model emulsions. *Carbohydrate polymers, 98*(2), 1667-1676.

Shampine, L. F., & Reichelt, M. W. (1997). The MATLAB ODE suite. *SIAM Journal on Scientific Computing, 18*(1), 1-22.

Smuda, Mareen, & Glomb, Marcus A. (2013). Fragmentation pathways during Maillard-induced carbohydrate degradation. *Journal of Agricultural and Food Chemistry, 6*(43), 10198-10208.

van Boekel, M. A. J. S. (2002). On the use of the Weibull model to describe thermal inactivation of microbial vegetative cells. *International Journal of Food Microbiology, 74*(1-2), 139-159.

van Boekel, M. A. J. S. (2009a). Multiresponse kinetic modeling of chemical reactions. *Kinetic modeling of reactions in foods. Boca Raton, FL: Taylor & Francis*.

van Boekel, M. A. J. S. (2009b). Modelling the food matrix. *Kinetic modeling of reactions in foods. Boca Raton, FL: Taylor & Francis*.

Yaylayan, V. A. (1997). Classification of the Maillard reaction: a conceptual approach. *Trends in Food Science and Technology, 8*(1), 13-18.

Yaylayan, V. A. (2003). Recent advances in the chemistry of Strecker degradation and Amadori rearrangement: implications to aroma and color formation. *Food Science and Technology Research, 9*(1), 1-6.