Sebastien Gaucel

Ioan Cristian Trelea

Steven Le Feunteun

Yan Li

David Julian Mcclements

Correction to New Mathematical Model for Interpreting pH-Stat Digestion Profiles: Impact of Lipid Droplet Characteristics on in Vitro Digestibility

) should be 2/3 instead of 3/2, leading to several changes in subsequent equations. In this note, we derive the corrected equations and compare them to the original ones, which show similar dynamical behavior.

Assumptions and equations from the original article

In order to make this note self-sufficient, we first recall the assumptions and equations from the original article. Numbering was kept identical to the original article, using an additional prime symbol ( ' ) for the corrected equations.

Assuming that two free fatty acid (FFA) molecules are produced per triacylglycerol molecule, the percentage of released FFA, %𝐹𝐹𝐴, is given by

%𝐹𝐹𝐴 = 100 × � 𝑉 𝑁𝑎𝑂𝐻 × 𝑚 𝑁𝑎𝑂𝐻 × 𝑀 𝑙𝑖𝑝𝑖𝑑 2 × 𝑤 𝑙𝑖𝑝𝑖𝑑 � (1)
where 𝑉 𝑁𝑎𝑂𝐻 is the volume of added sodium hydroxide solution (m 3 ), 𝑚 𝑁𝑎𝑂𝐻 is its molarity (kmol.m -3 ), 𝑤 𝑙𝑖𝑝𝑖𝑑 is the initial mass of the triacylglycerol oil (kg) and 𝑀 𝑙𝑖𝑝𝑖𝑑 is its molar mass (kg.kmol -1 ).

Assumption 1 The rate of FFA release from lipid droplets is proportional to the droplets surface area per unit volume of the emulsion, S d (m 2 .m -3 )

𝑑𝑚 𝑑 𝑑𝑡 = -𝑘 𝑆 𝑑 𝑉 𝐸 (2)
where 𝑚 𝑑 (𝑡) is the number of moles of digestible FFA in the droplets, 𝑑𝑚 𝑑 /𝑑𝑡 is the number of moles of FFA leaving the droplets per unit time, 𝑘 the rate constant (kmol.s -1 .m -2 ) and 𝑉 𝐸 the volume of the emulsion (m 3 ).

Remark: it is important to notice that equation (2) has a physical meaning if and only if 𝑚 𝑑 ≥ 0.

Assumption 2 The total number of lipid droplets remains constant and their diameter 𝑑(𝑡) decreases when the products of the lipolysis leave the droplet. The specific surface area of the droplets per unit volume of the emulsion is then given by

𝑆 𝑑 (𝑡) = 𝑛 0 𝜋 𝑑 2 (𝑡) (3)
where 𝑛 0 is the number of droplets per unit volume of emulsion (m -3 ). 𝑛 0 is computed from the initial state of the emulsion

𝑛 0 = 6 𝜙 0 𝜋 𝑑 0 3 (4)
where 𝜙 0 is the initial droplet volume fraction and 𝑑 0 the initial droplet diameter (m).

Assuming that all of the digested lipids leave the droplets, the ratio of the droplet volume during digestion to the initial droplet volume is given by

𝑉 𝑑 (𝑡) 𝑉 0 = 1 6 𝜋 𝑑 3 (𝑡) 1 6 𝜋 𝑑 0 3 = 𝑚 𝑑 𝑚 0 (5)
where 𝑚 0 is the initial number of moles of FFA in the droplets. Hence

𝑑(𝑡) = 𝑑 0 � 𝑚 𝑑 𝑚 0 3 (6)
Remark: from equation ( 6), 𝑑(𝑡) could become negative if the nonnegativity of 𝑚 𝑑 (𝑡) is not ensured, which would lead to inconsistent results.

From this point, a discrepancy occurs between the article and the actual solution, due to a misprint in the exponent. Indeed, substituting equations (3) to (6) in equation (2), the article gives

𝑑𝑚 𝑑 𝑑𝑡 = -𝑘 6 𝜙 0 𝑑 0 � 𝑚 𝑑 𝑚 0 � 𝟑 𝟐 � 𝑉 𝐸 (7)
while one should obtain

𝑑𝑚 𝑑 𝑑𝑡 = -𝑘 6 𝜙 0 𝑑 0 � 𝑚 𝑑 𝑚 0 � 𝟐 𝟑 � 𝑉 𝐸 (7')
Note the modification of the exponent between equations ( 7) and (7'): 3/2 versus 2/3.

Analytical solutions

Equation from the original article Corrected equation

𝑑𝑚 𝑑 𝑑𝑡 = -𝑘 6 𝜙 0 𝑑 0 � 𝑚 𝑑 𝑚 0 � 𝟑 𝟐 � 𝑉 𝐸 𝑑𝑚 𝑑 𝑑𝑡 = -𝑘 6 𝜙 0 𝑑 0 � 𝑚 𝑑 𝑚 0 � 𝟐 𝟑 � 𝑉 𝐸
The analytical solution writes

𝑚 𝑑 (𝑡) = � 3 𝑘 𝜙 0 𝑑 0 𝑚 0 3/2 𝑉 𝐸 𝑡 + 𝑚 0 -1/2 � -2 (8) 𝑚 𝑑 (𝑡) = �𝑚 0 1/3 - 2 𝑘 𝜙 0 𝑑 0 𝑚 0 2/3 𝑉 𝐸 𝑡� 3 (8')
which leads to

Φ = 1 - 𝑚 𝑑 𝑚 0 = 1 -�1 + 3 𝑘 𝜙 0 𝑉 𝐸 𝑡 𝑑 0 𝑚 0 � -2 (9) Φ = 1 - 𝑚 𝑑 𝑚 0 = 1 -�1 - 2 𝑘 𝜙 0 𝑉 𝐸 𝑡 𝑑 0 𝑚 0 � 3 (9')
Rewriting 𝜙 0 as

𝜙 0 = 𝑤 𝑑 𝑉 𝐸 𝜌 = 𝑚 0 𝑀 2 𝑉 𝐸 𝜌
with 𝜌 the density of the triacylglycerol oil (kg.m -3 ), it leads to

Φ = 𝜙 𝑚𝑎𝑥 �1 -�1 + 3 𝑘𝑀 2 𝑑 0 𝜌 𝑡� -2 � (10) Φ = 𝜙 𝑚𝑎𝑥 �1 -�1 - 𝑘 𝑀 𝑑 0 𝜌 𝑡� 3 � (10')
where 𝜙 𝑚𝑎𝑥 is the asymptotic maximal percentage of released FFA.

Dynamical behaviors of the solutions

A parameter estimation has been performed on equation (10') in order to fit the equation ( 10). We added a condition of nonnegativity for 𝑚 𝑑 in equation (10') which then became

Φ = 𝜙 𝑚𝑎𝑥 �1 -𝑚𝑎𝑥 �0, �1 - 𝑘 𝑀 𝑑 0 𝜌 𝑡� 3 �� (11')
The results are depicted in Figure 1 and show similar dynamics for both solutions (10) and (11').

Qualitatively, both solutions have a similar shape but the corrected one reaches a constant plateau more quickly. If nonnegativity is not preserved, equation (10') shows a totally different behavior. Indeed, mathematically, this solution is monotonically increasing, with an inflection point, but goes to infinity (Figure 2). In contrast, the original solution, equation (10), does not need an explicit condition of nonnegativity.

Conclusion

Both original and corrected equations are able to fit in vitro lipolysis data reasonably well. The original equation (10) reproduces the observed trend of experimental data somewhat better than the corrected equation (11') and has no need to explicitly enforce the nonnegativity of number of moles (or mass) to avoid unphysical results. Equation (10) does not correspond to the physical hypotheses described in the original article, however.

Figure 1 .

 1 Figure 1. Comparison between original (continuous line) and corrected (dashed line) solutions. The corrected solution included a condition for m d ≥0.

Figure 2 .

 2 Figure 2. Comparison between original (continuous line) and corrected (dashed line) solutions when nonnegativity of m d is not explicitly preserved.


Marie Noelle Leclercq-Perlat

Mariette Sicard

Nathalie Perrot
email: perlat@grignon.inra.fr

Ioan-Cristian C Trelea

Daniel Picque

Georges Corrieu

M.-N Leclercq-Perlat

M.-N Leclercq-Perlat

Temperature and relative humidity influence the ripening descriptors of Camembert-type cheeses throughout ripening

Keywords: descriptors of cheese ripening, temperature, relative humidity

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.

INTRODUCTION

Among environmental parameters affecting the Camembert-type cheese ripening, temperature, relative humidity (RH), and gaseous atmosphere composition of ripening chambers have the most effects on sensorial properties [START_REF] Choisy | The biochemistry of ripening[END_REF][START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF][START_REF] Weissenfluh | The effect of environmental conditions in the ripening room on the growth of Penicillium camemberti and the quality of Camembert cheese[END_REF].

Theses authors have highlighted the predominant role of temperature and have reported that low increase of temperature can accelerate cheese ripening but it can also produce off-flavor components or inconsumable products. This is also in agreement with the studies about Cheddar cheeses by [START_REF] Hannon | Application of descriptive sensory analysis and key chemical indices to assess the impact of elevated ripening temperatures on the acceleration of Cheddar cheese ripening[END_REF] and Reggianito Argentino cheeses by [START_REF] Sihufe | The influence of ripening temperature on the sensory characteristics of Reggianito Argentino cheese[END_REF]. To our knowledge, any systematic works about ripening conditions of Camembert-type cheeses have been realized on the influence of ripening temperature and RH on the cheese sensorial qualities. Consequently, it is very difficult to have a trend about the influence of temperature and RH on the sensorial qualities of Camembert-type cheeses throughout their ripening. However, it is admitted that sensorial properties like surface coating appearance, surface color or odor, are basic because they determine the choice and preferences of consumers [START_REF] Dufossé | Spectrocolorimetry in the CIE L*a*b* color as useful for monitoring the ripening process and the qualities of PDO red smear soft cheeses[END_REF].

According to [START_REF] Trystram | Computerized control systems in the food industry[END_REF], managing sensory properties at the manufacturing stage with the aim of controlling them is particularly difficult due to a few numbers of sensors available.

Knowledges about the links between microbiological and physicochemical phenomena and the changes in cheese sensorial properties are still lacking. According to [START_REF] Perrot | Decision support system design using the operator skill to control cheese ripening-application of the fuz-zy symbolic approach[END_REF] and [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF], operators often play an important role and, in practice, interact with automatic systems to assess the sensorial properties of the soft cheeses and to adjust the process. According to [START_REF] Lemoine | L'affinage, une étape décisive pour la qualité du fromage[END_REF], operator assessment and reasoning still play a major role in the cheese ripening process, especially concerning sensorial property control. In factories, the cheese-makers may use to instrumental measurements and empirical sensorial perceptions at the same time [START_REF] Perrot | Decision support system design using the operator skill to control cheese ripening-application of the fuz-zy symbolic approach[END_REF][START_REF] Picque | Effects of atmospheric composition on respiratory behavior, weight loss, and appearance of camembert-type cheeses during chamber ripening[END_REF]. Moreover, Bonaïti et al. (2004) have shown that temperature and RH had also an influence on creamy underrind thickness and on its consistency. Indeed, at 16°C under 100 % RH the underrind of a smear cheese became completely runny after 20 days of ripening. Recently, [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF] have shown that it is possible to study Camembert ripening process by macroscopic evolutions of cheeses evaluated from an expert's point of view. These authors established a correlation of 76 % between the microbiological, physicochemical, and biochemical data and the sensory phases measured according to expert knowledge.

The aim of this study was to evaluate the effects of temperature and relative humidity (RH) on developments of sensorial properties in the Camembert-type cheeses throughout their ripening (from d0 to d40) for the same set of cheese-making runs. These descriptors were studied in association to the microbial and physicochemical evolutions occurring throughout ripening.

MATERIAL AND METHODS

Camembert cheese-making and ripening

The microorganisms used (Kluyveromyces marxianus, Geotrichum candidum, Penicillium camemberti, and Brevibacterium aurantiacum) and the cheese-making have been described in detail previously by [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. Surface mould soft cheeses (100 cheeses per cheese-making run, each weighing 300 ± 20 g) were manufactured in a sterile environment. The cheeses were aseptically transferred to a previously sterilized ripening chamber. After 24 h at periodically renewed atmosphere. To that end, on d0, the ripening chamber, which contained air (CO 2 = 0%), was sealed. Given that during respiration the CO 2 concentration increased and when the CO 2 concentration reached 0.5%), it was automatically decreased to 0.1% by injecting humid sterile air [START_REF] Picque | Effects of atmospheric composition on respiratory behavior, weight loss, and appearance of camembert-type cheeses during chamber ripening[END_REF]. The cheeses were turned on the fifth day. On d13, they were kept for 24 h under 12 °C and 85 % RH. On d14, the cheeses were wrapped into a reference wrap (CdL, Amcor Flexibles, Barbézieux, France) and ripened at 4°C until d40.

Sensorial analyses performed on the cheeses

The sensorial analyses were performed on each cheese sample under the same light and the same temperature (14 ± 1 °C). For that, each cheese was put into a glass box with a cover and the sensory evaluation by the assessors was made immediately.

According to [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF] methodology, five expert assessors for analyzing each cheese were sufficient. Each expert assessed each cheese, taking into account five descriptors generally used by the cheese-makers, and scored on a scale from 1 to 5 (Table 1). Six sessions were devoted to training and evaluation before the validation of panel's performance [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF]. To determine the progress of Camembert-type cheese ripening, the appearance of the cheese (P. camemberti coating density and its surface covering) was assessed as previously described by Bonaïti et al. (2004) and [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF]. P. camemberti appearance was assessed in terms of mycelium density and uniformity, varying from 1 (less than 10 % of cheese surface were covered by P. camemberti mycelium) to 5 (more than 90 % of surface were covered by P. camemberti). In addition, for each cheese sampling, the overall odor, rind color under P. camemberti coating, underrind consistency (noted C UR ), and core hardness (noted H core ) on d14 and on d40 were considered. The higher differences of each sensorial descriptor between two panelists were lower than 0.2. The overall cheese odor varied from fresh cheese (score 1) to ammoniac (score 5). The typical odor of Camembert corresponded to the score of 3. The cheese color under the P. camemberti coating varied from white (score 1) to dark brown (score 5). The typical color (cream) of a Camembert corresponded to a score of 3. The creamy underrind thickness is a main factor allowing Camembert ripening according to Bonaïti et al. (2004) and Leclercq-Perlat et al (2004a, 2006[START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF], but this descriptor does not allow cheese-makers to define completely the quality of ripening due to the destruction of cheeses. Indeed, in practice, cheese-makers use another descriptor: the underrind consistency (C UR ) which characterizes the texture of ripened part of cheese as highlighted by [START_REF] Sicard | Toward the integration of expert knowledge and instrumental data to control food processes: Application to Camembert-type cheese ripening[END_REF]. The creamy underrind thickness (T UR , mm) was measured with a decimeter with an accuracy of 0.2 mm. The underrind consistency (Picque et al., 2011) varied from dried (score 1) to very runny (score 5). The core hardness, characterizing the texture of the core part, was estimated by the ease with which the cheese was ground with a mortar and pestle. This hardness varied from a very soft (score 1) to very hard (score 5).

To simplify the analysis, two characteristic days of ripening were only considered: the wrapping day (on d14) from which the cheeses were marketed and the end of ripening (on d40) corresponded to the end of shelf-life of the cheeses.

Experimental design

The effects of temperature and RH on the cheese sensory descriptors were examined using a two-factor, three-level complete factorial experimental design (3²) as described by [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. The 9 combinations of temperature and RH are shown in Tables 2 and3.

The levels of each factor used were 8, 12, and 16 °C for temperature and 88, 92, and 98 % for RH. These levels were chosen according to the range of temperature and RH generally used for Camembert cheese ripening. The 9 combinations of temperature and RH are shown in Tables 2 and3 (left part). Because the experiment took place over a lengthy period, the runs corresponding to the central point of the experimental design (12°C, 92 % RH) were quintupled (run 6 to 10). The runs under 8°C and 98 % RH (runs 3 and 4) and those under 16°C and 88 % RH (runs 12 and 13) were simply duplicated. The remaining runs were done just once.

Statistical analyses

For the runs made under the same conditions (runs 3 and 4; runs 6 to 10, and runs 12 and 13), sensorial descriptors were compared by a 2 factor (time and trial) analysis of variance (ANOVA). The hypotheses examined were the equality of the trials and the absence of any interactions. This test was significant for the risk α (α = 1 -p(F obs <F crit )) ≤ 0.01). To control the equality of trials, the test power (1-β), representing the risk of a false interpretation, was determined according to Mann and Whitney's method [START_REF] Cohen | A power primer[END_REF]. This power risk must be higher than 0.9 to consider that the interpretation of this test is accurate.

The general linear model within Statistica Software was performed to calculate quadratic models and to determine the influence of the two factors (θ, RH) on each sensorial descriptor. The nonsignificant terms were omitted one by one, using the procedure Stepwise Backwards (Statistica), consequently only the terms significant at 99 % of confidence level (p-value < 0.01) were considered. The three-dimensional response surfaces of each descriptor in relation to θ and RH factors were plotted to illustrate the main and interactive effects. These effects can be linear, quadratic and/or interactive.

RESULTS AND DISCUSSION

Whatever the run, on d0, P. camemberti appearance (absence of visual mycelium), overall odor (fresh cheese) and rind color (white) scored 1. Core hardness scored 3 (medium). For these descriptors, the standard deviations were less than 0.2. The underrind was non-existent.

Description of ripening under 12°C and 92 % RH

Statistical run reproducibility.

The reproducibility of cheese ripening descriptors was studied for runs 6 to 10 carried out under RH 92% and 12°C (Table 2). The hypotheses of equality of the means were greatly satisfactory (1 -α = 0.99 and the test power (1-β) = 0.94). Although the run period had lasted for three years, the runs carried out in the same conditions were statistically identical.

The time (t PCm in d), for which the first mycelia were observed in the rind or edge of cheeses, defined the growth of P. camemberti mycelium and the appearance of P. camemberti (noted N PC ) covering was defined as mycelium density and uniformity. Whatever the central point runs (6 to 10), the P. camemberti mycelia became visible on d5. From d14 to d40, all the cheeses were covered with a dense and white thick mycelium, giving the highest notation (Table 2).

On d14, the overall odor and the rind color were, respectively, "camembert" and cream while on d40 they scored, respectively, "camembert-ammoniac" and dark cream (Table 3). From d14 to d40, the core hardness was medium and the underrind consistency was qualified to soft. The underrind thickness (T UR) was close to 3.3 mm on d14 and formed all the cheese thickness on d40 (Table 3). For the 5 runs carried out under the standard ripening, all the descriptors were statistically identical with 99 % confidence levels.

Influence of temperature and RH on ripening descriptors

The scores of ripening descriptors are reported in Table 2 for P. camemberti appearance and Table 3 for the other descriptors.

Appearance of Penicillium camemberti

From the descriptors previously defined (t PCm and N PC ), t PCm was independent to RH and only related to temperature (Table 2). On d14, the appearance of P. camemberti mycelium (N PC14 ) was a function of temperature (linear and quadratic terms), RH (linear term) and temperature and RH (interactive term), and on d40 (N PC40 ) it was related to temperature and RH (linear and quadratic terms) (Table 4). The three-dimensional response surfaces of appearance of P. camemberti (N PC ) on d14 and d40 are given in Figures 1A and1B. The N PC14 surface response (Fig. 1A) shows a main effect of temperature. Indeed, from 8°C up to about 12°C, N PC14 increased with temperature and decreased afterwards (quadratic effect). The relative humidity effect was relatively smaller and depended on the temperature, pointing out an interactive effect: N PC14 increased with RH at low temperature and decreased at high temperature. The N PC40 surface response (Fig. 1B) shows the importance of the quadratic effects of temperature and RH: N PC40 reached its maximum value for a temperature near 12°C and a RH near 93%.

Ripening in chamber at 16°C under 98 % of RH promoted P. camemberti sporulation as previously shown by [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. [START_REF] Lenoir | La maîtrise du bioréacteur Fromage[END_REF] have pointed out that the optimal mycelium growth is obtained under 90 -94 % RH for a temperature between 10 and 14°C. Moreover, P. camemberti mycelium is not fond of higher RH [START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF]. Indeed, white and down coating of the surface of Camembert disappears when the RH of ripening room is higher than 95 % for a temperature of 12-14°C even if these conditions last only 2 or 3 days. However, neither P. camemberti visual notation nor P. camemberti spore counts give a real idea for mycelium growth and development, the first due to the limitation of human eyes (saturation when P. camemberti coating was complete) and the second due to the destruction of mycelium. Indeed, P. camemberti sporulation was accelerated throughout the ripening by the temperature [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF][START_REF] Leclercq-Perlat | Dynamics of Penicillium camemberti growth quantified by real-time PCR on Camembert-type cheeses at different conditions of temperature and relative humidity[END_REF]. Higher the temperature was, earlier its growth occurred and earlier its concentration became maximal. Moreover, at 16 °C and 98 % RH, P. camemberti mycelium appeared damaged due to disruption of the equilibrium between this mold and G. candidum in favor to this yeast while P. camemberti spore concentration was higher.

Influence of temperature and RH on ripening descriptors

Rind odor. On d14, the rind odor was musty at 8°C, Camembert at 12°C and Camembertammonia at 16°C (Table 3). This descriptor was related to temperature (individual and quadratic terms), RH (individual term), and temperature and RH (interactive term) (Table 4). In addition, the odor uniformly increased with temperature between 8 and 16°C while the RH effect was almost negligible at 8°C and positive at 16°C (Fig. 1C). On d40, the rind odor remained musty at 8°C under 88 and 92 % but became "camembert" under 98%. It was "camembert-ammoniac" at 12°C whatever the RH and at 16°C and 88 %. It was "ammoniac" at 16°C under 92 and 98 %.

This descriptor was a function of temperature (individual and quadratic terms), and RH (linear term) (Table 4). The odor increased with temperature and RH had a constant positive linear effect whatever the temperature (Fig. 1D).

An increasing of the temperature induced a better growth of the ripening microorganisms, mainly G. candidum and P. camemberti mycelium and when the temperature was higher than 10°C B. aurantiacum maximal exponential growth rate took place on d15 (12°C and 92 % RH, just after the wrapping) or on d12.6 at 16°C and 88 % RH [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. A better expression of these enzymatic activities can also be induced by RH [START_REF] Choisy | The biochemistry of ripening[END_REF].

Moreover, an increasing of temperature and RH can allow an increasing of aroma compound production and an increasing in the variety of flavor compounds [START_REF] Le Quéré | Cheese flavor[END_REF][START_REF] Mcsweeney | Biochemical pathways for the production of flavor compounds in cheese during ripening: a review[END_REF][START_REF] Sablé | Current knowledge of soft cheeses' flavor and related compounds[END_REF].

Rind color. On d14, the rind color increased with the temperature from white, corresponding to a cheese with a very little ripening (8°C under 88 and 92 % or 12°C under 88 %) to brown (16°C under 92 and 98 %). On d40, this color rose from white (8°C under 88 or 92 %) to dark brown, corresponding to an unacceptable color for a Camembert cheese (16°C under 92 or 98 %). On d14, rind color was a function of temperature and RH (individual and quadratic terms) while on d40, it was only related to temperature and RH (individual terms) (Table 4). The fig. 1E shows the importance of quadratic terms on the color on d14 while the fig. 1F shows that on 40 the color uniformly increased with temperature and RH (linear effects).

The rind color depends directly on production of pigments by Brevibacterium species [START_REF] Dufossé | Assessment of the coloring strength of Brevibacte-rium linens strains: spectrocolorimetry versus total carotenoid extraction/quantifycation[END_REF]. At 8°C whatever the RH or at 12°C and 88 % RH, a delay to the yeast installation and growths [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF] and P. camemberti mycelium was observed as shown by t PCm (Table 2). These main consequences were a delay to rind deacidification and to B. aurantiacum growth [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. Indeed, on d40, these authors showed that 1) at 8°C under 92 or 98 %, the bacterium concentrations were the weakest (between 6x10 6 and 2x10 8 CFU per g). On the contrary, at 16°C under 92 or 98 %, the rind deacidification became earlier (from d2.5) and B. aurantiacum grew better, reaching 8x10 9 CFU/g; and 2) at 16°C under 98 % RH, P. camemberti mycelium disappeared after d8 [START_REF] Leclercq-Perlat | Dynamics of Penicillium camemberti growth quantified by real-time PCR on Camembert-type cheeses at different conditions of temperature and relative humidity[END_REF] while the mold sporulation increased. These results highlighted the dependence of RH on rind color due to a better development of B. aurantiacum and a bad down coating of P. camemberti mycelium when the RH was higher than 95 %. Moreover, indirect interactions between yeasts, molds and B. aurantiacum can allow a better development of color (Leclercq-Perlat et al., 2004b;[START_REF] Lecocq | Importance de l'association entre Geotrichum candidum et Brevibacterium linens pendant l'affinage des fromages à croûte lavée[END_REF], this could be due to an important proteolysis [START_REF] Bockelmann | Color development by red-smear surface bacteria[END_REF].

Underrind thickness. On d14, T UR,14 rose around 50 % when the temperature increased from 8 °C to 12°C while it rose between 32 % and 40 % when temperature increased from 12°C to 16°C. T UR,14 increased from 10 % (8-12°C) to 15 % (16°C) when the RH increased from 88 to 92 % while it rose from 9 % to 17 % when the RH increased from 92 to 98 %. These results are in agreement with the ones of Gomes et al. (1998) and Gomes and Malcata (1998). Indeed, these authors have shown that temperature had a more important effect than RH for different proteolysis levels. T UR,14 (in mm) was related to temperature (linear term) and temperature and RH (interactive term) (Table 4) and T UR,14 increased along the temperature axis from 8 to 16°C and the temperature effect was stronger at higher RH (Fig. 2A).

On d40, at 8°C whatever the RH, the underrind thickness (T UR,40 in mm) reached around onethird of cheese thickness (Table 3). This can be explained by a slowdown of all ripening microorganism growths and as a result by their enzymatic reactions due to the excessive drying of cheese [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. Indeed, Lesage-Meessen and Cahagnier (1998) have highlighted that when the water activity, strongly linked up to RH during ripening [START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF], is progressively lowered the microbial growth is repressed. Whatever the temperature, for the ripening carried out under 88 % RH the underrind thickness did not increase as well as the ripening made under 92 or 98 % RH. For the other ripening conditions, all the thickness of cheese was ripened. This descriptor was related to temperature (quadratic term) and RH (linear and quadratic term) and their interactive term. These effects were shown by Fig. 2B which exhibited a maximum towards 14 ± 1°C and 93 ± 1 %.

Underrind consistency and core hardness

Underrind consistency (C UR ). On d14 or d40 at 8 °C and 88 % RH (run 1) the underrind was very dried (Table 3). Whatever the temperature under 88 % RH the score of underrind consistency was lower than the ones obtained for the two other RH. Only the central points (runs 6 to 10) had a soft and creamy underrind typical of a Camembert cheese (Table 3).

Whatever the RH at 16°C the cheeses did not have an underrind consistency acceptable because the underrind had runny edges (88 % RH) or completely runny (92 and 98 % RH). The underrind consistency which showed almost linear increases along the temperature and RH axes in the considered range (prevailing linear effects without interaction), confirmed the influences of temperature and RH an d14, (Fig. 2C) as well as on d40 (Fig. 2D). Moreover, the temperature effect was relatively stronger, especially on d14. Whatever the ripening time, the underrind consistency was correlated with temperature and RH (linear terms).

Core hardness. On d14 under 88 % RH the core was hard (at 8°C) or very hard (at 12 °C or 16 °C) (Table 3). Whatever the temperature under 92 % or 98 %, the hardness remained medium On d40, under 88 % RH, it remained hard at 8 °C or very hard at 12 °C but it scored medium at 16°C. At 12°C or 16°C under 92 or 98 %, all the cheese thickness was ripened and the core hardness cannot be scored by the assessors.

On d14, the core hardness was a function of temperature (individual term) and RH (interactive term) (Fig. 2E). This descriptor reached a minimum for a temperature that depended on the RH (interaction): it was around 13°C for 98 % RH and around 11°C for 88% RH as shown in fig. 1K.

Core hardness decreased with RH increase and the RH effect was slightly stronger for higher temperatures. On d40, it was impossible to determine a relation between core hardness and the ripening conditions because all the cheese was ripened.

The underrind consistency and the core hardness are two descriptors drawing the texture of cheeses. The cheese texture mainly depends on its pH, the entirety of caseins and lipids, and the moisture [START_REF] Lawrence | Texture development during cheese ripening[END_REF]. Indeed, whatever the ripening conditions, the pH of underrind was close to the one of rind [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. This can be explained by the lactate and ammonium concentration gradients which exist between the core and the rind (Leclercq-Perlat et al., 2004a;Vassal et al., 1986). The cheese pH is also defined by the production of microbial enzymes and theirs activities [START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF]. Moreover, whatever the temperature under 88 % RH, there was no soft and creamy underrind due to 1) an excessive drying of cheese, 2) a delay of K. marxianus and G. candidum growths, 3) a lower growth of P. camemberti and B. aurantiacum [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF][START_REF] Leclercq-Perlat | Dynamics of Penicillium camemberti growth quantified by real-time PCR on Camembert-type cheeses at different conditions of temperature and relative humidity[END_REF], and 4) lower enzymatic activities as shown by overall rind odor. Indeed, at 16°C and 92 % RH and under 98 % RH at 12

or 16 °C, the cheese underrind was completely runny, noticed an over-ripening. With the rapid increase of pH and its consequence on the texture, the underrind "liquefaction" can be explained by the casein breakdown and enzymatic activities associated due to a better and swift microorganism growths [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF][START_REF] Leclercq-Perlat | Dynamics of Penicillium camemberti growth quantified by real-time PCR on Camembert-type cheeses at different conditions of temperature and relative humidity[END_REF]. Indeed, the hydrolyze reactions which take place in cheese involve the water fixation and in a cheese core, the breakdown products of caseins are largely water soluble and proteolytic activities in cheese is mainly determined by the level of moisture, temperature of ripening, and changes in pH throughout the ripening [START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF][START_REF] Lawrence | Texture development during cheese ripening[END_REF]. Moreover, the changes in underrind consistency can be explained by enzymatic activities of surface ripening flora though the enzymes remained near their production place, i.e. the rind [START_REF] Noomen | The role of the surface flora in the softening of cheeses with a low initial pH[END_REF]. K. marxianus has only an endocellular proteolytic activity [START_REF] Klein | Peptidase activity of four yeast species frequently encountered in dairy products: comparison with seven dairy bacteria[END_REF], G. candidum strongly contributes to the proteolysis by its huge enzymatic system [START_REF] Boutrou | Interest in Geotrichum candidum for cheese technology[END_REF][START_REF] Boutrou | Contribution for Geotrichum candidum to the proteolysis of soft cheese[END_REF], P. camemberti also owns a complex proteolytic system which is more active when the rind pH is up to 6.0 [START_REF] Choisy | The biochemistry of ripening[END_REF]Lenoir et al., 2005), and B. aurantiacum is also well-known for its important proteolytic activities [START_REF] Lecocq | Importance de l'association entre Geotrichum candidum et Brevibacterium linens pendant l'affinage des fromages à croûte lavée[END_REF][START_REF] Rattray | Aspects of enzymology and biochemical properties of B. linens relevant to cheese ripening: a review[END_REF]. Under 98 % RH at 12°C and under 92 or 98 % RH at 16°C, [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF] have shown that 1) the death of K. marxianus cells began earlier and was more intensive than under the other ripening conditions; 2) G. candidum growth also began earlier (1-2 days), reaching its maximal concentration faster; 3) P. camemberti spore concentrations were the highest, reaching around 5.5x10 6 spores per g on d40 at 16°C; and 4) under 98 % RH at 12 or 16 °C, B. aurantiacum also grew quickly, reaching a concentration higher than 1x10 9 CFU per g. Moreover, a synergy has been found between G. candidum and P. camemberti, the first allowed an easier assimilation of peptides and amino acids as carbon sources by P. camemberti, involving a more intense proteolysis in the underrind of cheeses [START_REF] Aziza | Evidences for synergistic effects of Geotrichum candidum on Penicillium camemberti growing in cheese juice[END_REF][START_REF] Aziza | Commensalism during submerged mixed culture of Geotrichum candidum and Penicillium camemberti on glutamate and lactate[END_REF]. The underrind consistency is also related to lipid breakdown mainly due to G. candidum and P. camemberti [START_REF] Lenoir | La maîtrise du bioréacteur Fromage[END_REF] while K. marxianus and B. aurantiacum have no significant lipolytic activities [START_REF] Leclercq-Perlat | Controlled Production of Camembert-Type Cheeses: Part III. Free fatty acid evolutions[END_REF]. Lipolytic enzyme activities mainly depend on temperature and cheese moisture which is dependent on RH [START_REF] Choisy | The biochemistry of ripening[END_REF][START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF][START_REF] Ramet | Comparing ripening technology of the various types of cheese[END_REF]. Indeed, a low RH led to much dehydration of the cheese rind by water losses Water was vaporized from the surface toward the atmosphere, while at the same time water diffused from the core toward the surface [START_REF] Simal | Water and salt diffusion during cheese ripening: Effect of the external and internal resistances to mass transfer[END_REF]Bonaïti et al., 2004). This dehydration increased with the temperature: under 88 % RH at 16°C, the dry matter (DM) on d14 and on d40 was the most important, 47 % and 53 ± 1 % respectively, leading the stopping of ripening. However, whatever the RH at 8°C, the core remained hard on d40 although its DM (42-44 %) did not differ between the other ripening conditions. This could be explained by a slowdown of enzymatic reactions occurring in the cheese probably due to the lower temperature. By consequence of higher enzymatic activities and a lower dehydration of cheese, under 98 % RH at 12°C or 16°C, the underrind lost its softening to become liquid. These two factors involved a cheese moisture variable in relation to conditions of the ripening. Whatever the temperature under 88 % RH, the water cheese mass losses and the DM of the core were the highest [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF]. Under the same RH, the water mass loss increased from 54 % at 8 °C (88 % RH) and 68 % at 16 °C (92 or 98 % RH) when the temperature increased from 8°C to 16°C. This showed cheese moisture importance and, by consequence, the importance of RH on core hardness.

CONCLUSIONS

Because of all cheeses were wrapped with the same film and ripened under wrapping in the same conditions, the conditions of ripening in room govern the Camembert-type cheeses happening during their wrapping. In this work, the cheese properties assessed by cheese ripening experts were highly dependent on temperature and relative humidity used throughout the ripening. All sensorial descriptor changes (rind odor and color, underrind thickness, underrind consistency, and core harness) can be explained by microorganism growth, pH, carbon substrate diffusion and assimilation, cheese moisture as well as microbial enzymatic activities. All these factors were also linked to temperature and RH used throughout the ripening. Under 88% RH at 8°C, on d40, all sensorial descriptors scored the worst by the assessors: the cheese was too dry, its odor and its color were similar to the ones of the cheeses defectively ripened, the underrind was the driest, and the core was the hardest. On the contrary, under 98 % RH at 16 °C, the odor and the color of the cheeses were ammonia and dark brown, the underrind represented all the thickness of the cheeses but it was completely runny. These descriptors noticed an over-ripened cheese. As previously shown for microbial and physicochemical kinetic descriptors [START_REF] Leclercq-Perlat | Temperature and relative humidity influence the microbial and physicochemical characteristics of Camemberttype cheese ripening[END_REF], the best ripening conditions to obtain an optimum between sensorial qualities and cheese ripening properties appear to be 13 °C and 94 % RH. Such a study will provide main elements 1) to develop a mechanistic model allowing monitoring the Camembert-type cheese ripening process in relation to the ripening conditions, 2) to slow down or to accelerate the ripening process in relation to market demand and, 3) to correct some flaws. T UR , underrind thickness in mm (standard deviation = 0.5 mm).

6

The overall cheese odor varied from fresh cheese (score 1) to ammoniac (score 5). The cheese 7 color under the P. camemberti coating varied from white (score of 1) to dark brown (score 5).

8

The underrind consistency (noted C UR ) varied from a dried (score 1) to a very runny underrind 9 (score 5). The cheese core hardness (noted H) was estimated by the ease with which the 10 cheese can be ground with a mortar and pestle and varied from very hard (score of 1) to very 11 soft (score 5). -, the core was totally ripened. Table 4. Best-fit equations for the effects of temperature and relative humidity on Camembert-type cheese ripening descriptors. All factors of the equations for temperature (θ, °C) and relative humidity (RH, %) are given: *, p-value value lesser than 0.01 or **, p-value lesser than 0.001. R, determinant coefficient. df, degree of freedom. SE, standard error. p-value for the relation. C UR -16 (± 2)** + 0.33 (± 0.02)** x θ + 0.17 (± 0.02)** x RH 0.96 12 0.294 < 10

-7

H core NS N PC , appearance notation of Penicillium camemberti. Overall rind cheese odor varied from fresh cheese (score 1) to ammoniac (score 5).

The cheese color under the P. camemberti coating varied from white (score 1) to dark brown (score 5). The cheese core hardness (Hard core ) was estimated by the ease with which the cheese could be ground with a mortar and pestle: from very soft (score 1) to very hard (score 5). T UR , underrind thickness (mm). The underrind consistency (C UR ) varied from a dried (score 1) to a very runny underrind (score 5). NS, none significant.

  41)* -0.7 (± 0.02)* x θ + 3.3 (± 0.9)* x RH + 0.042 (± 0.007)** x θ² -0.017 (± 0.005)* x RH² 0.98 10 0.195 1.16 x 10 -7 T UR (mm) -1.8 (± 0.3)** -0.7 (± 0.1)** x θ + 0.012 (± 0.001)** x θ x RH 0± 2)** + 0.41 (± 0.05)** x θ + 0.10 (± 0.02)** x RH 0± 1)* + 0.057 (± 0.007)** x θ² -0.014 (± 0.002)** x RH x θ 0.91 12 0.384 1.17 x 10 -5 At the end of ripening (on d 40) N PC40 -365 (± 76)** + 3.1 (± 0.3)** x θ + 8.0 (± 0.9)** x RH -0.13 (± 0.01)** x θ² -0.042 (± 0.009)** x RH² 097)** + 10 (± 2)** x RH -0.08 (± 0.02)** x θ² -0.06 (± 0.01)** x RH² + 0.025 (± 0.004)** x θ x RH 0.97 10 0.453 1.01 x 10 -6

