*Article* **A Population Balance Model to Describe the Evolution of Sublethal Injury**

**Simen Akkermans 1,2,3, Davy Verheyen 1,2,3, Cindy Smet 1,2,3 and Jan F. M. Van Impe 1,2,3,\***


**Abstract:** The detection and quantification of sublethal injury (SI) of pathogenic microorganisms has become a common procedure when assessing the efficiency of microbial inactivation treatments. However, while a plethora of studies investigates SI in function of time, no suitable modelling procedure for SI data has been proposed thus far. In this study, a new SI model structure was developed that relies on existing microbial inactivation models. This model is based on the description of inactivation kinetics between the subpopulations of healthy, sublethally injured and dead cells. The model was validated by means of case studies on previously published results, modelled by different inactivation models, i.e., (i) log-linear inactivation; (ii) biphasic inactivation; and (iii) loglinear inactivation with tailing. Results were compared to those obtained by the traditional method that relies on calculating SI from independent inactivation models on non-selective and selective media. The log-linear inactivation case study demonstrated that the SI model is equivalent to the use of independent models when there can be no mistake in calculating SI. The biphasic inactivation case study illustrated how the SI model avoids unrealistic calculations of SI that would otherwise occur. The final case study on log-linear inactivation with tailing clarified that the SI model provides a more mechanistic description than the independent models, in this case allowing the reduction of the number of model parameters. As such, this paper provides a comprehensive overview of the potential and applications for the newly presented SI model.

**Keywords:** food safety; predictive microbiology; mathematical models; microbial inactivation; sublethal injury

#### **1. Introduction**

In predictive microbiology, the effects of intrinsic and/or extrinsic factors on the microbial behaviour (e.g., growth, inactivation) in foods are studied and quantified to predict the effect of varying environmental conditions on the microbial response. Accumulated knowledge on microbial behaviour is distilled into mathematical models, which can be integrated into user-friendly software tools. These tools can be used by food producers, governments and scientists to determine food safety and quality aspects [1,2].

Due to the health consequences related to the ingestion of pathogenic foodborne microorganisms, food industries design processing treatments to inactivate microorganisms that may be present in food products [3]. Exposing food products containing bacterial populations to an inactivation treatment (e.g., heating, irradiation, antimicrobial agents) leads to the formation of three bacterial subpopulations, i.e., (i) healthy cells that are uninjured; (ii) sublethally injured cells that can potentially still recover from their injuries; and (iii) dead cells that are completely inactivated/killed [4]. Sublethal injury (SI) is generally defined as "a consequence of exposure to a chemical or physical process that damages but does not kill a microorganism" [5–7]. The injury can involve a loss of the

**Citation:** Akkermans, S.; Verheyen, D.; Smet, C.; Van Impe, J.F.M. A Population Balance Model to Describe the Evolution of Sublethal Injury. *Foods* **2021**, *10*, 1674. https://doi.org/10.3390/foods 10071674

Academic Editors: Carlos Vilas, Míriam R. García and Jose A. Egea

Received: 17 June 2021 Accepted: 14 July 2021 Published: 20 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

permeability barrier in the cell wall and/or membrane (i.e., structural damage), damage to functional cell components, such as ribosomes and structural DNA (i.e., metabolic damage), or a combination of both types [8]. Sublethally injured cells are sensitive to selective components to which uninjured cells show resistance, making them unable to grow on selective media commonly used for the detection of foodborne pathogens in the food industry [9]. Consequently, the challenges related to the occurrence of sublethally injured cells in foods are twofold. The first challenge relates to food diagnostics, as the number of microorganisms is underestimated when solely enumerated by plating on selective media. As a second challenge, sublethally injured microorganisms might recover from their damage if exposed to optimal conditions for a sufficient amount of time, for example during food storage [10–13].

The degree of SI of a microbial population is traditionally assessed by the difference in plate counts on non-selective and selective media [7,8,14]. Non-selective media allow the growth of the total population of all viable culturable cells (i.e., uninjured and sublethally injured), while selective media only allow the growth of uninjured cells [15]. Hence, SI can be quantified as a percentage of the entire population by means of formulas similar to Equation (1), introduced by the likes of Busch et al. and Dykes [16,17].

$$\text{SI} = \left(1 - \frac{\text{Counts on selective medium}}{\text{Counts on nonselectronic medium}}\right) \cdot 100\% \tag{1}$$

Following the increased attention to SI in the last decades, the scientific literature experienced a rise in the number of predictive microbiology studies that include a representation of the SI evolution of microorganisms during inactivation. In most studies, SI is calculated at different time points, regardless of the inactivation behaviour of the cells [18–32]. To investigate the SI of cells during inactivation treatments in more depth, some researchers have calculated the SI evolution over time based on the modelled inactivation of the total culturable population and the healthy subpopulation. Apart from visually illustrating the SI evolution in a continuous manner, calculating the SI evolution as a function of the inactivation treatment time allows a comparison of total amounts of SI among different treatment conditions, e.g., using the time-averaged injured cells coefficient (TICC) [33].

SI cells should not be confused with so-called viable but non-culturable (VBNC) cells. VBNC cells are not regarded as dead cells because their cell membrane and genetic material is intact and they are metabolically active [34]. SI cells are not dead either. They have sustained injuries that cause them not to be viable on stress-containing selective media, but they are still viable on non-selective media. As such, the difference between VBNC and SI cells is that the former cannot be cultured on any media, while the latter can still be cultured on non-selective media [35]. The use of a combination of optimal rich media and selective media therefore allows differentiating between SI cells, on the one hand, and VBNC and dead cells on the other hand.