Figure 1 :

 1 Figure 1: Estimated response surface plots of ripening descriptors as a function of temperature (θ, °C) and RH (%). A) PC appearance notation on d14; B) PC appearance notation on d40; C) Odor on d14; D) Odor on d40; E) Color on d14;FA) Color on d40.

Figure 2 :

 2 Figure 2: Estimated response surface plots of ripening descriptors as a function of temperature (θ, °C) and RH (%). A) Underrind thickness (T UR ) on d14 (mm); B) Underrind thickness on d40; C) Underrind consistency (C UR ) on d14; D) Underrind consistency on d40; E) Core hardness on d14; F) Core hardness on d40.

Table 1 .

 1 Definition and score of the descriptors used to evaluate the Camembert-type cheeses

	1
	2	ripening progress.
						Scale and reference
		Descriptor	Notation	1	2	3	4	5
		PC appearance	N PC	< 10 %	25 %	50 %	75 %	> 90 %
		Overall odor	Odor	Fresh cheese	Musty	Camembert	Camembert ammoniac	Ammoniac
		Rind color	Color	White	Cream Dark cream Light brown Dark brown
		Underrind consistency	C UR	Very dried	Dried	Soft	Soft with runny edge	Runny
		Core hardness	H core	Very soft	Soft	Medium	Hard	Very hard
	3
	4