Thus far, different approaches for calculating the SI have been used with no consensus on a standard method having been reached. In the most common method, the SI evolution is calculated by using the non-log transformed model outputs corresponding to non-selective and selective media as input data for the SI equation [36–43]. The major disadvantage related to this methodology is the frequent occurrence of SI trends that are (in part) artefacts of the methodology, rather than accurate representations of physical phenomena, e.g., (i) extended periods with constant values of SI caused by consecutive differences in model fits on the non-selective and selective medium, e.g., during tailing; (ii) strange SI trends caused by differences in the model output behaviour on the non-selective and selective medium; and (iii) negative values of SI caused by intersecting model fits on the nonselective and selective media. Examples of these artefacts are illustrated in Figure 1. Given these artefacts, the methodology requires the assumption of zero SI when the model output on the selective medium is higher than on the non-selective medium. To avoid some of the aforementioned artefacts, Noriega et al. used a slightly different methodology; they directly used the log transformed model fit values on the non-selective and selective medium as

input for the SI equation [12]. On the one hand, the adapted approach resulted in more smooth SI evolutions characterised by less abruptly changing SI evolutions. On the other hand, the subpopulation of sublethally injured cells was not represented as a percentage of the total cell population, essentially not corresponding to the SI equation. Moreover, SI evolutions not explained by physical phenomena still occurred due to the SI evolution solely being based on the difference between model fits on non-selective and selective media. Therefore, Verheyen et al. used the raw cell count data directly as an input for the SI equation and fitted a third degree polynomial to the resulting SI data points [44,45]. Log transformed cell count data were used because fitting a model to the absolute values of an exponential process would lead to a large variance in SI values over the course of inactivation treatments. While this methodology was able to describe the SI-behaviour of microbial inactivation fairly well, some artefacts of using the third-degree polynomial as an SI model were observed, e.g., a sudden increase in SI at the end of inactivation treatments. In addition, the subpopulation of sublethally injured cells was not represented as a percentage of the total cell population due to the use of log transformed cell counts, similar to what was the case for the methodology of Noriega et al. [12]. Moreover, this method is a black box approach to describing the evolution of SI, ignoring any mechanistic knowledge that is available on the process. As such, there is also no guarantee that the selected polynomial is appropriate to represent the true underlying phenomena.

**Figure 1.** Examples of artefacts occurring due to calculating SI directly from non-log transformed model fit values. The modelled evolution of the population density on non-selective (—) and selective (- - -) medium is presented in the top row and the resulting SI in the bottom row. Both are presented as a function of unitless time. Left: long periods of constant SI caused by consecutive differences in model fits on the non-selective and selective medium. Middle: strange SI evolution caused by differences in model output shapes on non-selective and selective medium. Right: negative values in the SI evolution caused by intersecting model fits on the non-selective and selective media. For each case, top figures represent the model fits to non-selective and selective media, while bottom figures represent the SI calculated according to Equation (1).

The objective of this study was to develop a SI modelling method specifically tailored to use in combination with existing microbial inactivation models, to avoid disadvantages related to the approaches discussed above. In essence, all these disadvantages can be avoided if a more mechanistic modelling approach is used that can guarantee an appropriate relationship between the model fits to the non-selective and selective media. Therefore, a method was proposed to model the evolution of experimental data both on non-selective and selective media by a single model that describes the relationship between the corresponding populations of all culturable and healthy culturable cells. This methodology assumes that cells are first sublethally injured prior to inactivation; an assumption that can

be incorporated directly into primary inactivation models. The newly developed modelling concept was validated by means of three parameter estimation case studies from literature, comparing integrated SI modelling results to the results obtained with the traditional methodology that relies on independent models, as used in the respective papers.

#### **2. Materials and Methods**

#### *2.1. Datasets*

Three datasets from previous research by the authors' research group were used to evaluate the implementation of the novel sublethal injury (SI) model (presented in Section 3.1). For each dataset, the SI model was compared with the classical method of using independent models. The three datasets were selected specifically because they display different inactivation kinetics and therefore require different inactivation models as well. Dataset 1 was published by Smet et al. and displays simple log-linear inactivation kinetics [42]. Specifically, this dataset was obtained by applying cold atmospheric plasma treatments to *Salmonella* typhimurium that was grown as planktonic cells in an environment at pH 5.5 with 60 g/L NaCl. The reported inactivation curves were the result of the storage of samples at 8 ◦C after they had been treated with this plasma in a liquid carrier. Samples were collected over a period of 770 h and there were 22 samples for both the non-selective medium (TSA; Tryptic Soy Agar; Oxoid, Basingstoke, UK) and selective medium (XLD; Xylose Lysine Deoxycholate agar; Merck & Co, Rahway, NJ, USA). Dataset 2 and Dataset 3 have been published by Noriega et al. in the context of studying the effect of cell immobilization on SI during mild heat treatments (54 ◦C) [12]. Dataset 2 was from the inactivation of low-density colonies of *Listeria innocua* and displays a biphasic log-linear inactivation. There were 104 datapoints on non-selective medium and 87 on selective medium. The non-selective medium was TSA supplemented with 6 g/L yeast extract (Merck, Darmstadt, Germany) and the selective medium was the same supplemented with an additional 65 g/L NaCl (VWR, Leuven, Belgium). It should be noted that *L. innocua* was used as an innocuous surrogate organism for the food pathogen *L. monocytogenes*. Although this is a commonly used surrogate, it is not as tolerant to stress factors as *L. monocytogenes*. Dataset 3 was obtained from the inactivation of high-density colonies of *Salmonella* Typhimurium and displays log-linear inactivation with a tail. This dataset contained 130 datapoints on non-selective medium (TSA) and 89 datapoints on selective medium (TSA supplemented with 50 g/L NaCl). In this work, all data were converted to be expressed in the natural logarithm instead of the common logarithm of the cell density.