Table 2 .

 2 Descriptors of P. camemberti appearance obtained for each temperature (θ, °C) and

	6
		N°	θ (°C)	RH (%)	t PCm	N PC14	N PC40
		1	8	88	10	2	2
		2	8	92	10	2	3
		3	8	98	9	3	2
		4	8	98	9	3	2
		3 & 4	8	98	9	3	2
		5	12	88	6	4	5
		6	12	92	4.5	5	5
		7	12	92	5	5	5
		8	12	92	5	5	5
		9	12	92	4.5	5	5
		10	12	92	5	5	5
		6 to 10	12	92	5	5	5
		11	12	98	4	4	3
		12	16	88	3	4	2
		13	16	88	3	4	2
		12 & 13	16	88	3	4	2
		14	16	92	3	3	3
		15	16	98	3	2	1
	7
	8	t PCm , day for which P. camemberti mycelium is visually detected by the panel (standard
	9	deviation = 0.5 d). N PC14 and N PC40 , the appearance notation on d14 (wrapping) and on d40 (end
		of ripening), respectively.

5

relative humidity (RH, %) used during Camembert-type cheese ripening.

Table 3 .

 3 Scores of descriptors of Camembert-type cheese ripening on d14 (before wrapping)

	1
	2	and on d40 (end of ripening) obtained for each temperature (θ, °C) and relative humidity (RH,
	3	%).
	4
						On d14 (before wrapping)			On d40 (at the end of ripening)
	Run	θ (°C)	RH (%)	Odor	Color	T UR	C UR	H core	Odor	Color	T UR	C UR	H core
	1		8	88	2	1	1.3	1	4	2	1	9.5	1	4
	2		8	92	2	1	1.8	1	4	2	1	10.5	2	4
	3		8	98	2	1.5	2.0	2	3	3	2	10	3	4
	4		8	98	2	1.5	2.0	2	3	3	2	10	3	4
	3 & 4	8	98	2	1.5	2.0	2	3	3	2	10	3	4
	5		12	88	3	1	2.7	2	5	4	2	11	2	5
	6		12	92	3.5	2	3.3	3	3	4	3	14	3	-
	7		12	92	3.5	2	3.3	3	3	4	3	14	3	-
	8		12	92	3.5	2	3.5	3	3	4	3	14	3	-
	9		12	92	3.5	2	3.4	3	3	4	3	14	3	-
	10		12	92	3.5	2	3.3	3	3	4	3	14	3	-
	6 to 10	12	92	3.5	2	3.3	3	3	4	3	14	3	-
	11		12	98	3.5	2	4.1	3	2	4.5	3	14	4	-
	12		16	88	4	3	3.9	4	5	4	4	12.5	4	3
	13		16	88	4	3	3.9	4	5	4	4	12	4	3
	12 & 13	16	88	4	3	3.9	4	5	4	4	12	4	3
	14		16	92	5	4	5.7	5	4	5	5	14	5	-
	15		16	98	5	4	6.0	5	3	5	5	14	5	-
	5

ACKNOWLEDGMENTS

This work was supported by the French Research Ministry through a contract of National Program of Agronomic Research (PNRA). The authors are deeply grateful to F. Lecornué and T. Cattenoz (UMR GMPA) for their technical assistance and they are gratefully acknowledged the five assessors for their expertise.


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

Fernanda Fonseca

Stéphanie Passot

Denis Flick

A binary gas transport model improves the prediction of mass transfer in freeze-drying Short title: Binary gas transport model in freeze-drying

Keywords: lyophilization, dynamic model, mass transfer, convection, diffusion

Monitoring partial vapor pressure in the freeze-drying chamber is a cheap, global and non-intrusive way to assess the end of the primary drying stage.

Most existing dynamic freeze-drying models which predict this partial pressure describe mass transfer between the product and the condenser via a mass transfer resistance or a mass transfer coefficient. Experimental evidence suggests that such models can be significantly in error for some values of the sublimation flux, leading to physically inconsistent predictions and possibly incorrect assessment of primary drying termination, with potential risk of product damage if moving to secondary drying and increasing shelf temperature while some ice is still present.

Assuming a binary gas transport model for vapor and inert gas leads to improved and consistent predictions and explains the apparent variation of the mass transfer resistance with total pressure, shelf temperature and product sublimation area.