#### *2.2. Parameter Estimation and Uncertainty*

Model parameters were estimated using the function *lsqnonlin* of MATLAB R2019b (The MathWorks). The parameter estimation was based on the minimization of the sum of squared errors SSE:

$$\text{SSE} = \sum\_{\mathbf{i}=1}^{\mathbf{v}\_{\mathbf{m}}} \left( \mathbf{n}\_{\mathbf{m},\mathbf{i}}(\mathbf{t}\_{\mathbf{i}}) - \mathbf{n}\_{\mathbf{p},\mathbf{i}}(\mathbf{t}\_{\mathbf{i}},\mathbf{p}) \right)^{2} \tag{2}$$

with ν<sup>m</sup> the number of measurements in the dataset, nm,i the logarithm of the measured cell concentration, np,i the logarithm of the predicted cell concentration, ti the time point corresponding to a specific measurement and p the vector of model parameters. The different models will be described in the results and discussion section when used. The 95% confidence bounds of the estimated parameters were estimated using Equation (3) [46]:

$$\left[\mathbf{p}\_{\mathrm{i}} \pm \mathbf{t}\_{0.975, \ \mathrm{v}\_{\mathrm{m}} - \mathbf{v}\_{\mathrm{p}}} \cdot \mathbf{o}\_{\mathrm{P}\_{\mathrm{i}}}\right] \tag{3}$$

with t as the inverse Student's t-distribution, ν<sup>p</sup> the number of model parameters, and σpi the standard deviation of a model parameter pi . This standard deviation was calculated through Equations (4)–(7) [47]:

$$
\sigma\_{\rm Fe} = \sqrt{\rm V(i, i)}\tag{4}
$$

$$\mathbf{V} = \mathbf{F}^{-1} \tag{5}$$

$$\mathbf{F} = \frac{1}{\mathbf{M} \mathbf{S} \mathbf{E}} \mathbf{J}' \mathbf{J} \tag{6}$$

$$\text{MSE} = \frac{\text{SSE}}{\mathbf{v}\_{\text{m}} - \mathbf{v}\_{\text{P}}} \tag{7}$$

where V is the variance-covariance matrix of the model parameters, F is the Fisher Information Matrix, J is the Jacobian matrix and MSE is the mean sum of squared errors.

#### *2.3. Sublethal Injury*

Sublethal injury (SI) is defined here as the percentage of cells that cannot grow on selective media, but are able to grow on non-selective media:

$$\text{SI} = \frac{\text{N}\_{\text{non}-\text{selective}} - \text{N}\_{\text{selective}}}{\text{N}\_{\text{non}-\text{selective}}} \cdot 100\% \tag{8}$$

with Nnon−selective and Nselective the quantity of cells as detected on non-selective and selective media. It should be noted that the actual quantity of cells is used and not a logarithmic transformation of this quantity. As such, SI can be expressed as a percentage.

#### **3. Results and Discussion**

The current study presents a new modelling concept to describe sublethal injury (SI) as a function of time through sets of differential equations. A variety of inactivation models is described in literature. Together, these models have the ability to produce a wide variety of inactivation kinetics and they are selected based on the requirements of the specific scenario that is under study. Consequently, the goal of this work is not to present a specific model, but rather to present a method that allows the user to develop a new SI model based on any existing inactivation model that is available as a (set of) differential equation(s). The global concept of this modelling method is unveiled in Section 3.1. Then, this concept is applied to three different scenarios based on different datasets that required different inactivation models (Sections 3.2–3.4). For each scenario, the SI model is composed, the model parameters to be estimated are identified, the resulting model is tested on an available dataset and the results are compared with what would have been obtained if the data would be treated independently as done in the conventional method.

#### *3.1. The Sublethal Injury Model*

The main lack of the conventional approach to describe the evolution of SI is that the evolution of all culturable cells and of the healthy culturable cells as a function of time are treated as independent phenomena. Therefore, the main requirement for the new modelling concept is to link these inactivation kinetics. Composing this model concept is done by defining different subpopulations within the model. The distinct subpopulations are the healthy culturable cells NH, the sublethally injured (SI) culturable cells NSI and the dead (and non-culturable) cells ND. Additionally, the healthy cells and injured cells can be grouped together in the class of all culturable cells NC. This classification is presented in the Venn diagram of Figure 2. Based on these subpopulations, the inactivation of the microbial population can be described as the transition of cells from one subpopulation to another. A healthy cell can become a SI cell due to the inactivation method or survival conditions that it undergoes. It can be assumed that, under the studied conditions that cause inactivation of the cells, injured cells will be unable to become healthy cells again. Therefore, the rate of cells becoming injured is equal to the change in the quantity of healthy cells. This is described mathematically by a differential equation for dNH/dt. Experimentally, the evolution of the healthy population of cells is monitored through viable plate counting on selective media. Consequently, if injured cells are harmed further, their injuries become irreparable and they are thus converted to dead cells. If one assumes that the healthy cells are always first injured before they can be killed, the rate of injured