INTRODUCTION

Freeze-drying (lyophilization) is widely used for long term preservation of thermosensitive biological material and pharmaceuticals, such as proteins, vaccines, bacteria, mammalian cells and high quality food. [START_REF] Tang | Design of freeze-drying processes for pharmaceuticals: practical advice[END_REF][START_REF] Sadikoglu | Freeze-drying of pharmaceutical products: Research and development needs[END_REF] It removes water or an organic solvent in a way that minimizes the modification of the molecular structure of the active ingredient and creates a porous matrix with high rehydration properties. Freeze-drying remains a costly process, however, and lot of research was devoted to its optimization, often based on mathematical modeling. [START_REF] Sadikoglu | Freeze-drying of pharmaceutical products: Research and development needs[END_REF][START_REF] Sadikoglu | Mathematical modelling of the primary and secondary drying stages of bulk solution freeze-drying in trays: parameter estimation and model discrimination by comparison of theoretical results with experimental data[END_REF][START_REF] Sadikoglu | Optimal control of the primary and secondary drying stages of bulk solution freeze drying in trays[END_REF][START_REF] Sadikoglu | Optimal control of the primary drying stage of freeze drying of solutions in vials using variational calculus[END_REF][START_REF] Gan | Heating policies during the primary and secondary drying stages of the lyophilization process in vials: effects of the arrangement of vials in clusters of square and hexagonal arrays on trays[END_REF][START_REF] Chouvenc | Optimization of the freeze-drying cycle: a new model for pressure rise analysis[END_REF][START_REF] Boss | Freeze drying process: real time model and optimization[END_REF][START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF][START_REF] Velardi | Development of simplified models for the freezedrying process and investigation of the optimal operating conditions[END_REF][START_REF] Pisano | In-line optimization and control of an industrial freezedrying process for pharmaceuticals[END_REF][START_REF] Antelo | Toward optimal operation conditions of freeze-drying processes via a multilevel approach[END_REF] The process consists of three main steps: freezing, ice removal by sublimation (primary drying) and unfrozen water removal by desorption from the solid matrix (secondary drying).

The primary drying stage is usually the most time and energy consuming part of a freezedrying cycle. A lot of research has been devoted to the development of methods to precisely and consistently identify the end of the primary drying. Most of these methods rely on the monitoring of the vapor partial pressure in the freeze-drying chamber, either directly (moisture sensors such as aluminum oxide probes) and/or indirectly (Pirani gauge, mass spectrometer, Tunable Diode Laser Absorption Spectroscopy (TDLAS)). [START_REF] Gieseler | Evaluation of tunable diode laser absorption spectroscopy for in-process water vapor mass flux measurements during freeze drying[END_REF][START_REF] Patel | Determination of end point of primary drying in freeze-drying process control[END_REF] Moreover, abundant literature is available concerning modeling of primary freeze-drying for optimizing this step of the process off line and also more recently in real time. [START_REF] Pisano | In-line optimization and control of an industrial freezedrying process for pharmaceuticals[END_REF][START_REF] Antelo | Toward optimal operation conditions of freeze-drying processes via a multilevel approach[END_REF][START_REF] Fissore | In-line control of a freeze-drying process in vials[END_REF] Models often predict the product temperature and the sublimation flux. The rate of solvent vapor removal from the product is usually described as being governed by three barriers or resistances: resistance of the dried-product layer, resistance of the containers (vials and stoppers), and resistance of the chamber to condenser pathway. Since the gas phase is mainly vapor during the primary drying and gradients of concentrations are small, true diffusion, or flow under the influence of a concentration gradient, is often considered of minor importance and models focus on bulk flow mechanisms. [START_REF] Pikal | Heat and mass transfer in low-pressure gases: applications to freeze-drying[END_REF] The resistance of the dried product layer and of the containers has been studied both theoretically and experimentally by many groups. These studies present two limitations, however. Firstly, the chamber resistance is usually neglected because of measurement difficulties, or it is estimated by rather complex approaches involving Computational Fluid Dynamics (CFD) simulations. [START_REF] Alexeenko | Computational analysis of fluid dynamics in pharmaceutical freeze-drying[END_REF][START_REF] Rasetto | On the use of a dual-scale model to improve understanding of a pharmaceutical freeze-drying process[END_REF][START_REF] Ganguly | Modeling and measurements of water-vapor flow and icing at low pressures with application to pharmaceutical freeze-drying[END_REF][START_REF] Petitti | Cfd modelling of condensers for freezedrying processes[END_REF] Secondly, the partial vapor pressure is usually not predicted but taken as fixed and equal to the total chamber pressure, despite its importance for real time optimization and freeze-drying control. [START_REF] Sadikoglu | Mathematical modelling of the primary and secondary drying stages of bulk solution freeze-drying in trays: parameter estimation and model discrimination by comparison of theoretical results with experimental data[END_REF][START_REF] Velardi | Development of simplified models for the freezedrying process and investigation of the optimal operating conditions[END_REF][START_REF] Mascarenhas | A computational model for finite element analysis of the freeze-drying process[END_REF][START_REF] Sheehan | Modeling of the primary and secondary drying stages of the freeze drying of pharmaceutical products in vials: Numerical results obtained from the solution of a dynamic and spatially multi-dimensional lyophilization model for different operational policies[END_REF][START_REF] Song | A finite volume analysis of vacuum freeze drying processes of skim milk solution in trays and vials[END_REF][START_REF] Hottot | Freeze-drying of pharmaceutical proteins in vials: modeling of freezing and sublimation steps[END_REF] Our group previously developed a one-dimensional heat and mass transfer model, which can accurately represent both the primary and secondary drying stages and the gradual transition between them. [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] When adapting it to other freeze-dryers and to bacterial suspensions, however, it turned out that the determined chamber resistance strongly varied with total pressure and shelf temperature. Moreover, for a sufficiently high sublimation flux, e.g. generated by a high shelf temperature at the beginning of the primary drying, the model can predict a partial vapor pressure higher than the total one, which is obviously incorrect.

In the present study the chamber resistance was experimentally determined and a model of mass transfer between the freeze-drying chamber and the condenser is proposed for the solvent vapor and inert gas pair. The model is based on the mass transfer theory in binary gas mixtures and includes bulk flow and mutual diffusion terms. This model improved the prediction of the chamber vapor pressure and of the primary drying termination. It also accounted for apparent variations of the effective mass transfer resistance with total pressure, shelf temperature and product sublimation area.

MATERIALS AND METHODS

Determination of the mass transfer resistance (vapor flow) of the chamber to condenser pathway

Experiments were carried out on a LyoBeta special freeze-dryer (Telstar, Terrassa, Spain) equipped with 5 thermocouples, a capacitive manometer and a Pirani gauge. Purified water and bacterial cell suspensions were successively considered.

A stainless steel tray was filled with purified water to a thickness of 3 cm. The evaporative flow was determined using the gravimetric method. After freezing the water at -50°C (cooling rate of 0.6 °C/min), and sublimation time between 4 and 8 hours depending on the experimental conditions, the process was stopped and the tray was weighted to determine the amount of water removed and thus the evaporative flux (F V , in kg s -1 ). The mass transfer coefficient (k m , kg s -1 Pa -1 ), or, equivalently, resistance of the chamber to condenser pathway (r m , Pa s kg -1 ) was determined from Eq. 1:

Eq. 1

(

)

m V V V V m V r P P P P k F 2 1 2 1 - = - =
with P V1 , the water vapor partial pressure just above the product (Pa), estimated from the ice temperature and P V2 , the water vapor partial pressure at the condenser (Pa), estimated from the condenser temperature.

The effect of three process variables on the effective mass resistance (r m ) value was investigated: sublimation area (0.033, 0.074 and 0.15 m²), chamber pressure (10, 20, 40 and 60 Pa) and shelf temperature (-15, 0 and 15°C).

Freeze-drying experiments with bacterial suspensions

Lactic acid bacteria were produced by fermentation in controlled conditions of pH and temperature. [START_REF] Passot | Critical water activity and amorphous state for optimal preservation of lyophilized lactic acid bacteria[END_REF] After concentration, the bacterial cells were re-suspended in a 1:2 cells/protective medium ratio. The protective medium was composed of 200 g/L of sucrose and 0.15 M of NaCl. A stainless steel tray was filled with 450g of bacterial suspension and the following freeze-drying protocol was applied: freezing at -50°C (cooling rate of 0.6°C/min); primary drying at -20, 0 or 25°C; secondary drying at 25°C for 8 hours. The total chamber pressure was controlled at either 20 or 60Pa. Thermocouples were placed at the bottom of the product in the centre of the tray.

Estimation of partial vapor pressure from the Pirani gauge

The reading of a Pirani gauge in a pure gas is proportional to the pressure of the gas and to the molecular heat conductivity of the gas. [START_REF] Jennings | Lyophilization: Introduction and basic principles[END_REF] The main gases present in the freeze-drying chamber are solvent vapor (usually water) and an inert gas such as nitrogen or air. Let Eq. 2

N V λ λ α =
be the ratio of the molecular heat conductivities. This ratio is equal to 1.6 for the water vapor and nitrogen pair. If the Pirani gauge is calibrated with inert gas, nitrogen in the present case, then the reading in pure solvent vapor will be:

Eq. 3

V P V P P α =
In a gas mixture, the Pirani gauge will measure the sum of heat fluxes due to conduction in each gas. The reading will be: Eq. 4

N V P P P P + = α
In contrast, the capacitive sensor will give the true total pressure: Eq. 5

N V t P P P + =
From these two readings, one can determine the partial pressures of vapor and inert gas:

Eq. 6

MASS TRANSFER MODELS BETWEEN CHAMBER AND CONDENSER

Mass transfer resistance model

In most previously developed freeze-drying models, the mass flux between the freeze-drying chamber containing the product (location 1) and the condenser (location 2) was either ignored, considering the vapor pressure in the freeze-drying chamber as given, [START_REF] Sadikoglu | Mathematical modelling of the primary and secondary drying stages of bulk solution freeze-drying in trays: parameter estimation and model discrimination by comparison of theoretical results with experimental data[END_REF][START_REF] Velardi | Development of simplified models for the freezedrying process and investigation of the optimal operating conditions[END_REF][START_REF] Mascarenhas | A computational model for finite element analysis of the freeze-drying process[END_REF][START_REF] Sheehan | Modeling of the primary and secondary drying stages of the freeze drying of pharmaceutical products in vials: Numerical results obtained from the solution of a dynamic and spatially multi-dimensional lyophilization model for different operational policies[END_REF][START_REF] Song | A finite volume analysis of vacuum freeze drying processes of skim milk solution in trays and vials[END_REF] or expressed via a mass transfer coefficient (k m ) or, equivalently, via a mass transfer resistance (r m ), [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF][START_REF] Antelo | Toward optimal operation conditions of freeze-drying processes via a multilevel approach[END_REF][START_REF] George | Development and validation of heat and mass transfer models for freeze-drying of vegetable slices[END_REF][START_REF] Chouvenc | Optimization of pharmaceuticals freeze-drying cycles: characterization of annealing effects by the pressure rise analysis method[END_REF][START_REF] Hottot | Determination of mass and heat transfer parameters during freeze-drying cycles of pharmaceutical products[END_REF] as given by Eq. 1.

Note that in Eq. 1 the mass flux is proportional to the partial vapor pressure difference between the two locations. Vapor pressure P V2 is usually fixed by the condenser temperature, while the mass flux F V is mostly determined by the heat transfer flux from the temperaturecontrolled shelf to the product, i.e. by the shelf temperature and the product sublimation area.

If the mass transfer resistance was constant, as it is usually assumed, then the partial vapor pressure in the chamber P V1 would result from Eq. 1. If the mass flux was high enough, due for example to high shelf temperature, a model based on Eq. 1 could even give a vapor pressure higher that the total chamber pressure, which is obviously impossible. A more consistent mass transfer model is developed in the next section.

Binary gas transport model

Consider a mixture of two gases at low pressure; in the case of freeze-drying these gases would be solvent vapor, usually water, released from the product and captured by the condenser and the inert gas, usually nitrogen, used to control the total chamber pressure. To determine the relevant mass transport theory, the Knudsen number, defined as the ratio of the free mean path of the molecules and the characteristic dimension of the system, was estimated. In the considered conditions, the free mean path is of order of 0.5mm at 20Pa, while the typical diameter of the chamber to condenser duct is 100mm. Knudsen diffusion can thus be safely neglected for the chamber to condenser pathway considered here, since K n ≈ 5×10 -3 . Note that this is usually not the case in the porous product layer where the pore size can be much less than 0.5mm, but mass transport in the porous layer is out of the scope of the present study.

The classical theory of binary diffusion in gases states that the molar flux density of solvent vapor, in stationary coordinates, is given by: [START_REF] Bird | Transport phenomena[END_REF] Eq. 8

(

) ( )

V VN N V N V V V x D c c N N x N ∇ + - + =
where N V and N N are molar flux densities, c V and c N molar concentrations of vapor and nitrogen respectively, x V is the molar fraction of vapor and D VN is the mutual diffusion coefficient. The first term on the right-hand side of Eq. 8 results from the bulk motion of the fluid while the second one expresses the diffusion of vapor due to its molar fraction gradient.

Assuming pseudo-stationary conditions, compared to the duration of a freeze-drying cycle which is of one or several days, and assuming insignificant leaks in the freeze-drying chamber, the net inert gas flux is negligible compared to the vapor flux:

Eq. 9 Eq. 11

RT P RT P P c c t N V N V = + = + Eq. 12 V V V N M f =
With these substitutions, and assuming that the variation of total pressure between freezedrying chamber and condenser is small compared to the variation in vapour partial pressure, Eq. 8 becomes:

Eq. 13 ( )

V VN t V V t V P D P RT M P P f ∇ - = -
To solve Eq. 13, it is assumed that the vapour flux between the product and the condenser is essentially one-dimensional. This assumption is relatively sound for freeze-driers equipped with a condenser separated from the main freeze-drying chamber by a pipe, as is the case in the present study, but it is questionable when the product and the condenser are situated in the same cavity. With this assumption, Eq. 13 can be integrated between the freeze-drying chamber (location 1) and the condenser (location 2) giving:

Eq. 14 1 2 ln V t V t VN t eff V V V P P P P D P RT l AM Af F - - = =
where F V is the vapor mass flux, while A and l eff are the diffusion cross-section area and effective duct length respectively. The effective length accounts for the presence of inlet and outlet sections, bends, valve and other elements that can introduce additional resistances.

The mutual gas diffusion coefficient D VN is inversely proportional to the total pressure, thus the product P t D VN is constant. [START_REF] Bird | Transport phenomena[END_REF] For a given freeze-dryer geometry (l eff , A) a constant β can be introduced: Eq. 15

VN t V eff D P AM R l = β so that Eq. 16 1 2 ln 1 V t V t V P P P P T F - - = β
The dependence of the mass flux on the partial vapor pressures given by Eq. 16 appears quite different from Eq. 1. It can be shown, however, using a first order Taylor series expansion of the logarithm function, that Eq. 1 is an approximation of Eq. 16 in the case when the vapor pressures are much less than the total pressure: Eq. 17 It should be noted, however, than in freeze-drying the condition t V P P << 1 is usually not satisfied as the gas in the freeze-drying chamber during primary drying is mostly solvent vapor. The approximation given by Eq. 17 is thus not expected to be valid. This can explain the mentioned inconsistency with Eq. 1 and the apparent variation of the effective mass transfer resistance with total pressure and mass flux, as further discussed below. For now, note that Eq. 16 can accommodate arbitrarily large vapor fluxes without P V1 exceeding P t , provided that the difference 1 V t P P -becomes small enough, while remaining positive.

RESULTS AND DISCUSSION

Experimental determination of chamber to condenser mass transfer resistance