cells dying off is equal to the change in the population of all culturable cells. This is mathematically expressed by a differential equation for dNL/dt and the process can be monitored experimentally by using viable plate counts on non-selective medium. This assumption has been confirmed, e.g., in the research of Perni et al. on the inactivation of *E. coli* by pulsed electric fields [48]. If cells first have to become injured before being killed, this also means that the size of dNC/dt is proportional to the quantity of sublethally injured cells, not to the total quantity of culturable cells. The change of the population of injured cells is then a combination of the transition of healthy cells to injured cells and of injured cells to dead cells (dNSI/dt = −dNH/dt + dNC/dt). Within this modelling framework, two differential equations are needed to describe the evolution of the three subpopulations: (i) the change in the number of healthy cells dNH/dt and (ii) the change in the number of culturable cells dNC/dt. These two populations are monitored experimentally by plating respectively on selective and non-selective medium. As such, it is possible to identify the model equations on the data that are available out of the same types of experiments that are commonly being done to study SI. The data that are already being collected is sufficient, but the model itself needs to link the subpopulations as described above. By linking these subpopulations in a semi-empirical population balance model, there is a guarantee that the model will deliver reasonable predictions of the phenomena where, e.g., the number of healthy cells can never be higher than the total number of culturable cells. To make the proposed modelling concept more tangible, the next three subsections will present how to compose the model equations for the SI model based on three different inactivation models that are taken from literature.

**Figure 2.** Overview of the microbial subpopulations that are considered in the sublethal injury model. All cells can be divided in culturable and dead cells. The culturable cells can be subdivided in healthy and injured cells. The mathematical notations and conversions are marked in blue and the experimental quantifications are written in green.

#### *3.2. Case Study 1: Log-Linear Inactivation*

The first SI model that is composed here is based on simple log-linear inactivation for both the healthy and sublethally injured populations. As such, the model of Bigelow and Esty can be used to describe the evolution of the healthy population as a function of time [49]:

$$\frac{d\mathbf{N}\_{\rm H}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm H,max} \cdot \mathbf{N}\_{\rm H}(\mathbf{t})\tag{9}$$

with kH,max the maximum specific inactivation rate of the healthy population, which is expressed in an inverse time unit (e.g., 1/min, 1/s). For computational purposes, it is better to work with the logarithmic transformation of the population size:

$$\frac{d\mathbf{n}\_{\rm H}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm H,max} \tag{10}$$

with nH the logarithm of the size of the population of healthy cells. The next equation that can be composed is that for the total population of culturable cells:

$$\frac{d\mathbf{N}\_{\rm C}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm SL,max} \cdot \mathbf{N}\_{\rm SI}(\mathbf{t})\tag{11}$$

with kSI,max the maximum specific inactivation rate of the population of injured cells. This equation demonstrates that the decrease of the total quantity of culturable cells is only a function of the quantity of injured cells. This is a consequence of the assumption that cells first have to become injured before they can be killed. The quantity of injured cells can be calculated as the difference between the total population of culturable cells and the population of healthy cells (NSI(t) = NC(t) − NH(t)). As such, on a logarithmic scale, Equation (11) becomes:

$$\frac{d\mathbf{n}\_{\rm C}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm SL,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm H}(\mathbf{t}) - \mathbf{n}\_{\rm C}(\mathbf{t})) \right) \tag{12}$$

with nC the logarithm of the size of the population of culturable cells. The differential equation for the evolution of the population of injured cells itself is a combination of the increase due to the healthy cells becoming injured (Equation (9)) and a decrease due to the injured cells dying off (Equation (11)). The Equation (14) for the evolution of the injured population can therefore be composed as:

$$\frac{\text{dN}\_{\text{SI}}(\mathbf{t})}{\text{dt}} = -\frac{\text{dN}\_{\text{H}}(\mathbf{t})}{\text{dt}} + \frac{\text{dN}\_{\text{C}}(\mathbf{t})}{\text{dt}} \tag{13}$$

$$\frac{d\mathbf{N\_{SI}(t)}}{dt} = \mathbf{k\_{HI,max}} \cdot \mathbf{N\_{H}(t)} - \mathbf{k\_{SI,max}} \cdot \mathbf{N\_{SI}(t)}\tag{14}$$

Converting Equation (14) to the logarithmic scale results in the following equation:

$$\frac{d\mathbf{n}\_{\rm SI}(\mathbf{t})}{d\mathbf{t}} = \mathbf{k}\_{\rm HI,max} \cdot \frac{\exp(\mathbf{n}\_{\rm HI}(\mathbf{t}))}{\exp(\mathbf{n}\_{\rm SI}(\mathbf{t}))} - \mathbf{k}\_{\rm SI,max} \tag{15}$$

Given (i) that the quantity of injured cells at any point in time can be calculated from the difference between the amount of culturable and healthy cells and (ii) that the experimental data of viable plate counts on non-selective and selective media provides direct information about the quantity of culturable and healthy cells, it is not necessary to include Equation (15) on the evolution of the injured population when solving the differential equations that describe this SI model. Instead, the quantity of injured cells can be calculated from the results of the culturable and healthy populations. As such, the SI model for log-linear inactivation can be described by Equations (10) and (12). To solve this system of differential equations, the initial conditions are required. These initial conditions are:

$$\mathfrak{n}\_{\mathbb{H}}(\mathfrak{t}=0) = \mathfrak{n}\_{\mathbb{H},0} \tag{16}$$

$$
\mathfrak{n}\_{\mathbb{C}}(\mathfrak{t}=0) = \mathfrak{n}\_{\mathbb{C},0} \tag{17}
$$

with nH,0 and nC,0 the logarithm of the population sizes at the initial time point. These initial conditions cannot take any value as the initial quantity of all culturable cells should of course be larger than its subpopulation of initial healthy cells. This constraint can be added into the model by converting the first initial condition to:

$$\mathbf{n}\_{\rm H}(\mathbf{t}=\mathbf{0}) = \mathbf{f}\_{\rm H,0} \cdot \mathbf{n}\_{\rm C,0} \quad \text{with} \quad \mathbf{f}\_{\rm H,0} \in [0;1] \tag{18}$$

with fH,0 the fraction of the total culturable population of cells that are healthy cells. The full model can thus be summarised as:

$$\frac{d\mathbf{n}\_{\rm H}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm H,max} \quad \text{with} \quad \mathbf{n}\_{\rm H}(\mathbf{t}=0) = \mathbf{f}\_{\rm H,0} \cdot \mathbf{n}\_{\rm C,0} \quad \text{and} \quad \mathbf{f}\_{\rm H,0} \in [0;1] \tag{19}$$

$$\frac{d\mathbf{n}\_{\mathbb{C}}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\mathbb{S}\mathbf{l},\max} \cdot \left(1 - \exp(\mathbf{n}\_{\mathbb{H}}(\mathbf{t}) - \mathbf{n}\_{\mathbb{C}}(\mathbf{t})) \right) \quad \text{with} \quad \mathbf{n}\_{\mathbb{C}}(\mathbf{t} = \mathbf{0}) = \mathbf{n}\_{\mathbb{L},0} \tag{20}$$

with the following model parameters to be estimated: fH,0, nC,0, kH,max and kSI,max. In contrast, when considering the healthy and culturable populations as independent populations, they would be described by the following independent models:

$$\frac{d\mathbf{n}\_{\rm H}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm H,max} \quad \text{with} \quad \mathbf{n}\_{\rm H}(\mathbf{t}=0) = \mathbf{n}\_{\rm H,0} \tag{21}$$

$$\frac{d\mathbf{n}\_{\mathbb{C}}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\mathbb{C},\max} \quad \text{with} \quad \mathbf{n}\_{\mathbb{C}}(\mathbf{t}=0) = \mathbf{n}\_{\mathbb{C},0} \tag{22}$$

where kC,max is the maximum specific inactivation rate of the population of culturable cells. The model parameters in these independent models are nH,0, nC,0, kH,max and kC,max.

The SI model of Equations (19) and (20) is compared with the independent model of Equations (21) and (22) based on Dataset 1 (see Section 2.1). For this comparison, the model parameters of each set of equations were estimated and the results are presented in Table 1. The parameters nC,0 and kH,max, which appear in both models, have identical values and confidence bounds. The parameters kSI,max and kC,max of respectively the SI model and the independent models also have identical values and confidence bounds. These parameters are essentially each other's equivalent for the dataset that was used in this example. In contrast, the initial number of cells in the independent models is described by the parameter nH,0, whereas for the SI model the parameter fH,0 determines the fraction of the nC,0 that are healthy cells. Relative to their sizes, the uncertainty on the parameter fH,0 appears to be high compared to that of nH,0. However, the parameter fH,0 expresses the fraction of the population size on a linear scale whereas the parameter nH,0 follows a logarithmic scale. On the other hand, there is indeed additional uncertainty in the calculation of the initial population of healthy cells in the SI model because its value depends on both fH,0 and nC,0, which are both marked by uncertainty. The model predictions of the SI model and the independent models are visually indistinguishable and are therefore both represented by the graphs in Figure 3a. The sublethal injury (SI) is calculated as a function of time based on the following equation:

$$\text{SI}(\mathbf{t}) = \frac{\exp(\mathbf{n}\_{\text{C}}(\mathbf{t})) - \exp(\mathbf{n}\_{\text{H}}(\mathbf{t}))}{\exp(\mathbf{n}\_{\text{C}}(\mathbf{t}))} \tag{23}$$


**Table 1.** Parameter estimation results from fitting the independent models (Equations (21) and (22)) and sublethal injury model (Equations (19) and (20)) based on log-linear inactivation on Dataset 1.

**Figure 3.** (**a**) Log-linear inactivation of all culturable (—) and healthy culturable (-·-) cells as described by the independent models or sublethal injury model that is constructed with the model of Bigelow and Esty [49]. Based on these modelling results, the evolution of the data are sublethally injured population (···) and (**b**) the percentage sublethal injury (**—**) is calculated. The based on counts of culturable (x) and healthy (o) cells on non-selective and selective media as published in Smet et al. [42]. The dataset was constructed by monitoring the survival of *Salmonella* Typhimurium, which was grown as planktonic cells in an environment at pH 5.5 with 60 g/L NaCl, after being treated in a liquid carrier with cold atmospheric plasma and stored at 8 ◦C (Dataset 1).

The evolution of SI with time as calculated by both models is illustrated in Figure 3b.

This comparison of the model parameters and predictions of the SI model and independent models based on Dataset 1 demonstrates that the SI model is in a way equivalent to the independent models in cases where there can be no mistake in the calculation of the SI. Even though the equations and the model parameter constraints of both models are different, the model predictions are quasi-identical and the model parameters are equivalent. As will be demonstrated in the next section, these models start to differ when the independent model would yield wrongful predictions of the SI.