Table 1 summarizes the experimental values of the mass transfer resistance of chamber to condenser pathway measured for various process conditions.

The effective mass transfer resistance appeared to increase when increasing chamber pressure, decreasing shelf temperature and decreasing sublimation area. The mass transfer resistance was thus inversely correlated with the sublimation flux. Increasing the sublimation rate resulted in an apparent decrease of the effective mass transfer resistance of the chamber to condenser pathway. These observed variations support the idea that the mass transfer resistance is only an effective model parameter and a different perspective is needed to get a more physically meaningful description of the transfer phenomena. This is further discussed in the next sections.

Prediction of partial vapor pressure in the chamber by the two modeling approaches: resistance vs. binary gas transport

The predictions of the previously developed freeze-drying model, [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] adapted for bacterial suspension conditioned in trays, were compared to experimental measurements acquired in freeze-drying experiments performed as described in the "Materials and methods" section.

Two versions of the model were tested, one based on a mass transfer resistance between the chamber and the condenser (Eq. 1) and the other based on the binary gas transport assumption (Eq. 16). Details of the existing model [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] adaptation to the new freeze-dryer, container and product, as well as implementation details of the new chamber-to-condenser mass transport model based on Eq. 16 are given in the Appendix. The model parameters (notably k m and β) were determined by fitting the models to some of the experimental data obtained in bacterial suspension experiments (Appendix and Table 2) An example of the predictions of the two versions of the model are shown in Figure 1, for two experiments performed at 0°C shelf temperature/60Pa chamber pressure and 25°C/20Pa, respectively. These are validation experiments not used for model parameter identification. Dotted lines represent calculations based on the mass transfer resistance model and solid lines correspond to the binary gas transport model.

The experiment at 0°C/60Pa (Figure 1, left panels) illustrates a situation where the prediction of the primary drying termination differs by about 2 hours between the considered models.

The end of the primary drying was arbitrarily defined here as the moment when the sublimation flux (Figure 1C) decreased below 5% of its maximum value. The resistance model predicts the end of the primary drying at about 11h, in agreement with measured product temperature approaching the shelf temperature (Figure 1A), while the binary gas model predicts the end at bout 13h, in agreement with the decrease of vapor pressure in the chamber (Figure 1B). Note that product temperature measurement is a local one, potentially disturbed by the presence of the probe, while partial vapor pressure measurement is a global one, representative of the whole product in the freeze-dryer. Vapor pressure should thus be preferred to assess the end of the primary drying. [START_REF] Patel | Determination of end point of primary drying in freeze-drying process control[END_REF] It is also non intrusive and usually more conservative, minimizing the risk of premature temperature increase and product damage, if some ice was still present when starting the secondary drying step.

The experiment at 25°C/20Pa (Figure 1, right panels) illustrates a situation where both models predict the same end point of the primary drying (11h), but the vapor pressure predicted by the resistance model is completely in error. Compared to the previously described experiment (0°C/60Pa), in this case the heat transfer resistance from the shelf to the product is higher due to lower chamber pressure (lower contribution of the gas conduction), but this is nearly compensated by a higher shelf temperature, leading to a similar heat flux. The net result is a similar primary drying time and sublimation flux (Figure 1C) for both experimental conditions. As expected from Eq. 1, a similar sublimation flux implies a similar vapor pressure calculated with the resistance model, around 50Pa in this case (Figure 1B), but this value is 2.5 times higher than the total pressure, which is obviously incorrect. In contrast, the binary gas model based on Eq. 16 consistently predicts a vapor pressure close to, but slightly less than the total one, as physically expected.

Both these examples suggest that the binary gas transport model should be preferred to the resistance model. Similar trends were observed in all performed experiments with bacterial suspensions (data not shown). The properties of the proposed binary gas transport model are further examined in the following sections.

Effective mass transfer resistance of the chamber to condenser pathway

Since Eq. 1 and Eq. 14 are two different ways of expressing the same vapor flux based on different physical assumptions, these equations can be used to define an effective mass transfer resistance of chamber to condenser pathway as: Eq. 18

(
) It readily appears from Eq. 18 that this resistance is not constant but is an effective model coefficient that depends on the actual values of the total and vapor pressure.

1 2 2 1 1 2 2 1 ln ln 1 V t V t V V V t V t VN t V V V
The effective mass transfer resistance given by Eq. 18 is plotted in Figure 2 (solid line) as a function of the vapor pressure near the product (P V1 ), for two different total pressures (P t = 20 and 60Pa) used in our experiments. The vapor pressure near the condenser was fixed to a low value P V2 = 0.055 Pa, corresponding to a typical condenser temperature of -80°C. For comparison, the value of r m determined for the mass transfer resistance model in the same conditions is also given in Figure 2 (dotted line).

Firstly, Figure 2 indicates that the effective mass transfer resistance strongly increases with the total pressure. This was expected from the experimental data in Table 1 and confirms the usual practice of performing freeze-drying at low pressure to improve mass transfer.

Secondly, the effective mass transfer resistance predicted by the binary gas model approaches zero when the vapor pressure in the chamber approaches the total pressure. This was expected from the mathematical analysis of Eq. 18 and is consistent with the fact that partial vapor pressure in the chamber cannot exceed the total one, whatever the sublimation flux. Note that this is not the case with a constant mass transfer resistance, when the vapor pressure in the chamber can formally exceed the total pressure for a sufficiently high sublimation flux in Eq.

1.

Finally, Figure 2 shows that the value of the mass transfer resistance determined using the resistance model based on Eq. 1 is situated between the bounds given by Eq. 18, as expected.

The effective resistance given by the binary gas transport model is close to this value for vapor pressures very close to the total pressure. This means that the vapor pressure calculated by the binary gas model for the primary drying is most of the time very close to the total pressure, as observed in Figure 1 and well known from freeze-drying practice. This is further illustrated in the following section.

Sublimation flux and predicted vapor pressure in the chamber

In usual freeze-drying conditions the sublimation mass flux is the result of several process parameters (Table 1), such as total chamber pressure (that changes both heat and mass transfer), shelf temperature and sublimation interface area, etc. For the mass transfer resistance model given by Eq. 1, the predicted vapor pressure in the freeze-drying chamber varies linearly with the mass flux (dotted line in Figure 3) and can, as already mentioned, exceed the total pressure, for example in the case of a sufficiently high shelf temperature (Figure 1, right panels). In contrast, the vapor pressure predicted by the gas diffusion model given by Eq. 16 asymptotically tends towards the total pressure for high sublimation mass fluxes, which is physically consistent (solid lines in Figure 3). This also corroborates the well known fact that gas composition in the freeze-drying chamber is dominated by vapor during most of the primary drying stage, when sublimation fluxes are high (approaching 2×10 -5 kg s -1 in the considered experimental setting).

Sublimation flux and effective mass transfer resistance of the chamber to condenser pathway

Considering the effect of the sublimation mass flux in the binary gas transport model gives a common framework to understand the apparent variations of the effective mass transfer resistances observed in Table 1 for different shelf temperatures and sublimation interface areas. Increasing both these variables increases the sublimation flux and this, in turn, causes an apparent reduction of the effective mass transfer resistance (Figure 4). For sublimation fluxes approaching 2×10 -5 kg s -1 , representative of the primary drying in the considered experiments (Figure 1), the mass transfer resistance determined for the resistance model is intermediate between those calculated with the binary gas model for 20 and 60Pa total pressure. This reflects the compromise that the model fitting procedure finds to accommodate experiments performed at both these pressures, but the achieved compromise underestimates the vapor pressure at 60Pa and overestimates it at 20Pa, as already shown in Figure 1.

The binary gas transport model needs a unique constant parameter

As an additional difference between the two models, note that in the mass transfer resistance model the value of the resistance between the chamber and the condenser pathway is only an effective model parameter that appears to be a function of several process variables (Table 1): total pressure, shelf temperature and product sublimation area. Accounting for this would require varying the resistance in a tricky and ad-hoc manner. Considering a constant value of the mass transfer resistance (r m = 2.78×10 6 Pa s kg -1 in our case) leads to poor predictions and inconsistencies, such as calculated partial vapor pressure in the chamber exceeding the total one (Figure 1 and Figure 3). On the contrary, binary gas diffusion theory states that the product P t D VN is constant and therefore in the binary gas transport model the analogous parameter β, given by Eq. 15, is constant. All solid plots in Figures 1234were obtained with the same value of β = 2.66×10 3 s kg -1 K -1 .

In-process mass flux estimation and detection of the primary drying termination

Eq. 16 could be used in principle for obtaining vapor flow rate estimations based on total and partial pressure measurements, and thus act as an in-process mass flow meter. In practice, however, estimations of "high" vapor flow rates (i.e. when the partial vapor pressure approaches the total one) are expected to be inaccurate due to the vapor saturation of the chamber illustrated in Figure 3: a wide range of sublimation fluxes corresponds to almost identical partial vapor pressures, very close to the total chamber pressure. For example, with the considered values of the parameters (Appendix) at a total pressure of 60Pa, a vapor pressure of 59.99Pa corresponds to a sublimation flux of 1.14×10 -5 kg s -1 , while a vapor pressure of 58.99Pa corresponds to a sublimation flux of 0.53×10 -5 kg s -1 . In other words, an uncertainty of only 1Pa in the measurement of the vapor pressure would induce a factor 2 error on the flux value, which is clearly inappropriate for practical use.

On the contrary, Figure 3 suggests that the estimation of "low" vapor fluxes (i.e. when the partial vapor pressure is well below the total one) based on pressure measurements could be quite accurate. For example, vapor pressures of 30 and 29Pa correspond to vapor fluxes of 0.091×10 -5 kg s -1 and 0.086×10 -5 kg s -1 respectively; a measurement uncertainty of 1Pa would induce only about 5% error on the vapor flux.

Above considerations suggest that the proposed model is appropriate for a safe detection of the primary drying termination because a decrease of the vapor pressure well below the total one indicates a strong decrease of the sublimation flux. Continuing the above numerical examples, a vapor pressure in the chamber decreasing to the half of the total pressure corresponds to a factor 10 reduction of the sublimation flux. Of course, exact significance of "high" and "low" mass fluxes and the relevant decrease of vapor pressure depend on numerical values grouped in the parameter β defined by Eq. 15, and also involved in Eq. 16.