#### *3.3. Case Study 2: Biphasic Inactivation*

Dataset 2 portrays inactivation data that follows biphasic inactivation kinetics, as can be seen from Figure 4. Biphasic microbial inactivation is typically modelled with the model of Cerf [50]. Independent models for the healthy and culturable cells can be defined based on the model of Cerf as follows:

$$\frac{d\mathbf{n}\_{\rm HI}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm HI,max} \quad \text{with} \quad \mathbf{n}\_{\rm HI}(\mathbf{t}=0) = \mathbf{n}\_{\rm HI,0} \tag{24}$$

$$\frac{d\mathbf{n}\_{\rm H2}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm H2,max} \quad \text{with} \quad \mathbf{n}\_{\rm H2}(\mathbf{t}=0) = \mathbf{n}\_{\rm H2,0} \tag{25}$$

$$\frac{d\mathbf{n}\_{\rm C1}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm C1,max} \quad \text{with} \quad \mathbf{n}\_{\rm C1}(\mathbf{t}=0) = \mathbf{n}\_{\rm C1,0} \tag{26}$$

$$\frac{d\mathbf{n}\_{\rm{C2}}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm{C2},\rm{max}}\quad\text{with}\quad\mathbf{n}\_{\rm{C2}}(\mathbf{t}=0) = \mathbf{n}\_{\rm{C2},0} \tag{27}$$

with nH1 and nH2 the healthy cells and nC1 and nC2 the culturable cells of, respectively, the first and second subpopulation. If one subpopulation has both a lower initial size and a lower inactivation rate than the other, the typical biphasic inactivation is observed [51]. The parameters of the independent models that have to be estimated are nH1,0, nH2,0, nC1,0, nC2,0, kH1,max, kH2,max, kC1,max, and kC2,max. When constructing the SI model based on the model of Cerf, the only change in the equations for the evolution of healthy cells in each subpopulation compared to Equations (24) and (25) is that the initial number of

cells is defined as a fraction of the initial number of total culturable cells in the respective subpopulation (fH1,0 and fH2,0):

$$\frac{d\mathbf{n}\_{\rm HI}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm HI,max} \quad \text{with} \quad \mathbf{n}\_{\rm HI}(\mathbf{t}=0) = \mathbf{f}\_{\rm HI,0} \cdot \mathbf{n}\_{\rm Cl,0} \text{ and } \mathbf{f}\_{\rm HI,0} \in [0;1] \tag{28}$$

dnH2(t) dt <sup>=</sup> <sup>−</sup>kH2,max with nH2(<sup>t</sup> <sup>=</sup> <sup>0</sup>) <sup>=</sup> fH2,0 · nC2,0 and fH2,0 [0; 1] (29)

**Figure 4.** Biphasic inactivation of all culturable (—) and healthy culturable (-·-) cells as described by (**a**) the independent models and (**c**) the sublethal injury model that are constructed with the model of Cerf [50]. Based on these modelling results, the evolution of the sublethally injured population (···) and the percentage sublethal injury (—, (**b**,**d**) are calculated. The data are based on counts of culturable (x) and healthy (o) cells on non-selective and selective media as published by Noriega et al. [12]. The data were obtained by heat-treating low-density colonies of *Listeria innocua* at 54 ◦C (Dataset 2).

The differential equations for the two subpopulations of culturable cells are analogous to the differential equation for the culturable cells in the SI model for log-linear inactivation (Equation (20)):

$$\frac{d\mathbf{n\_{C1}(t)}}{dt} = -\mathbf{k\_{S1,max}} \cdot \left(1 - \exp(\mathbf{n\_{H1}(t)} - \mathbf{n\_{C1}(t)})\right) \quad \text{with} \quad \mathbf{n\_{C1}(t=0)} = \mathbf{n\_{C1,0}} \tag{30}$$

$$\frac{d\mathbf{n}\_{\rm C2}(\mathbf{t})}{d\mathbf{t}} = -\mathbf{k}\_{\rm S2,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm H2}(\mathbf{t}) - \mathbf{n}\_{\rm C2}(\mathbf{t})) \right) \quad \text{with} \quad \mathbf{n}\_{\rm C2}(\mathbf{t} = 0) = \mathbf{n}\_{\rm LC,0} \tag{31}$$

The model parameters of the SI model in Equations (28)–(31) are fH1,0, fH2,0, nC1,0, nC2,0, kH1,max, kH2,max, kC1,max and kC2,max. When comparing the model structure and parameterization between the two models, the differences in the SI model (Equations (28)–(31)) compared to the independent model (Equations (23)–(27)) are (i) that the initial condition of the healthy subpopulation is defined as a fraction of the initial culturable subpopulation and (ii) that the inactivation rate of the culturable cells depends on the number of injured cells instead of on the number of culturable cells.

Both models were fitted on the Dataset 2. The resulting model outputs and the calculated SI evolution are illustrated in Figure 4. In the results for the independent models, the model output for the culturable cells decreases below the model output for the healthy cells. In practice, it is of course not possible that the number of healthy cells would ever be higher than the number of culturable cells since the former is a subpopulation of the latter. However, due to this error, the number of injured cells that is calculated becomes negative and its logarithm a complex number. This is all because the independent models do not consider that the data are meant to describe the phenomenon of SI. The model output of the SI model on the other hand provides a good approximation of the experimental data and a reasonable calculation of the evolution of SI. The parameter estimation results for the independent models in Table 2 do not reveal the problems that occur when calculating the SI. Since the independent models fit the respective datasets for culturable and healthy cells well, the model parameter uncertainty is relatively low and the model fit is decent (as seen from the MSE). The problems in this approach only occurs when using the model outputs for calculating the SI. When comparing the parameters of the independent models and the SI model, the large difference in the initial populations catches the eye. These values of the initial populations are part of predicting the biphasic behaviour. In the SI model, the decrease of the culturable population is bound by the number of healthy cells. As such, there is no need for the biphasic behaviour of the inactivation model to predict the levelling off of the culturable cells towards the healthy cells. Consequently, when fitting the model parameters of the SI, the flexibility of the biphasic model is used to predict faster inactivation in the first 10 h compared to the time interval between 10 and 70 h (see Figure 4c). This phenomenon can also be observed from comparing the inactivation rate of the culturable population at time zero in the SI model (as determined by kSI2,max) with that of the independent models (as determined by kC2,max). This comparison shows that the initial inactivation rate in the SI model is much higher than in the independent models. These differences in the model structure also resulted in a better model fit of the SI model, as indicated from the MSE. Due to its constraints, the SI model will not always provide a better approximation of experimental data than the independent models. On the other hand, these constraints cause the SI model to be able to describe some phenomena with less model complexity than the independent models. As an example, in the current modelling results for the SI model, the tail in the culturable population of cells is a result from the tail in the healthy population of cells. As such, this tailing is described by the model without the need for it to be included separately in the differential equations of the culturable population of cells. In this case, the additional model complexity of the two phases is used to obtain a better approximation of the microbial inactivation in the initial stage of the experiment. The next section will further demonstrate the possibility to work with a decreased model complexity in the SI model.


**Table 2.** Parameter estimation results from fitting the independent models (Equations (24) to (27)) and sublethal injury model (Equations (28) to (31)) based on biphasic inactivation on Dataset 2.


**Table 2.** *Cont.*

This example demonstrates that the proposed SI model provides a guarantee for reasonable predictions on the evolution of SI by incorporating basic knowledge on the relationship of the different microbial subpopulation into the model structure.

#### *3.4. Case Study 3: Log-Linear Inactivation with Tailing*

Figure 5a illustrates the log-linear inactivation behaviour with tailing of both the population of culturable cells and the subpopulation of healthy cells in Dataset 3. Based on the inactivation model of Geeraerd et al., the following independent models are proposed to describe the populations of culturable and healthy cells [52]:

$$\frac{d\mathbf{n}\_{\rm HI}(t)}{dt} = \mathbf{k}\_{\rm H,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm H,res} - \mathbf{n}\_{\rm H}(t))\right) \quad \text{with} \quad \mathbf{n}\_{\rm H}(t=0) = \mathbf{n}\_{\rm H,0} \tag{32}$$

$$\frac{d\mathbf{n}\_{\rm C}(\mathbf{t})}{d\mathbf{t}} = \mathbf{k}\_{\rm C,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm C,res} - \mathbf{n}\_{\rm C}(\mathbf{t})) \right) \quad \text{with} \quad \mathbf{n}\_{\rm C}(\mathbf{t} = 0) = \mathbf{n}\_{\rm C,0} \tag{33}$$

with nH,res and nC,res the resistant subpopulations of the healthy and culturable cells. Converting these models into a linked SI model, results in the following set of equations:

$$\frac{d\mathbf{n}\_{\rm H(t)}}{dt} = \mathbf{k}\_{\rm H,max} \cdot \left(1 - \exp\left(\mathbf{n}\_{\rm H,res} - \mathbf{n}\_{\rm H}(\mathbf{t})\right)\right) \tag{34}$$

$$\frac{\text{dn}\_{\text{C}(\text{t})}}{\text{dt}} = \text{k}\_{\text{SI},\text{max}} \cdot \left(1 - \frac{\exp(\text{n}\_{\text{SI},\text{res}})}{\exp(\text{n}\_{\text{SI}}(\text{t}))}\right) \cdot \frac{\exp(\text{n}\_{\text{SI}}(\text{t}))}{\exp(\text{n}\_{\text{C}}(\text{t}))} \tag{35}$$

with nSI,res the resistant population of injured cells. To obtain a set of differential equations that is only a function of the healthy and culturable cells, Equation (35) is rewritten to:

$$\frac{d\mathbf{n}\_{\mathbb{C}(t)}}{dt} = \mathbf{k}\_{\mathrm{l,max}} \cdot \left(1 - \frac{\exp(\mathbf{n}\_{\mathbb{C},\text{res}}) - \exp(\mathbf{n}\_{\mathbb{H},\text{res}})}{\exp(\mathbf{n}\_{\mathbb{C}}(\mathbf{t})) - \exp(\mathbf{n}\_{\mathbb{H}}(\mathbf{t}))}\right) \cdot \frac{\exp(\mathbf{n}\_{\mathbb{C}}(\mathbf{t})) - \exp(\mathbf{n}\_{\mathbb{C}}(\mathbf{t}))}{\exp(\mathbf{n}\_{\mathbb{C}}(\mathbf{t}))} \tag{36}$$

If the model in Equation (36) is used, the resistant population of healthy cells needs to be described as a fraction of the resistant population of culturable cells fH,res:

$$\mathbf{n}\_{\rm H,res} = \mathbf{f}\_{\rm H,res} \cdot \mathbf{n}\_{\rm C,res} \tag{37}$$