Finally, note that the proposed in-process flow estimation method based on total and partial pressure measurements in the chamber is complementary to the method described by Patel et al. [START_REF] Patel | Choked flow and importance of Mach I in freeze-drying process design[END_REF] , based on differential pressure measurement between the chamber and the condenser. Indeed, the former works well for low sublimation fluxes while the latter is expected to be more accurate when the sublimation flux, and hence the measured pressure difference, are large.

CONCLUSIONS

Results presented in this paper suggest that using one-dimensional mass transfer theory in binary gas mixtures to describe vapor transfer between the freeze-drying chamber and the condenser should be preferred to the mass transfer resistance approach. Advantages of the binary gas transport model include: (i) prediction of the vapor pressure in the freeze-drying chamber closer to the experimental one, leading to more reliable prediction of the end of the primary drying; (ii) calculated values of the vapor pressure physically consistent with the total pressure, whatever the sublimation vapor flux; (iii) need to determine a single model parameter for a given freeze-dryer geometry, independent on the operating conditions (total pressure, shelf temperature, ice sublimation area).

The main drawback of the mass diffusion model is the non-linear dependence between the vapor flux and vapor pressure, which can somewhat complicate model simulations. Solution of non-linear differential and algebraic equation sets is already part of most freeze-drying models, however, and the additional difficulty introduced by the diffusion model is expected to be minor, as confirmed by the authors' experience.

The assumption of one-dimensional mass transfer between the freeze-drying chamber and the condenser could be questioned in the future, especially for freeze-dryers with the condenser situated in the freeze-drying chamber. The challenge would be to end up with a reasonably simple closed-form relationship between mass flux and pressure, in order to maintain the overall model minimalism and short simulation time.

APPENDIX: EMBEDDING THE BINARY GAS TRANSPORT EQUATION IN A COMPLETE FREEZE-DRYING MODEL

The binary gas transport equation proposed in this paper (Eq. 16) describes the mass flux between the freeze-drying chamber and the condenser, in replacement of the more classical equation involving a mass transfer resistance or a mass transfer coefficient (Eq. 1). Both equations were tested in conjunction with an existing freeze-drying model, described in full detail previously. [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] In a first step, the parameters of the existing model [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] based on Eq. 1 were re-estimated to account for a different freeze-dryer (LyoBeta special freeze-dryer from Telstar, Terrassa, Spain, instead of SMH15 freeze dryer from Usifroid, Maurepas, France), a different container (metallic tray instead of glass vial) and a different product (lactic acid bacterial suspensions in sucrose medium instead of polyvinylpyrrolidone). Parameters related to the sorption isotherm and to the glass transition of the product were determined with the methodology described previously. [START_REF] Passot | Critical water activity and amorphous state for optimal preservation of lyophilized lactic acid bacteria[END_REF] Parameters related to the heat and mass transfer in the freeze-dryer and in the product were estimated by fitting the product temperature and vapor pressure predicted by the model to experimental measurements. Four experiments were used simultaneously for model parameter estimation, performed in the following combinations of shelf temperature and total chamber pressure: -20°C/20Pa, -20°C/60Pa, 25°C/20Pa and 25°C/60Pa. A genetic optimization algorithm (Global Optimization Toolbox for Matlab, Natick, MA) was used for parameter identification. Since the considered optimization algorithm is stochastic, the parameter set with the best fit from 8 independent optimization runs was selected. Most optimization runs consistently converged to similar sets of parameters, usually within ±10%.

Model parameters used in the present study are listed in Table 2.

In a second step, the original mass transfer equation between the chamber and the condenser (Eq. 1) was replaced by the newly proposed one (Eq. 16) and the parameter estimation procedure was repeated using the same data and the same algorithm. The estimated parameter values are also given in Table 2. As expected, the estimated values of the parameters are similar between the two models, with the exception of the parameter describing mass transfer between the freeze-drying chamber and the condenser (k m for the resistance model and β for the binary gas transport model), which have different physical meanings.

Since the newly introduced equation is nonlinear and the vapor pressure in the chamber is unknown a priori, the sublimation flux cannot be determined explicitly as in the resistance version of the model. [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] The vapor pressure in the chamber (P V1 ) was considered as an additional state variable and the corresponding mass balance equation in the freeze-drying chamber was added to the existing set of differential equations:

Eq. 19

V V V V F F dt dP RT V M - = ⋅ 31 1 1
Here V 1 is the freeze-drying chamber volume, F V is the chamber to condenser vapor flux given by Eq. 16 and F V31 is the sublimation vapor flux between the ice sublimation front (location 3) and the freeze-drying chamber (location 1), given by: Eq. 20

(
)

1 3 41 34 31 1 1 1 V V m m V P P k k F - + =
where k m34 and k m41 are the mass transfer coefficients between the sublimation front and the product top (location 4), and between the product top and the chamber, respectively. Thus, both mass fluxes involved in Eq. 19 can be written explicitly in terms of known (fixed or state) variables (P V1 , P V2 , and P V3 ).

Note that in the original freeze-drying model vapor accumulation in the freeze-drying chamber was neglected, based on a time scale analysis. [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF] This is formally equivalent to setting the chamber volume 0 1 ≈ V in Eq. 19, which would thus become an algebraic instead of a differential equation. This version of the model was also tested but not retained, due to occasional failures in solving the equation set during parameter identification. Observed difficulties are likely to be related to the specific differential-algebraic equations (DAE) solver available in Matlab (ode15s) which had numerical problems with some sets of parameters tested by the optimisation algorithm. Parameters specific for the product formulation, determined in separate experiments [START_REF] Passot | Critical water activity and amorphous state for optimal preservation of lyophilized lactic acid bacteria[END_REF] Sorption isotherm, wet basis (kg kg -1 ) ( )( ) Heat transfer coefficient by gas conduction, between shelf and product bottom (W m -2 K -1 Pa -1 ) 0.400 0.366

Heat conductivity between sublimation front and product top (W m -1 K -1 ) 0.0957 0.0996 Heat transfer coefficient between product top and chamber walls (W m -2 K -1 ) 3.79 3.03

Vapor conductivity between sublimation front and product top (kg s -1 m -1 Pa -1 ) 2.07×10 -8 1.98×10 -8 Vapor transfer coefficient between product top and chamber (kg s -1 m -2 Pa -1 ) 8.12×10 -5 3.68×10 -4 Vapor transfer coefficient between chamber and condenser k m (kg s -1 Pa -1 ) 3.60×10 -7 Not applicable

Vapor transfer parameter between chamber and condenser β (s kg -1 K -1 ) Not applicable 2.66×10 3

Residual frozen layer thickness for gradual transition between primary and secondary drying (m) 0.00173 0.00170

Effective chamber wall temperature (°C) 15.3 13.9

Complete model equations and exact definitions of the model parameters are given in ref. [START_REF] Trelea | An interactive tool for freeze-drying cycle optimisation including quality criteria[END_REF]

  Calculations were performed with Matlab™ 8 software (The MathWorks Inc., Natick, MA) equipped with the Statistics Toolbox and the Global Optimization Toolbox.

N

  Usual variables in freeze-drying are mass fluxes instead of molar fluxes and partial pressures instead of molar fractions and concentrations. Using the ideal gas law one gets:

  to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement CAFE n° KBBE-212754.

Figure 1 .

 1 Figure 1. Comparison between experimental measurements and simulations with the binary gas transport and resistance models, for two experiments with lactic acid bacteria suspensions (left panel 0°C/60Pa, right panel 25°C/20Pa). (A): Shelf and product temperatures, (B): Total and vapour pressures, (C): Sublimation fluxes. Bold: shelf temperature and total pressure, symbols: measurements, solid: binary gas model, dotted: resistance model, vertical line: end of primary drying.

Figure 2 .

 2 Figure 2. Variation of the effective mass transfer resistance with the vapor and total pressure.Numeric values are for lactic acid bacteria suspension. Solid: binary gas model, bold: total pressure 60Pa, thin: total pressure 20Pa, dotted: resistance model.

Figure 3 .

 3 Figure 3. Predicted vapor pressure as a function of the sublimation mass flux. Numeric values are for lactic acid bacteria suspension. Solid: binary gas transport model, bold: total pressure 60Pa, thin: total pressure 20Pa, dotted: resistance model.

Figure 4 .

 4 Figure 4. Apparent variation of the effective mass transfer resistance of the chamber to condenser pathway as a function of the sublimation mass flux. Numeric values are for lactic acid bacteria suspension. Solid: binary gas transport model, bold: total pressure 60Pa, thin: total pressure 20Pa, dotted: resistance model.

Figure 1 .Figure 2 . 1 ) 7 Figure 3 .Figure 4 .

 121734 Figure 1. Comparison between experimental measurements and simulations with the binary 606 gas transport and resistance models, for two experiments with lactic acid bacteria suspensions 607 (left panel 0°C/60Pa, right panel 25°C/20Pa). (A): Shelf and product temperatures, (B): Total 608 and vapour pressures, (C): Sublimation fluxes. Bold: shelf temperature and total pressure, 609 symbols: measurements, solid: binary gas model, dotted: resistance model, vertical line: end 610 of primary drying. 611

Table 1 .

 1 Chamber to condenser mass transfer resistance measured with pure ice.

	Total chamber pressure (Pa)	Shelf temperature (°C)	Sublimation interface area (m²)	Chamber to condenser resistance (10 6 Pa s/kg)
	10			0.65
	20 40	0	0.074	1.29 2.64
	60			3.72
		-15		4.31
	40	0	0.074	2.64
		15		1.61
			0.033	2.34
	20	0	0.074	1.29
			0.15	0.67

Table 2 .

 2 Values of the freeze-drying model parameters 600

	Parameters defined by the experimental setup
	Number of trays (-) Sublimation area (m 2 )	1 0.0745
	Mass of dry product (kg)	0.0736
	Mass of ice (kg)	0.3763
	Product height in the tray (m) Volume of the freeze-drying chamber (m 3 )	0.0067 0.202