When assuming all cells that can become injured can also be killed, there would be no residual population of injured cells. Consequently, the SI model for log-linear inactivation with tailing could be written as:

$$\frac{d\mathbf{n}\_{\rm H(t)}}{dt} = \mathbf{k}\_{\rm H,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm H,res} - \mathbf{n}\_{\rm H}(\mathbf{t})) \right) \quad \text{with} \quad \mathbf{n}\_{\rm H}(\mathbf{t} = 0) = \mathbf{n}\_{\rm H0} \tag{38}$$

$$\frac{d\mathbf{n}\_{\rm C(t)}}{dt} = \mathbf{k}\_{\rm SL,max} \cdot \left(1 - \exp(\mathbf{n}\_{\rm H}(\mathbf{t}) - \mathbf{n}\_{\rm C}(\mathbf{t})) \right) \quad \text{with} \quad \mathbf{n}\_{\rm C}(\mathbf{t} = 0) = \mathbf{n}\_{\rm C0} \tag{39}$$

**Figure 5.** Log-linear inactivation with tailing of all culturable (—) and healthy culturable (-·-) cells as described by (**a**) the independent models and (**c**) the sublethal injury model that are constructed with the model of Geeraerd et al. [52]. Based on these modelling results, the evolution of the sublethally injured population (···) and the percentage sublethal injury (—), (**b**,**d**) are calculated. The data are based on counts of culturable (x) and healthy (o) cells on non-selective and selective media as published by Noriega et al. [12]. The data were obtained by heat-treating high-density colonies of *Salmonella* Typhimurium at 54 ◦C (Dataset 3).

These two equations basically describe log-linear inactivation with tailing for the subpopulation of healthy culturable cells, while the decrease of the total culturable population depends on the log-linear inactivation of the injured subpopulation. In this model, the decrease of the population of culturable cells will asymptotically approach the resistant population of healthy cells. The independent models of Equations (32) and (33) and the SI model of Equations (38) and (39) are fitted to Dataset 3. The model outputs of both models and the calculated SI are presented in Figure 5. For the independent models, the resistant population of the culturable cells is slightly higher than that of the healthy cells, suggesting that a fraction of the injured cells would be resistant as well. On the other hand, according to the SI model, the total population of culturable cells approaches the resistant population of healthy cells under the assumption that all injured cells are killed. Based on the current data and modelling results, there is of course no way of telling which model predicts the correct inactivation of each subpopulation. However, it would be perfectly possible to use the extended SI model of Equations (34) and (36), which would give practically the same model output as the simplified model of Equations (38) and (39). Because the estimation of the resistant population often comes with quite some uncertainty, it often happens that the resistant population of the culturable cells drops below that of the healthy cells, resulting in a negative SI percentage (similar to Case study 2). This problem will in any case be avoided when using the SI model that is presented here. The model parameters that are estimated are listed in Table 3. For the main part, the parameters are again equivalent between

both models, with the main difference that the SI model now has one parameter less than the independent models. Even though the number of model parameters is reduced from six to five, the model accuracy only reduces slightly as shown by the MSE that changes from 1.59 to 1.63.

**Table 3.** Parameter estimation results from fitting the independent models (Equations (32) and (33)) and sublethal injury model (Equations (38) and (39)) based on log-linear inactivation with tailing on Dataset 3.


#### **4. Conclusions**

A new method was proposed to model the evolution of sublethal injury (SI) during microbial inactivation. This method relies on defining the two subpopulations of healthy and injured cells within the culturable cells and describing inactivation as a two-step mechanism: (i) injuring healthy cells and (ii) killing injured cells. Based on three case studies, it was illustrated that this model concept can be applied to any existing dynamic inactivation model that is available in literature. Moreover, these case studies demonstrated, respectively, three important properties of the SI model: (i) the model is equivalent to the conventional method of using independent models when there is sufficient difference between the total quantity of culturable cells and the subpopulation of healthy cells, i.e., when there can be no mistake in the calculated SI; (ii) the new modelling method avoids unrealistic calculations of the SI by using its mechanistic definition of this process; and (iii) in some cases, it is possible to reduce the number of model parameters in the SI model because there is no need for additional model parameters to describe phenomena that are already hard coded into the more mechanistic description of the model. These properties make the SI model suitable for describing SI as a function of time during food processes. As such, this model can contribute to the optimization and control of the safety of food processes by considering SI dynamics.

**Author Contributions:** Conceptualization, S.A. and D.V.; methodology, S.A.; validation, S.A.; investigation, S.A.; data curation, D.V. and C.S.; writing—original draft preparation, S.A., D.V., C.S. and J.F.M.V.I.; writing—review and editing, S.A., D.V., C.S. and J.F.M.V.I.; visualization, S.A. and C.S.; supervision, J.F.M.V.I.; project administration, S.A. and J.F.M.V.I.; funding acquisition, S.A., D.V. and J.F.M.V.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** Authors Simen Akkermans and Davy Verheyen were supported by the Research Foundation Flanders (FWO, grants 1224620N and 1254421N, respectively). This work was supported by the KU Leuven Research Fund through project C24/18/046, by the Research Foundation - Flanders through project G0B4121N and by the EU H2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 956126.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data is available without restrictions upon request through contact with the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


52. Geeraerd, A.H.; Herremans, C.H.; Van Impe, J.F. Structural model requirements to describe microbial inactivation during a mild heat treatment. *Int. J. Food Microbiol.* **2000**, *59*, 185–209.
