*Article* **Mechanistic Modelling of Biomass Growth, Glucose Consumption and Ethanol Production by** *Kluyveromyces marxianus* **in Batch Fermentation**

**Yolocuauhtli Salazar <sup>1</sup> , Paul A. Valle 2,\* , Emmanuel Rodríguez <sup>2</sup> , Nicolás O. Soto-Cruz <sup>3</sup> , Jesús B. Páez-Lerma <sup>3</sup> and Francisco J. Reyes-Sánchez <sup>3</sup>**


**Abstract:** This paper presents results concerning mechanistic modeling to describe the dynamics and interactions between biomass growth, glucose consumption and ethanol production in batch culture fermentation by *Kluyveromyces marxianus* (*K. marxianus*). The mathematical model was formulated based on the biological assumptions underlying each variable and is given by a set of three coupled nonlinear first-order Ordinary Differential Equations. The model has ten parameters, and their values were fitted from the experimental data of 17 *K. marxianus* strains by means of a computational algorithm design in Matlab. The latter allowed us to determine that seven of these parameters share the same value among all the strains, while three parameters concerning biomass maximum growth rate, and ethanol production due to biomass and glucose had specific values for each strain. These values are presented with their corresponding standard error and 95% confidence interval. The goodness of fit of our system was evaluated both qualitatively by in silico experimentation and quantitative by means of the coefficient of determination and the Akaike Information Criterion. Results regarding the fitting capabilities were compared with the classic model given by the logistic, Pirt, and Luedeking–Piret Equations. Further, nonlinear theories were applied to investigate local and global dynamics of the system, the Localization of Compact Invariant Sets Method was applied to determine the so-called localizing domain, i.e., lower and upper bounds for each variable; whilst Lyapunov's stability theories allowed to establish sufficient conditions to ensure asymptotic stability in the nonnegative octant, i.e., **R** 3 <sup>+</sup>,0. Finally, the predictive ability of our mechanistic model was explored through several numerical simulations with expected results according to microbiology literature on batch fermentation.

**Keywords:** asymptotic stability; batch fermentation; in silico experimentation; *Kluyveromyces marxianus*; nonlinear data fitting; nonlinear mechanistic model

### **1. Introduction**

Alcoholic fermentation is an anaerobic process that transforms sugars like glucose or fructose into ethanol and carbon dioxide. Several yeast species are used commonly in this process, e.g., *Kloeckera, Hanseniaspora, Candida, Pichia, Kluyveromyces*, and *Saccharomyces* among others. The growth rate of these microorganisms has an ultimate effect on the sensorial characteristics of the final product, which can be positive or negative depending on the yeast used [1].

Overall, yeasts are indispensable for biotechnological processes such as wine and beer production [2]. In this research, we focus on investigating glucose consumption and

**Citation:** Salazar, Y.; Valle, P.A.; Rodríguez, E.; Soto-Cruz, N.O.; Páez-Lerma, J.B.; Reyes-Sánchez, F.J. Mechanistic Modelling of Biomass Growth, Glucose Consumption and Ethanol Production by *Kluyveromyces marxianus* in Batch Fermentation. *Entropy* **2023**, *25*, 497. https:// doi.org/10.3390/e25030497

Academic Editor: Pavel Kraikivski

Received: 12 January 2023 Revised: 8 March 2023 Accepted: 10 March 2023 Published: 14 March 2023

**Copyright:** © 2023 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/).

ethanol production from several strains of *Kluyveromyces marxianus* (*K. marxianus*). This yeast has a great potential for alcoholic fermentation due to its intraspecific characteristics such as higher specific growth rates, the ability to grow on various substrates, and tolerance to high temperatures [3–5]. Further, *Kluyveromyces* sp. produces aromatic compounds such as fruity esters, carboxylic acids, ketones, furans, and other alcohols in liquid fermentation such as 2-phenyl ethanol whose sensorial characteristics can influence the quality of wine, distilled drinks, and fermented foods [6]; refer to Fonseca et al. for an extensive review on the biotechnological potentials of *K. marxianus* [4,6].

Concerning industrial applications, fermentation is commonly performed in batch culture, which brings certain advantages such as the reduction of contamination risk, in addition to the fact that a large capital investment is not necessary since high-priced production equipment is not required compared to a continuous culture process [7]. Batch process implies that yeasts are incubated in a closed container under controlled conditions with a culture medium composed of the necessary nutrients [8]. Hence, biomass cannot grow indefinitely and four phases have been identified in its dynamics, i.e., lag phase, exponential growth, stationary state, and death phase. While this process is carried out, the substrate is consumed and converted into the product, e.g., ethanol produced by sugars such as glucose [9]. Therefore, properly identifying the time interval of these phases as well as predicting the maximum product concentration that could be produced from the initial concentrations of both substrate and biomass may help to optimize production costs on the resulting product of several applications. The latter may be achieved by mechanistic modeling through predictive microbiology, which can be considered a powerful tool to investigate and summarize the overall effects of varying conditions and environmental factors within food formulation and processing [10]. Further, mathematical models could aid in gaining insights concerning microbial food safety and quality assurance of increasingly complex food products [11,12], as well as estimating shelf life and forecasting food spoilage [13,14].

Mathematical models in predictive microbiology can be classified according to different criteria, uses, and functionalities that are not mutually exclusive. Based on the type and number of variables, models are classified into primary, secondary, and tertiary; they can also be differentiated on the basis of their mathematical background as mechanistic or empirical [15], and they can be categorized into structured and unstructured conforming to the complexity of the chemical compounds of the biomass [16]. Primary models are those that represent biomass growth dynamics as a function of time, the main equations in the literature are the exponential functions of Gompertz [17] and Vazquez-Murado [18], the logarithmic function of Baranyi et al. [19] and the cubic model of Garcia et al. [20]. All models are described by parameters such as maximum growth rate [*µ*max], lag time [*L*], and both initial [*X*0] and maximum biomass [*X*max] Concentrations, while secondary models relate to the latter with environmental conditions such as temperature and pH, and other variables such as substrate and product concentrations over time, e.g., equations of Monod [21], Teisser [22], Haldane [23], and Moser [24], which aim to describe biomass growth dynamics as a function of substrate concentration and have been widely used to investigate bacterial growth [25]. Tertiary models are the result of combining primary and secondary models through the use of computer tools that allow predictions regarding the growth or death of microorganisms in food when different environmental conditions are combined [26]. Concerning the second classification mentioned, mechanistic models are formulated by means of theoretical bases and provide an interpretation of microbial growth in terms of known processes and empirical models are usually composed of polynomials of the first or second degree and pragmatically describe the data with convenient mathematical relationships, this does not usually give information on precise responses of microorganisms, because they do not take into account known processes [27]. Finally, according to the third category described, unstructured models consider biomass only as a chemical compound in a culture and its dynamics is described by simple models, while structured models also take into account changes in the internal cellular structure of

biomass in terms such as the content of RNA, enzymes, reagents, metabolism and products [28]. The Gompertz, Vazquez-Murado, Baranyi and Garcia models, mentioned above, are also classified as unstructured models since biomass is considered a variable described only by its concentration. Mathematical models used by Sansonetti [29], Lei [30], and Steinmeyer [31] are classified as structured because they describe the growth of biomass considering the intracellular reactions produced by its metabolism.

Thus, it is important to highlight that in a batch fermentation process, multiple reactions occur, so the adaptability and evolution of microorganisms in short periods and changes in environmental conditions usually characterize this type of process, consequently, the modeling of these systems is complex due to time-varying characteristics of biological systems, resulting in nonlinear systems dynamics [28]. Hence, a mathematical model formulated from a system of nonlinear differential equations will allow the application of nonlinear systems control methods to optimize the process so that the characteristics of the final product can be predicted when the environmental conditions of the culture are controlled and the initial conditions of biomass, substrate, and product values are known. It is worth mentioning that most of the models found in the literature focus on the yeast *Saccharomyces cerevisiae* since it is one of the most used in the industry; however, biotechnological opportunities have been found in non-Saccharomyces yeasts because they have metabolic characteristics that lead to the production of compounds of interest. Therefore, it is important to model the growth of *K. marxianus* because of the great potential in the production of esters compounds of industrial importance [32]. Thus, in this paper, we applied mechanistic and computational modeling to formulate a system of three coupled nonlinear first-order Ordinary Differential Equations (ODEs) that describe dynamics between biomass, glucose (substrate), and ethanol (product) concentrations over time. Mechanistic modeling allowed us to provide both qualitative and quantitative descriptions concerning the relationships of biomass growth, glucose consumption, and ethanol production from 17 strains of *K. marxianus*, while computational modeling was used to fit experimental data from these three variables and establish numerical values for each parameter of the mathematical model. Further, nonlinear theories such as the Localization of Compact Invariant Sets (LCIS) method and Lyapunov's Stability Theory were applied to provide a complete analysis of the local and global dynamics of our proposed biological system [33].

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

This section provides all the information concerning the experimental data of biomass growth, substrate consumption and ethanol production, i.e., karyotypes of the *K. marxianus* strains with identifiable chromosomal differences among them, environmental conditions, chemical characteristics of the medium, lab equipment used for measurements, and periods for each measurement, then the mathematical model is formulated and each equation as well as values and units of parameters are described. This section concludes by describing the procedure to approximate the experimental data and to fit the numerical values of each parameter by designing an algorithm in Matlab.

### *2.1. Experimental Data: Culture Medium and Analytical Techniques*

Experimental data was obtained from alcoholic fermentation in batch culture by *K. marxianus*, 17 strains with different genetic profiles were incubated in 20 g/L of yeast extract peptone dextrose agar at 30 ◦C in order to study their kinetic growth, glucose consumption and ethanol production. Codification and origin of studied karyotypes of *K. marxianus* are identified by Páez et al. in [34], where 15 strains were obtained from different places of México, and they were isolated from agave fermentation for mezcal production, in addition to 2 reference strains that were isolated from pozol (CBS6556) in México, and from yoghurt (CBS397) in Netherlands.

Characteristics of the chemically defined medium are given as follows: glucose 20 g/L, KH2PO<sup>4</sup> 3 g/L, (NH4)2SO<sup>4</sup> 3 g/L, Na2HPO<sup>4</sup> 1.49 g/L, glutamic acid 1 g/L, MgCl<sup>2</sup> heptahydrate 0.41 g/L, ZnCl<sup>2</sup> 0.0192 g/L, CuCl<sup>2</sup> 0.0006 g/L, MnCl<sup>2</sup> 0.044 g/L, CoCl<sup>2</sup> 0.0005 g/L, CaCl 0.0117 g/L, FeCl<sup>2</sup> 0.011 g/L, (NH4)6Mo7O<sup>24</sup> 0.004 g/L, H3BO<sup>4</sup> 0.0030 g/L, aminobenzoic acid 0.0010 g/L, myo-inositol 0.1250 g/L, nicotinic acid 0.0050 g/L, pantothenic acid 0.005 g/L, pyridoxine 0.0050 g/L, thiamin HCl 0.005 g/L, biotin 0.000024 g/L [35]. This medium was used to culture the strains for biomass development with agitation, for the conservation of the strains, plates with the same medium with 20% agar were used and stored at 4 ◦C.

Biomass concentration was measured with a spectrophotometer UV-VIS DR 6000 (HACH, Loveland, CO, USA) by optical density at 600 nm, values in g/L were obtained relating optical density with a calibration curve of the dry weight of *K. marxianus*. For glucose consumption and ethanol production by High-Performance Liquid Chromatography (HPLC series 1200, Agilent Technology, Palo Alto, CA, USA), a BIORAD HP-87*H*+(8%) ion exchange column was used, in an AGILENT® 1200 series equipment, with H2SO<sup>4</sup> 0.005 N as mobile phase, at a flow of 0.5 mL/min, the column temperature was 60 ◦C, and the Refractive Index detector temperature was 60 ◦C. The injection volume of 5 µL, calibration curves were made with glucose and ethanol Sigma Aldrich at 99% purity or higher, and a determination coefficient higher than 0.99 for each compound [36,37].

Fermentation was made in duplicate for every strain and samples were taken each hour for 13 consecutive hours. Two samples were taken every hour for each variable in the time interval of the process where *t* goes from 0 to 13, then the average value of the two measurements was computed. Therefore, each variable, i.e., biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)], has 14 observations with a total of 42 experimental data points (*n*) for each *K. marxianus* strain.

### *2.2. The KM Mechanistic Model*

The KM mechanistic model is proposed to describe the dynamics of alcoholic fermentation. This is a biochemical process carried out by yeasts (also known as biomass), to transform sugars such as glucose into ethyl alcohol, otherwise known as ethanol (main product) and other byproducts. In this case, the alcoholic fermentation is taken in a batch fermentation process with established laboratory conditions of temperature and known initial glucose concentrations (substrate). Our mathematical model describes the relationships between biomass concentration [*x*(*t*)], glucose consumption [*y*(*t*)], and ethanol production [*z*(*t*)] over time. The set of three first-order ODEs is presented below

$$
\dot{\mathfrak{x}}\_{\text{ }} = \frac{\rho\_1 \mathfrak{x} y}{\rho\_2 + y} - \rho\_3 \mathfrak{x} \underline{z} - \rho\_4 \mathfrak{x}\_{\text{ }} \tag{1}
$$

$$
\dot{y}\_{\quad} = -\rho\_5 xy - \rho\_6 yz - \rho \gamma y\_{\quad} \tag{2}
$$

$$
\dot{z}\_{\quad} = \rho\_8 xz + \rho\_9 yz - \rho\_{10} z\_{\quad} \tag{3}
$$

where each state variable *x*(*t*), *y*(*t*) and *z*(*t*) are measured in g/L, and the time unit is given in *hours*. Now, by considering results from Leenheer and Aeyels (see Section II.A in [38]), all solutions with nonnegative initial conditions [*x*(0), *y*(0), *z*(0) ≥ 0, ] will be located in the nonnegative octant as indicated below

$$\mathbf{R}\_{+,0}^3 = \{ \mathfrak{x}(t), y(t), z(t) \ge 0 \},$$

i.e., each positive half trajectory of the system will be positively forward invariant in **R** 3 <sup>+</sup>,0. The latter also considers the biological sense of each variable as there is no meaning for negative values of biomass, glucose or ethanol in the scope of the KM system (1)–(3). It is important to mention that variables cannot grow exponentially indefinitely, and they must have biologically feasible limits which will be discussed in the next section. Values and units of each parameter of the KM system (1)–(3) are given in Table 1.


**Table 1.** Description, values, and units of variables and parameters for the KM mechanistic model.

Now, let us describe our mechanistic model based on the experimental data described in the previous section and the following biological assumptions. Biomass growth dynamics is described by Equation (1), where the first term uses the classical Monod form for the growth of microorganisms [21], where *ρ*<sup>1</sup> is the biomass maximum growth rate (also found in the literature as *µ*max), and *ρ*<sup>2</sup> is the affinity or half-velocity constant for glucose consumption. The second term describes biomass death due to ethanol accumulation toxicity by the law of mass action (see Section 2.3 in [39]) with a rate *ρ*3. This term is negative because ethanol accumulation increases the membrane fluidity and negatively affects the membrane protein's function, which can lead to cell growth inhibition or even death [40,41]. The third term represents the natural yeast death rate [*ρ*4] mainly due to environmental resources depletion [42]. Glucose dynamics is formulated in Equation (2) as a decrescent function where the law of mass action gives the first two terms. The first one represents glucose consumption to support biomass growth. In contrast, the second term accounts for the glucose consumption used for ethanol production, with rates *ρ*<sup>5</sup> and *ρ*6, respectively. The third term represents the spontaneous decomposition rate of glucose [*ρ*7] [43]. The latter stems from the fact that the culture medium is placed in a sealed container in batch fermentation, and no other nutrients (primarily glucose) are supplied into the system. Ethanol dynamics is described in Equation (3). The first term represents ethanol production associated with biomass growth. It is due to ethanol being recognized as a primary metabolite, a product obtained from reactions required or cellular growth [44,45]. The second term represents the glucose conversion to ethanol not directly linked to cellular growth, attributed to the need for Gibbs's free energy for cellular maintenance [44,46]. In both cases, terms are formulated by the law of mass action with respective rates *ρ*<sup>8</sup> and *ρ*9. Finally, the third term represents ethanol degradation with a rate *ρ*10. The flow diagram shown in Figure 1 was constructed to illustrate the dynamics of the system.

**Figure 1.** Flow diagram describing interactions between each variable and their corresponding relationship with each parameter.

It should be noted that fixed parameters values were estimated for the 17 *K. marxianus* strains, particularly for *ρ*2, *ρ*3, *ρ*4, *ρ*5, *ρ*6, *ρ*<sup>7</sup> and *ρ*10. Further, concerning the death rate of biomass [*ρ*4], spontaneous decomposition rate of glucose [*ρ*7], and degradation rate of ethanol [*ρ*10], one can see that they are in a different order of magnitude and the following constraint is formulated for these three parameters:

$$
\rho\_{10} > \rho\_4 > \rho\_7. \tag{4}
$$

Now, concerning the equilibrium points of the KM system, Equations (1)–(3) have a unique biologically meaningful equilibrium point in the domain **R** 3 <sup>+</sup>,0 given by

$$(\mathfrak{x}\_{0\prime}^\* y\_{0\prime}^\* z\_0^\*) = (0,0,0). \tag{5}$$

Another set of five equilibria with at least one negative value is shown in Appendix A; therefore, these equilibrium points are discarded from the biologically meaningful dynamics of the system. Further, from the biological characteristics of each variable the following can be stated with respect to each solution as time increases

$$\lim\_{t \to \infty} x(t) = \lim\_{t \to \infty} y(t) = \lim\_{t \to \infty} z(t) = 0\_\prime$$

due to the eventual death of microorganisms, glucose consumption and ethanol degradation [47], asymptotic stability of trajectories is discussed in the next section.

### *2.3. Parameter Value Estimation*

First, let us compute the glucose decomposition rate [*ρ*7] by assuming a first-order kinetics [48] for glucose dynamics, and considering a half-life [*t*1/2] of 96 years [43]. Then, *ρ*<sup>7</sup> can be computed from the next equation

$$y(t) = y\_0 e^{-\rho\_7 t} \rho$$

as follows

$$\frac{y\_0}{2} = y\_0 e^{-\rho \gamma t\_{1/2}} \nu$$

where *y*<sup>0</sup> is the glucose initial concentration, i.e., *y*(0); hence

$$\rho\_7 = \frac{\ln 2}{t\_{1/2}} = 824.233 \times 10^{-9} \, h^{-1}.$$

Now, in order to determine the numerical values of parameters *ρ<sup>i</sup>* , *i* = 1, 2, 3, 4, 5, 6, 8, 9, 10; the computational model of Equations (1)–(3) was formulated as follows

$$\mathbf{x}\_{i+1} = \mathbf{x}\_{i} + \left(\frac{\rho\_{1}\mathbf{x}\_{i}y\_{i}}{\rho\_{2} + y\_{i}} - \rho\_{3}\mathbf{x}\_{i}z\_{i} - \rho\_{4}\mathbf{x}\_{i}\right)\Delta\_{t\prime} \tag{6}$$

$$y\_{i+1} = \\_y\_i + (-\rho\_5 x\_i y\_i - \rho\_6 y\_i z\_i - \rho\_7 y\_i)\Delta\_{t\prime} \tag{7}$$

$$z\_{i+1} = \;z\_i + (\rho\_8 \mathbf{x}\_i \mathbf{z}\_i + \rho\_9 y\_i \mathbf{z}\_i - \rho\_{10} \mathbf{z}\_i)\Delta\_{tr} \tag{8}$$

by applying Euler's method (see Section 1.7 in [49]) where <sup>∆</sup>*<sup>t</sup>* was set to <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> . Then, an algorithm was formulated in Matlab 2022b with the *lsqcurvefit* function from the optimization toolbox as its core [50] (initial points were set as <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>1</sup> for each parameter [*ρ*1, *ρ*8, *ρ*9], i.e., *<sup>x</sup>*<sup>0</sup> = [<sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>1</sup> ; 1 <sup>×</sup> <sup>10</sup>−<sup>1</sup> ; 1 <sup>×</sup> <sup>10</sup>−<sup>1</sup> ], and optimotions of the function were set as follows: Max Function Evaluations <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> , Max Iterations <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> , and Function Tolerance <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>9</sup> ). This allowed us to establish a fixed value for parameters *ρ<sup>j</sup>* , *j* = 2, 3, 4, 5, 6, 10; by averaging the corresponding values for each of the 17 strains, these results are shown in Table 1. However, this procedure was not applicable for parameters *ρ*1, *ρ*<sup>8</sup> and *ρ*9, as it was expected that each strain of *K. marxianus* will have its own biomass growth rate [*ρ*1], and its corresponding ethanol production rates (*ρ*<sup>8</sup> and *ρ*9), this is directly linked to the chromosomal differences among the strains affecting their growth kinetics. Hence, the main algorithm was redesigned to fit these three parameters and consider the others as fixed constants. Overall results are shown in Table 2 with their corresponding standard error (SE), and 95% confidence interval (CI). These two statistics allow us to establish that estimates for parameters *ρ*1, *ρ*8, and *ρ*<sup>9</sup> in the 17 strains are statistically significant. The latter follows from the fact that each *SE*(*ρ<sup>k</sup>* ) < *ρk*/2, *k* = 1, 8, 9; i.e., the value of the SE is less than half of the value fitted for each parameter, thus, the null hypothesis [*ρ<sup>k</sup>* = 0] can be rejected (see Section 5.2.8 from Koutsoyiannis [51]). Furthermore, both the lower and upper limit of the 95% CI of all fitted values are positive, hence, as there is no change in the sign of the bounds, this implies that the value of the null hypothesis is excluded, and one can conclude that all P-values are less than 0.05 (see Chapter 17 from Motulsky [52]).

**Table 2.** Fitted values, their standard error [*SE*(*ρ<sup>k</sup>* )], and 95% confidence intervals [*CI*(*ρ<sup>k</sup>* )] for the biomass growth rate [*ρ*1], and ethanol production rates [*ρ*8, *ρ*9], where all values are written with a magnitude of 10−<sup>3</sup> . Thus, it is possible to identify both lower and upper bounds for the values of the three fitted parameters as follows *<sup>ρ</sup>*<sup>1</sup> <sup>∈</sup> [289.385, 381.419] <sup>×</sup> <sup>10</sup>−<sup>3</sup> , *<sup>ρ</sup>*<sup>8</sup> <sup>∈</sup> [19.088, 49.816] <sup>×</sup> <sup>10</sup>−<sup>3</sup> , and *<sup>ρ</sup>*<sup>9</sup> <sup>∈</sup> [46.352, 70.349] <sup>×</sup> <sup>10</sup>−<sup>3</sup> .



**Table 2.** *Cont.*

Finally, it should be noted that the in silico experimentation performed in this research was done on a high-end desktop computer with a Ryzen 9 5950X CPU, 128 GB of RAM DDR4 CL18, a 12 GB GPU NVIDIA GeForce RTX 3080, and 1 TB Samsung 980 Pro Gen 4 NVMe M.2. The complete algorithm that was designed to fit the numerical values of parameters [*ρ*1, *ρ*8, *ρ*9], and to determine results concerning the statistics and goodness of fit can be found in the Supplementary Materials.

### **3. Results**

In this section, the in silico experimentation is performed by means of several numerical simulations, and results relating to the nonlinear analysis of the system are derived, i.e., bounds for the localizing domain, asymptotic stability, and existence and uniqueness for all solutions of our model in the nonnegative octant **R** 3 <sup>+</sup>,0.

### *3.1. In Silico Experimentation and Goodness of Fit*

First, qualitative results are illustrated by means of numerical simulations. For the sake of simplicity, the strains were clustered in groups of four from strain 1 to the 16 (see Figures 2–5, respectively), and results concerning only for the strain 17 are shown in Figure 6. In all panels, the × green marker represents the average value for the two experimental data measurements for each variable, i.e., biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)], whilst the blue continuous line represents the approximated value given by the KM system (1)–(3) when is solved by means of Equations (6)–(8) with <sup>∆</sup>*<sup>t</sup>* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> . The time units are given in *hours* and the concentration for each variable is measured in g/L as indicated in each axis. Values for all ten parameters corresponding to each strain are shown in Tables 1 and 2.

Now, let us provide a quantitative measure of the fitting capabilities of the KM mechanistic model (1)–(3), thus, the coefficient of determination - *R* 2 is calculated for each variable with results shown in Table 3.

] L=g[ )t( x

**Figure 2.** Each column from left to right (in landscape orientation) illustrates both the experimental data (× green marker), and the approximated values obtained with the KM system (continuous blue line) for each corresponding strain 1–4; the top row shows results for biomass [*x*(*t*)], the middle row for glucose [*y*(*t*)], and the lower row for ethanol [*z*(*t*)]. The × green marker represents the average value calculated from the two measurements that were made for each variable in every strain.

] L=g[ )t( z

] L=g[ )t( y

 x

**Figure 3.** Each column from left to right (in landscape orientation) illustrates both the experimental data (× green marker), and the approximated values obtained with the KM system (continuous blue line) for each corresponding strain 5–8; the top row shows results for biomass [*x*(*t*)], the middle row for glucose [*y*(*t*)], and the lower row for ethanol [*z*(*t*)]. The × green marker represents the average value calculated from the two measurements that were made for each variable in every strain.

] L=g[ )t( x

**Figure 4.** Each column from left to right (in landscape orientation) illustrates both the experimental data (× green marker), and the approximated values obtained with the KM system (continuous blue line) for each corresponding strain 9–12; the top row shows results for biomass [*x*(*t*)], the middle row for glucose [*y*(*t*)], and the lower row for ethanol [*z*(*t*)]. The × green marker represents the average value calculated from the two measurements that were made for each variable in every strain.

 y

 z

] L=g[ )t( x

**Figure 5.** Each column from left to right (in landscape orientation) illustrates both the experimental data (× green marker), and the approximated values obtained with the KM system (continuous blue line) for each corresponding strain 13–16; the top row shows results for biomass [*x*(*t*)], the middle row for glucose [*y*(*t*)], and the lower row for ethanol [*z*(*t*)]. The × green marker represents the average value calculated from the two measurements that were made for each variable in every strain.

 y

 z

**Figure 6.** Experimental data (×green marker), and approximated values obtained with the KM system (continuous blue line) for strain 17; the top panel shows results for biomass [*x*(*t*)], the panel row for glucose [*y*(*t*)], and the lower panel for ethanol [*z*(*t*)]. The × green marker represents the average value calculated from the two measurements that were made for each variable in every strain.

**Table 3.** The *R* <sup>2</sup> provides a measure of how well the experimental data are replicated by the KM mathematical model (1)–(3) for each strain. This coefficient was computed independently for each variable, i.e., biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)]. One can see that the values for *R* 2 ranges between 0.902 to 0.997 which allows us to conclude an overall well goodness of fit for the model.



**Table 3.** *Cont.*

Further, the Akaike Information Criterion (*AIC*) [53–55] was computed by considering a small sample relative to the number of parameters (*n*/*K* < 40) with a bias correction as indicated below

$$AIC = n \ln \left( \frac{\sum\_{i=1}^{n} [f\_e(i) - f\_a(i)]^2}{n} \right) + 2K + \frac{2K(K+1)}{n - K - 1} \lambda$$

where *n* is the total number of experimental data points; *f<sup>e</sup>* the experimental data and *f<sup>a</sup>* the approximated value for the residual sum of squares (*RSS*); and *K* the number of parameters of the system; therefore, *n*/*K* = (3 × 14)/10 = 4.2. Results, including *RSS*, *AIC* and *R* 2 , for the complete trajectory of the system, i.e., *φ*(*x*, *y*, *z*) for the total of 42 experimental points (14 for each variable) are summarized in Table 4.

**Table 4.** In order to provide overall measures for the fitting capabilities of our mathematical model, i.e., the KM system (1)–(3), values were calculated for the *RSS*, the *AIC*, and the *R* 2 to estimate and describe the dynamics between the three variables *φ*(*x*, *y*, *z*), where *x*(*t*), *y*(*t*) and *z*(*t*) represent, respectively, the evolution of biomass, glucose, and ethanol.



**Table 4.** *Cont.*

The *AIC* yields a value that relates the amount of information that our model loses when approximating the experimental data. Hence, one can compare the capabilities of the model to estimate the concentrations over time of biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)] among the 17 *K. marxianus* strains while providing a statistical measured for the quality of the KM system (1)–(3).

### *3.2. Nonlinear Analysis: Localizing Domain, Asymptotic Stability, Existence and Uniqueness*

The localizing domain can be determined by computing the upper bounds for all variables of the KM mechanistic model (1)–(3), the lower bounds are given by the boundary of the domain **R** 3 <sup>+</sup>,0, i.e., {*x*inf = 0, *y*inf = 0, *z*inf = 0}. The latter is achieved by means of integration and the LCIS method [56]. Within the localizing domain, one may find all biologically meaningful dynamics of the system, i.e., compact invariant sets such as equilibrium points, periodic orbits, limit cycles and chaotic attractors (see Section 3 in [57]), among others.

First, in order to find the upper bound for the glucose concentration [*y*(*t*)], Equation (2) is integrated as follows

$$\int \frac{dy}{y} = -\int\_0^t f(x, z)dt.$$

where

$$f(\mathbf{x}, z) = \rho\_5 \mathbf{x} + \rho\_6 z + \rho\_7 > 0\_r$$

by considering *ρ*<sup>7</sup> > 0 and *x*(*t*), *z*(*t*) ≥ 0 from the domain **R** 3 <sup>+</sup>,0. Then,

$$y(t) = y(0)\exp\left[-\int\_0^t f(x,z)dt\right],$$

with *y*(0) ∈ **R**+,0. Therefore, all solutions with nonnegative initial conditions will be bounded as indicated below

$$K\_y = \{0 \le y(t) \le y\_{\text{sup}} = y(0)\},$$

hence, any upper bound for *x*(*t*) and *z*(*t*) depending on *K<sup>y</sup>* will be directly related to the glucose initial concentration [*y*(0)], which is expected as biomass and ethanol production over time is directly related to glucose dynamics.

Now, let us provide the mathematical background that allows us to compute a localizing domain where all compact invariant sets of a nonlinear dynamical system are located. The General Theorem concerning the LCIS method was formalized by Krishchenko and Starkov (see Section 2 in [58]) and it states the following: *Each compact invariant set* Γ *of x*˙ = *f*(*x*) *is contained in the localizing domain:*

$$K(h) = \left\{ h\_{\inf} \le h(\mathfrak{x}) \le h\_{\sup} \right\}.$$

From the latter *f*(*x*) is a *C* <sup>∞</sup>−differentiable vector function where *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>** *n* is the state vector. *h*(*x*) : **R** *<sup>n</sup>* <sup>→</sup> **<sup>R</sup>** is a *<sup>C</sup>* <sup>∞</sup>−differentiable function called localizing function, *<sup>h</sup>*|*<sup>S</sup>* denotes the restriction of *h*(*x*) on a set *S* ⊂ **R** *<sup>n</sup>* with *<sup>S</sup>*(*h*) = <sup>n</sup> *x* ∈ **R** *n* | *L<sup>f</sup> h*(*x*) = 0 o , and *L<sup>f</sup> h*(*x*) = (*∂h*/*∂x*)*f*(*x*) is the Lie derivative of *f*(*x*). Hence, one can define *h*inf = inf{*h*(*x*) | *x* ∈ *S*(*h*)} and *h*sup = sup{*h*(*x*) | *x* ∈ *S*(*h*)}. Furthermore, if all compact invariant sets are contained in the set *K*(*hi*) and in the set *K hj* then they are contained in *K*(*hi*) ∩*K hj* as well. The nonexistence of compact invariant sets can be considered for a given set Λ ⊂ **R** *n* if Λ ∩ *K*(*h*) = ∅, then the system *x*˙ = *f*(*x*) has no compact invariant sets located in Λ.

Following the LCIS method, one can explore the next localizing function

$$h\_1 = \mathfrak{x} + \mathfrak{a}y; \mathfrak{a} > 0,$$

then, the Lie derivative may be written as follows

$$L\_f h\_1 = -\rho\_4 \text{x} - \alpha \rho\_7 y - \frac{\alpha \rho\_5 y + \alpha \rho\_2 \rho\_5 - \rho\_1}{\rho\_2 + y} \text{x} y - \rho\_3 \text{x} \text{z} - \alpha \rho\_6 y \text{z} \text{y}$$

and the set *S*(*h*1) = n *L<sup>f</sup> h*<sup>1</sup> = 0 o is given by

$$S(h\_1) = \left\{\rho\_4 x = -\alpha \rho\_7 y - \frac{\alpha \rho\_5 y + \alpha \rho\_2 \rho\_5 - \rho\_1}{\rho\_2 + y} xy - \rho\_3 xz - \alpha \rho\_6 yz\right\},$$

where *x* = *h*<sup>1</sup> − *αy*, therefore set *S*(*h*1) is rewritten as indicated below

$$S(h\_1) = \left\{ h\_1 = \frac{\alpha(\rho\_4 - \rho\_7)}{\rho\_4} y - \frac{\alpha \rho\_5 y + \alpha \rho\_2 \rho\_5 - \rho\_1}{\rho\_4 (\rho\_2 + y)} xy - \frac{\rho\_3}{\rho\_4} xz - \frac{\alpha \rho\_6}{\rho\_4} yz \right\},$$

and the next two conditions are formulated

$$
\rho\_4 - \rho\_7 \quad > \quad 0,\tag{9}
$$

$$a\_1 > \frac{\rho\_1}{\rho\_2 \rho\_5},\tag{10}$$

where (9) is directly fulfilled by (4). Now, let us apply the Iterative Theorem in order to find an upper bound for the localizing function

$$S(h\_1) \cap K\_{\mathcal{Y}} \subset \left\{ h\_1 \le \frac{\alpha(\rho\_4 - \rho\_7)}{\rho\_4} y\_{\sup} \right\},$$

then

$$K(h\_1) = \left\{ \mathfrak{x}(t) + \mathfrak{a}y(t) \le \frac{\mathfrak{a}(\rho\_4 - \rho\_7)}{\rho\_4} y\_{\text{sup}} \right\} \rho\_4$$

from the latter, the upper bound for the biomass concentration [*x*(*t*)] may be written in terms of the parameters and the initial glucose concentration [*y*(0)] as follows

$$K\_{\mathbf{x}} = \left\{ 0 \le \mathbf{x}(t) \le \mathbf{x}\_{\text{sup}} = \frac{\mathfrak{a}(\rho\_4 - \rho\_7)}{\rho\_4} y\_{\text{sup}} \right\}.$$

Now, an upper bound for the ethanol concentration [*z*(*t*)] can be determined by the following localizing function

$$h\_2 = \beta\_1 \mathfrak{x} + \beta\_2 \mathfrak{y} + \underline{z}; \ \beta\_1, \beta\_2 > 0,$$

whose Lie derivative is computed as indicated below

$$\mathbf{L}\_f \mathbf{h}\_2 = -\beta\_1 \rho\_4 \mathbf{x} - \beta\_2 \rho\_7 \mathbf{y} - \rho\_{10} \mathbf{z} - \frac{\beta\_2 \rho\_5 \mathbf{y} + \beta\_2 \rho\_2 \rho\_5 - \beta\_1 \rho\_1}{\rho\_2 + \mathbf{y}} \mathbf{x} \mathbf{y} - (\beta\_1 \rho\_3 - \rho\_8) \mathbf{x} \mathbf{z} - (\beta\_2 \rho\_6 - \rho\_9) \mathbf{y} \mathbf{z}$$

and at this step, the following conditions are formulated

$$
\begin{array}{ccccc}
\rho\_1 & > & \frac{\rho\_8}{\rho\_3} \\
& & . & . \\
\end{array}
\tag{11}
$$

$$\beta\_2 \quad > \quad \max \left\{ \frac{\rho\_9}{\rho\_6}, \frac{\beta\_1 \rho\_1}{\rho\_2 \rho\_5} \right\}\_{\prime} \tag{12}$$

then, set *S*(*h*2) = n *L<sup>f</sup> h*<sup>2</sup> = 0 o , can be written as follows

$$S(h\_2) = \left\{\rho\_{10}z = -\beta\_1\rho\_4\mathbf{x} - \beta\_2\rho\_7\mathbf{y} - \frac{\beta\_2\rho\_5\mathbf{y} + \beta\_2\rho\_2\rho\_5 - \beta\_1\rho\_1}{\rho\_2 + \mathbf{y}}\mathbf{x}\mathbf{y} - (\beta\_1\rho\_3 - \rho\_8)\mathbf{x}\mathbf{z} - (\beta\_2\rho\_6 - \rho\_9)\mathbf{y}\mathbf{z}\right\}\rho\_1$$
 
$$\text{hence } \text{ ac } z = \mathbf{b} \quad \text{ } \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ } \quad \text{ }$$

hence, as *z* = *h*<sup>2</sup> − *β*1*x* − *β*2*y*, then set *S*(*h*2) is rewritten as indicated below

$$S(h\_2) = \left\{ h\_2 = \frac{\beta\_1(\rho\_{10} - \rho\_4)}{\rho\_{10}} \mathbf{x} + \frac{\beta\_2(\rho\_{10} - \rho\_7)}{\rho\_{10}} \mathbf{y} - \frac{\beta\_2\rho\_5\mathbf{y} + \beta\_2\rho\_2\rho\_5 - \beta\_1\rho\_1}{\rho\_2 + \mathbf{y}} \mathbf{x} - (\beta\_1\rho\_3 - \rho\_8)\mathbf{x} - (\beta\_2\rho\_6 - \rho\_9)\mathbf{y} \right\}.$$

where the next condition is formulated

$$
\rho\_{10} > \max \{ \rho\_{4'} \rho\_7 \},
\tag{13}
$$

and it holds by (4). Then, the Iterative Theorem is applied to get the following result

$$S(h\_2) \cap K\_\mathbf{x} \cap K\_\mathbf{y} \subset \left\{ h\_2 \le \frac{\beta\_1(\rho\_{10} - \rho\_4)}{\rho\_{10}} \mathbf{x}\_{\text{sup}} + \frac{\beta\_2(\rho\_{10} - \rho\_7)}{\rho\_{10}} y\_{\text{sup}} \right\} \rho\_{\text{sup}}$$

then, the upper bound for the localizing function *h*<sup>2</sup> is derived as follows

$$K(h\_2) = \left\{ \beta\_1 x(t) + \beta\_2 y(t) + z(t) \le \frac{\beta\_1 (\rho\_{10} - \rho\_4)}{\rho\_{10}} x\_{\text{sup}} + \frac{\beta\_2 (\rho\_{10} - \rho\_7)}{\rho\_{10}} y\_{\text{sup}} \right\}.$$

now, from the latter one can get the upper bound for ethanol concentration [*z*(*t*)] over time in terms of the parameters, the initial glucose concentration [*y*(0)] and the upper bound of biomass - *x*sup as given below

$$K\_z = \left\{ 0 \le z(t) \le z\_{\rm sup} = \frac{\beta\_1(\rho\_{10} - \rho\_4)}{\rho\_{10}} x\_{\rm sup} + \frac{\beta\_2(\rho\_{10} - \rho\_7)}{\rho\_{10}} y\_{\rm sup} \right\}.$$

Results shown above allow us to conclude the following regarding the boundedness of the KM system (1)–(3) solutions:

**Theorem 1.** *Localizing domain. If conditions (9)–(13) are fulfilled, then all compact invariant sets of the KM mechanistic model (1)–(3) are located either at the boundaries or within the following domain*

$$\mathcal{K}\_{\Gamma} = \mathcal{K}\_{\mathfrak{x}} \cap \mathcal{K}\_{\mathfrak{y}} \cap \mathcal{K}\_{\mathfrak{z}\mathcal{A}}$$

*where K*<sup>Γ</sup> ⊂ **R** 3 <sup>+</sup>,0*, and the ultimate bounds for biomass* [*x*(*t*)]*, glucose* [*y*(*t*)]*, and ethanol* [*z*(*t*)] *concentrations over time are given below*

$$\begin{aligned} \mathsf{K}\_{\mathsf{X}} &= \left\{ 0 \le x(t) \le x\_{\text{sup}} = \frac{\mathfrak{a}(\rho\_{4} - \rho\_{7})}{\rho\_{4}} y\_{\text{sup}} \right\}, \\\ \mathsf{K}\_{\mathsf{Y}} &= \left\{ 0 \le y(t) \le y\_{\text{sup}} = y(0) \right\}, \\\ \mathsf{K}\_{\mathsf{Z}} &= \left\{ 0 \le z(t) \le z\_{\text{sup}} = \frac{\mathfrak{f}\_{1}(\rho\_{10} - \rho\_{4})}{\rho\_{10}} x\_{\text{sup}} + \frac{\mathfrak{f}\_{2}(\rho\_{10} - \rho\_{7})}{\rho\_{10}} y\_{\text{sup}} \right\}. \end{aligned}$$

Now, let us briefly provide the mathematical background concerning the stability theory in the sense of Lyapunov, particularly the direct method where it is necessary to formulate a Lyapunov candidate function, which is usually denoted as *V*(*x*) : **R** *<sup>n</sup>* <sup>→</sup> **<sup>R</sup>**, a continuously differentiable function whose temporal derivative is given by *V*˙ (*x*) = (*∂V*/*∂x*)*f*(*x*). This function must be positive definite, i.e., *V*(0) = 0 and *V*(*x*) > 0 for *x* 6= 0, whilst a negative definite function is also *V*(0) = 0 but *V*(*x*) < 0 for *x* 6= 0. Further, function *V*(*x*) is said to be radially unbounded if *V*(*x*) → ∞ as k*x*k → ∞. The latter allows the formulation of the Global Asymptotic Stability Theorem (see Chapter 4 in [59] and Chapter 2 in [60]) which states the following: *The equilibrium point x* ∗ *is globally asymptotically stable if there exists a function V*(*x*) *positive definite, radially unbounded and decrescent such that its temporal derivative V*˙ (*x*) *is negative definite.* A function *V*(*x*) satisfying properties of this theorem is called Lyapunov function.

Following the latter, let us explore the next Lyapunov candidate function

$$V(\mathbf{x}, y, z) = \gamma\_1 \mathbf{x} + \gamma\_2 y + z\_\star$$

with

$$\gamma\_1, \gamma\_2 > 0,$$

then, the time derivative is computed as shown below

$$\dot{V}(\mathbf{x}, y, z) = \gamma\_1 \left(\frac{\rho\_1 \mathbf{x} y}{\rho\_2 + y} - \rho\_3 \mathbf{x} z - \rho\_4 \mathbf{x}\right) - \gamma\_2 (\rho\_5 \mathbf{x} y + \rho\_6 y z + \rho\_7 y) + \rho\_8 \mathbf{x} z + \rho\_9 y z - \rho\_{10} z \mathbf{x}$$

which can be rewritten as follows

$$\mathcal{V}(\mathbf{x}, y, z) = -\gamma\_1 \rho\_4 \mathbf{x} - \gamma\_2 \rho\_7 y - \rho\_{10} z - (\gamma\_1 \rho\_5 - \rho\_8) \mathbf{x} \mathbf{z} - (\gamma\_2 \rho\_6 - \rho\_9) y \mathbf{z} - \frac{y \gamma\_2 \rho\_5 + \gamma\_2 \rho\_2 \rho\_5 - \gamma\_1 \rho\_1}{\rho\_2 + y} \mathbf{x} y$$

where it is evident that *V*˙ (0, 0, 0) = 0, therefore the following constraints on coefficients *γ*<sup>1</sup> and *γ*<sup>2</sup> are formulated to ensure *V*˙ (*x*, *y*, *z*) < 0

$$
\gamma\_1 \quad > \quad \frac{\rho\_8}{\rho\_3} \tag{14}
$$

$$\gamma\_2 \quad > \quad \max \left\{ \frac{\rho\_9}{\rho\_6}, \frac{\gamma\_1 \rho\_1}{\rho\_2 \rho\_5} \right\}\_{\prime} \tag{15}$$

thus, as parameters *ρ<sup>i</sup>* , *i* = 1, 2, 3, 5, 6, 8, 9; in both conditions are different for each term, then it is possible to assume that there exists a set of solutions that satisfies (14) and (15). Hence, the following result can be concluded:

**Theorem 2.** *Asymptotic stability. If conditions (14) and (15) are fulfilled, then the KM mechanistic model (1)–(3) is asymptotically stable and all trajectories will go to the equilibrium point* (*x* ∗ 0 , *y* ∗ 0 , *z* ∗ 0 ) = (0, 0, 0).

The latter implies that any given trajectory [*φ*(*x*(*t*), *y*(*t*), *z*(*t*))] with nonnegative initial conditions [*x*(0), *<sup>y</sup>*(0), *<sup>z</sup>*(0) <sup>≥</sup> <sup>0</sup>] passing through any point (*x*(*t*), *<sup>y</sup>*(*t*), *<sup>z</sup>*(*t*))*<sup>T</sup>* in **R** 3 <sup>+</sup>,0 its *ω*−limit set is not empty and it is a compact invariant set, i.e.,

$$\lim\_{t \to \infty} \phi(x(t), y(t), z(t)) = (0, 0, 0)^T \nu$$

see Lemma 4.1 by Khalil in [59] at Section 4.2 and Theorem 1 by Perko in [61] at Section 3.2.

Concerning the existence and uniqueness of solutions for the KM system (1)–(3), let us introduce the following notations for the sake of simplicity

$$\begin{array}{rcl} f\_1(t,x,y,z) &=& \frac{\rho\_1 x y}{\rho\_2 + y} - \rho\_3 x z - \rho\_4 x, \\ f\_2(t,x,y,z) &=& -\rho\_5 x y - \rho\_6 y z - \rho\_7 y, \\ f\_3(t,x,y,z) &=& \rho\_8 x z + \rho\_9 y z - \rho\_{10} z, \end{array}$$

and compute the Jacobian matrix [*∂ f* /*∂u*](*t*, *u*) (see [49] at Section 7.4) with results shown below for *fi*(*t*, *u*), *i* = 1, 2, 3; and *u* = [*x*, *y*, *z*] *T*

$$J = \begin{bmatrix} \frac{\rho\_1 y}{\rho\_2 + y} - \rho\_3 z - \rho\_4 & \frac{\rho\_1 \rho\_2 x}{\left(\rho\_2 + y\right)^2} & -\rho\_3 x \\ -\rho\_5 y & -\rho\_5 x - \rho\_6 z - \rho\_7 & -\rho\_6 y \\ \rho\_8 z & \rho\_9 z & \rho\_8 x + \rho\_9 y - \rho\_{10} \end{bmatrix},\tag{16}$$

and it is evident that *fi*(*t*, *u*) and [*∂ f* /*∂u*](*t*, *u*) are continuous and exist on the domain Ω = h *t*0, *t<sup>f</sup>* i <sup>×</sup> *<sup>K</sup>*<sup>Γ</sup> with <sup>h</sup> *t*0, *t<sup>f</sup>* i ∈ [*t*0, ∞] and *K*<sup>Γ</sup> ⊂ **R** 3 <sup>+</sup>,0 [33]. Hence, the latter implies that *fi*(*t*, *u*) is locally Lipschitz in *u* on Ω (see Lemma 3.2 by Khalil in [59] at Section 3.1). Further, each element of (16) is bounded by Theorem 1. Thus, the following can be concluded:

**Theorem 3.** *Existence and uniqueness. There is a Lipschitz constant* ` ≥ 0 *such that* k[*∂ f* /*∂u*](*t*, *u*)k ≤ ` *on* Ω. *Then, fi*(*t*, *u*) *satisfies the Lipschitz condition*

$$\|f(t, u\_1) - f(t, u\_2)\| \le \ell \|u\_1 - u\_2\|\_{\varkappa}$$

*and there exists some δ* > 0 *such that the KM mechanistic model (1)–(3), given as u*˙ = *fi*(*t*, *u*) *with u*(*t*0) = *u*0*, has a unique solution over* [*t*0, *t*<sup>0</sup> + *δ*]*.*

Although conditions for asymptotic stability of the equilibrium point(*x* ∗ 0 , *y* ∗ 0 , *z* ∗ 0 ) = (0, 0, 0) in **R** 3 <sup>+</sup>,0 were established in Theorem 2, it is straightforward to demonstrate its local asymptotic stability by evaluating (5) in (16) as follows

$$J\_{(0,0,0)} = \begin{bmatrix} -\rho\_4 & 0 & 0\\ 0 & -\rho\_7 & 0\\ 0 & 0 & -\rho\_{10} \end{bmatrix} / $$

where the eigenvalues [*λ<sup>i</sup>* , *i* = 1, 2, 3] are given by each element of the diagonal. Thus, *λ*<sup>1</sup> = −*ρ*4, *λ*<sup>2</sup> = −*ρ*7, and *λ*<sup>3</sup> = −*ρ*10. Therefore, Theorem 4.7 by Khalil in [59] allows us to conclude the next additional result to Theorem 2:

**Corollary 1.** *Local stability. The equilibrium point* (*x* ∗ 0 , *y* ∗ 0 , *z* ∗ 0 ) = (0, 0, 0) *of the KM mechanistic model (1)–(3) is locally asymptotically stable in* **R** 3 <sup>+</sup>,0.

### **4. Discussion**

The KM mechanistic model (1)–(3) was formulated by considering the biological relationships between each variable in a controlled batch fermentation where concentrations in g/L were measured for biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)] over 13 consecutive *hours*. Then, by means of the *lsqcurvefit* function, an algorithm was developed in Matlab to approximate the experimental data from the 17 *K. marxianus* strains discussed at Section 2; both qualitative (see Figures 2–6) and quantitative (see Tables 3 and 4) results were shown in Section 3. The in silico experimentation illustrates the capabilities of the system to approximate the experimental data of each strain, whilst both the *R* <sup>2</sup> and the *AIC* provide a value for the goodness of fit of the model to each set of data. In Table 4, one

can see that *R* <sup>2</sup> values range from 0.955 to 0.994, and *AIC* from <sup>−</sup>43.478 to 33.184, these values are for strains 7 and 9, respectively.

Now, it should be noted that the dynamics between biomass growth, substrate consumption and product generation have been modeled before by means of the logistic growth law [62], the Pirt Equation [63], and Luedeking–Piret Equation [64] as indicated below in Equations (17)–(19), respectively:

$$\dot{X} \quad = \quad \mu\_{\text{max}} X \Big(1 - \frac{X}{X\_{\text{max}}}\Big) \,\tag{17}$$

$$
\dot{S}\_{\perp} = -\frac{1}{Y\_{X/S}} \dot{X} - mX\_{\prime} \tag{18}
$$

$$
\dot{P}\_{\quad} = \ \ a \dot{X} + \beta X\_{\text{\\_}} \tag{19}
$$

where *µ*max is the biomass maximum growth rate, this parameter is equivalent to *ρ*<sup>1</sup> in our mathematical model; *X*max the maximum concentration value of biomass in the experimental data set for the time-interval of the process being observed; *YX*/*<sup>S</sup>* the biomass/substrate yield; *m* is the maintenance coefficient; *α* is the growth-associated coefficient for the product; and *β* is the non-growth-associated coefficient for the product. Our algorithm was applied to approximate the experimental data of the 17 *K. marxianus* strains with overall results shown in Table 5.

**Table 5.** The logistic, Pir, and Luedeking–Piret Equations (17)–(19) provides valuable information concerning biomass growth [*µ*max], biomass/substrate yield [*YX*/*S*], and product generation [*α*]; estimated numerical values are given in their respective columns. Concerning the goodness of fit, results regarding the *RSS*, *AIC*, and *R* <sup>2</sup> are provided in the following columns.


The main comparison between the KM system (1)–(3) and Equations (17)–(19) is performed with respect to the biomass maximum growth rate, given by *ρ*<sup>1</sup> and *µ*max, respectively. Tables 2 and 5 show that estimated values of *ρ*<sup>1</sup> are on average ∼0.717 smaller than those estimated for *µ*max. The latter is a direct consequence of the biological assumptions on which each mechanistic model was formulated. The KM system (1)–(3) was constructed

by considering interactions between the three variables as illustrated in the flow diagram of Figure 1, whilst the logistic, Pirt, and Luedeking–Piret Equations (17)–(19) are constructed by only assuming a logistic growth for biomass without taking into account the overall effect of ethanol production over the entire system as well as the death rate of biomass [*x*(*t*)], decomposition rate of glucose [*y*(*t*)], and degradation of ethanol [*z*(*t*)]. Further, the in silico experimentation concerning Equations (17)–(19) illustrated in Figures A1–A5 at Appendix B shows that approximated values for substrate [*S*(*t*)], i.e., glucose, becomes negative as time increases, which is not biologically possible for this variable. Further, one can see from the experimental data that ethanol production does not follow a smooth sigmoidal growth, the data even illustrates degradation among some strains, which is better approximated by our model as it is shown in the lower panels of Figures 2–6.

When comparing the goodness of fit by computing the *AIC* and *R* 2 , it is evident that the KM system (1)–(3) had overall better results than the logistic, Pirt, and Luedeking–Piret Equations (17)–(19). Although the latter has fewer parameters than ours (six and ten, respectively) and the *AIC* penalizes a model with more parameters to be fitted, results for the *RSS* were lower for the KM system which ultimately worked in our favor. Further, the capabilities of the KM mechanistic model may extend beyond its ability to approximate experimental data and estimate the biomass maximum growth rate, in Appendix C the in silico experimentation illustrates the dynamics for *t* ∈ [0, 39], i.e., three times the period for the experimental data. Figures A6–A10 show that as time increases and the substrate is no longer added into the system, then the death of biomass and degradation of ethanol begins to take over the system. The latter was expected from the asymptotic stability results of Section 3, particularly Theorem 2 and Corollary 1, as these state that the concentration of all variables will eventually be zero, i.e., both biomass [*x*(*t*)] and ethanol [*z*(*t*)] concentrations are going to be depleted. Additionally, it is important to note that all solutions of the KM system are bounded from above, which is consistent with the localizing domain results of Theorem 1.

Regarding the values of parameters *m* and *β*, our algorithm yielded results in the magnitude of 10−<sup>14</sup> for *m* in all strains; in fact, setting *m* to zero does not affect the ultimate results for the other parameters [*µ*max,*YX*/*S*, *α*, and *β*] which may allow us to completely disregard this term [−*mX*] from Equations (17)–(19). Concerning *β*, values for 12 strains were in the same order of magnitude - 10−14 , however, the following results were determined for strains 3, 4, 11, 12, and 17: 87.633 <sup>×</sup> <sup>10</sup>−<sup>3</sup> , 38.660 <sup>×</sup> <sup>10</sup>−<sup>3</sup> , 51.050 <sup>×</sup> <sup>10</sup>−<sup>3</sup> , 53.813 <sup>×</sup> <sup>10</sup>−<sup>3</sup> , 37.702 <sup>×</sup> <sup>10</sup>−<sup>3</sup> , respectively. Hence, the non-growth-associated coefficient for the product may influence the dynamics in some karyotypes of *K. marxianus*.

### **5. Conclusions**

Mechanistic modeling has proven to be a powerful tool capable of describing the relationships between different variables in the dynamics of biological systems when considering assumptions based on scientific principles underlying the phenomenon being modeled. In this work, a set of three coupled first-order ODEs was formulated which can approximate experimental changes over time of alcoholic fermentation in batch culture by 17 different strains of *K. marxianus*.

The KM mechanistic model (1)–(3) describes biomass growth [*x*(*t*)], glucose consumption [*y*(*t*)], and ethanol production [*z*(*t*)] in concentrations of g/L per hour. The parameter values of the system were estimated through a nonlinear curve-fitting algorithm in Matlab with the experimental data of each batch culture fermentation described in Section 2. The latter allowed us to conclude that seven parameters have the same numerical value for the dynamics observed in the 17 strains, particularly the affinity with substrate constant [*ρ*2], inhibition rate of biomass growth due to product accumulation [*ρ*3], biomass death rate [*ρ*4], consumption rates for biomass growth and ethanol production [*ρ*<sup>5</sup> and *ρ*6], glucose spontaneous decomposition rate [*ρ*7], and ethanol degradation rate [*ρ*10]; these values are shown in Table 1. However, the biomass maximum growth rate [*ρ*1], ethanol production associated with biomass growth [*ρ*8], and glucose converted in ethanol [*ρ*9] parameters

have specific values for each strain, results are shown in Table 2 with a 95% confidence interval that gives us the margin of error for each parameter value estimation.

As predictive microbiology establishes, mathematical models must be simplified until measurable parameters can be obtained, the KM mechanistic model successfully achieves this with *ρ*1, *ρ*8, and *ρ*<sup>9</sup> as the main parameters that describe the overall dynamics of the batch fermentation process under study in this research. The biomass growth rate is a very specific value for each strain that must be as high as possible. Ethanol production with respect to biomass growth represents the fermentative capacity of each strain, and the concentration of glucose converted to ethanol is directly related to these rates. It should be noted that in batch culture the latter requires high sugar concentrations to achieve alcoholic fermentation.

Further, the in silico experimentation illustrates that our model may be able to accurately predict the concentration of each variable as it is shown in Appendix C; nonetheless, further experimental data are needed to properly validate this assessment. One can see in Figures A6–A10 that when no more substrate is added to the culture, then biomass growth goes into the death phase, and ethanol degradation begins to happen in the system. This behavior is to be expected as the nonlinear analysis of the system allowed us to conclude that all concentrations will eventually go to zero in the absence of glucose, i.e., the asymptotic stability of the equilibrium (5) [(*x* ∗ 0 , *y* ∗ 0 , *z* ∗ 0 ) = (0, 0, 0)] by Theorem 2 and Corollary 1; further, concentrations over time of all variables are bounded by the Localizing Domain Theorem 1. The latter is illustrated in all panels for the predictions of biomass growth [*x*(*t*)], glucose consumption [*y*(*t*)], and ethanol production [*z*(*t*)].

Finally, the KM mechanistic model may be useful in the field of predictive microbiology, particularly in alcoholic fermentation through yeast and sugar, such as *K. marxianus* and glucose as only three parameters of our system needs to be fitted for different strains. Furthermore, when comparing the results of the biomass maximum growth rate of our model with the classic logistic, Pirt, and Luedeking–Piret Equations (17)–(19), our values are on average 71.7% smaller as the KM system (1)–(3) takes into account the effect of both substrate and product concentrations in the batch culture over the biomass growth phases.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/e25030497/s1. this algorithm can be applied to fit the numerical values of parameters [*ρ*<sup>1</sup> , *ρ*8, *ρ*9], and to determine results concerning the statistics and goodness of fit.

**Author Contributions:** Conceptualization, Y.S. and P.A.V.; methodology, J.B.P.-L., N.O.S.-C. and P.A.V.; software, P.A.V.; validation, J.B.P.-L. and N.O.S.-C.; formal analysis, P.A.V. and E.R.; investigation, Y.S. and E.R.; resources, J.B.P.-L., N.O.S.-C. and P.A.V.; data curation, J.B.P.-L., N.O.S.-C. and F.J.R.-S.; writing—original draft preparation, P.A.V., Y.S. and E.R.; writing—review and editing, P.A.V., Y.S., J.B.P.-L. and N.O.S.-C.; visualization, P.A.V. and Y.S.; supervision, Y.S. and P.A.V.; project administration, Y.S. and P.A.V.; funding acquisition, Y.S., P.A.V., J.B.P.-L. and N.O.S.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** Authors of this work were supported by different projects indicated as follows: Yolocuauhtli Salazar by the TecNM project 14166.22-P: Modelizado de la producción y vida de anaquel de microorganismos en productos fermentados; Paul A. Valle by the TecNM project: Gemelos digitales para el análisis y control de sistemas biológicos; Nicolás O. Soto-Cruz and Jesús B. Páez-Lerma by the CONACyT project Ciencia Básica 252465, and the TecNM project 10197-21-P.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data supporting the findings of the present work are available through authors Nicolás O. Soto-Cruz and Jesús B. Páez-Lerma upon reasonable request.

**Acknowledgments:** Emmanuel Rodríguez and Francisco J. Reyes-Sánchez express their gratitude to CONACyT for the financial support during their postgraduate studies under their scholarship agreements 813941 and 465646, respectively.

**Conflicts of Interest:** The authors declare no conflict of interest. Funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A. Equilibria**

In order to further discuss the nonlinear mathematical analysis of Section 3, all equilibrium points of the system are calculated and it becomes evident that the origin is the only biologically meaningful result. The set of equilibria of the KM system (1)–(3) is determined by solving the next system of equations

$$\begin{aligned} f\_1(t, \mathbf{x}, y, z) &= \, \, \varkappa \left( \frac{\rho\_1 y}{\rho\_2 + y} - \rho\_3 z - \rho\_4 \right) = 0, \\ f\_2(t, \mathbf{x}, y, z) &= \, \, \, -y(\rho\_5 \mathbf{x} + \rho\_6 z + \rho\_7) = 0, \\ f\_3(t, \mathbf{x}, y, z) &= \, \, \, z(\rho\_8 \mathbf{x} + \rho\_9 y - \rho\_{10}) = 0, \end{aligned}$$

from which one can compute the following

$$\begin{aligned} \left(x\_{0'}^\* y\_{0'}^\* z\_0^\*\right) &=& \left(0, 0, 0\right), \\ \left(x\_{1'}^\* y\_{1'}^\* z\_1^\*\right) &=& \left(0, \frac{\rho\_{10}}{\rho\_9}, -\frac{\rho\_7}{\rho\_6}\right), \\ \left(x\_{2'}^\* y\_{2'}^\* z\_2^\*\right) &=& \left(\frac{\rho\_{10}}{\rho\_8}, 0, -\frac{\rho\_4}{\rho\_3}\right), \\ \left(x\_{3'}^\* y\_{3'}^\* z\_3^\*\right) &=& \left(-\frac{\rho\_7}{\rho\_5}, \frac{\rho\_2 \rho\_4}{\rho\_1 - \rho\_4}, 0\right). \end{aligned}$$

and

$$\begin{aligned} (\mathfrak{x}\_{4\prime}^\* y\_{4\prime}^\* z\_4^\*)\_{\prime} \\ (\mathfrak{x}\_{5\prime}^\* y\_{5\prime}^\* z\_5^\*)\_{\prime} \end{aligned}$$

where

$$\begin{array}{rcl}x\_4^\* &=& \frac{-\rho\_1\rho\_6\rho\_8 + \rho\_3\rho\_5\rho\_{10} - \rho\_3\rho\_7\rho\_8 + \rho\_4\rho\_6\rho\_8 + \rho\_2\rho\_3\rho\_5\rho\_9 - \sqrt{\rho}}{2\rho\_3\rho\_5\rho\_8},\\y\_4^\* &=& \frac{\rho\_1\rho\_6\rho\_8 + \rho\_3\rho\_5\rho\_{10} + \rho\_3\rho\_7\rho\_8 - \rho\_4\rho\_6\rho\_8 - \rho\_2\rho\_3\rho\_5\rho\_9 + \sqrt{\rho}}{2\rho\_3\rho\_5\rho\_9},\\z\_4^\* &=& \frac{\rho\_1\rho\_6\rho\_8 - \rho\_3\rho\_5\rho\_{10} - \rho\_3\rho\_7\rho\_8 - \rho\_4\rho\_6\rho\_8 - \rho\_2\rho\_3\rho\_5\rho\_9 + \sqrt{\rho}}{2\rho\_3\rho\_6\rho\_8},\end{array}$$

$$\begin{aligned} x\_5^{\*} &= \frac{-\rho\_1 \rho\_6 \rho\_8 + \rho\_3 \rho\_5 \rho\_{10} - \rho\_3 \rho\_7 \rho\_8 + \rho\_4 \rho\_6 \rho\_8 + \rho\_2 \rho\_3 \rho\_5 \rho\_9 + \sqrt{\rho}}{2 \rho\_3 \rho\_5 \rho\_8}, \\ y\_5^{\*} &= \frac{\rho\_1 \rho\_6 \rho\_8 + \rho\_3 \rho\_5 \rho\_{10} + \rho\_3 \rho\_7 \rho\_8 - \rho\_4 \rho\_6 \rho\_8 - \rho\_2 \rho\_3 \rho\_5 \rho\_9 - \sqrt{\rho}}{2 \rho\_3 \rho\_5 \rho\_9}, \\ z\_5^{\*} &= \frac{\rho\_1 \rho\_6 \rho\_8 - \rho\_3 \rho\_5 \rho\_{10} - \rho\_3 \rho\_7 \rho\_8 - \rho\_4 \rho\_6 \rho\_8 - \rho\_2 \rho\_3 \rho\_5 \rho\_9 - \sqrt{\rho}}{2 \rho\_3 \rho\_6 \rho\_8}, \end{aligned}$$

with

*ρ* = (*ρ*3*ρ*7*ρ*<sup>8</sup> + *ρ*1*ρ*6*ρ*<sup>8</sup> − *ρ*4*ρ*6*ρ*<sup>8</sup> − *ρ*3*ρ*5*ρ*<sup>10</sup> − *ρ*2*ρ*3*ρ*5*ρ*9) <sup>2</sup> <sup>+</sup> <sup>4</sup>*ρ*3*ρ*5*ρ*8(*ρ*1*ρ*6*ρ*<sup>10</sup> <sup>+</sup> (*ρ*<sup>10</sup> <sup>+</sup> *<sup>ρ</sup>*2*ρ*9)(*ρ*3*ρ*<sup>7</sup> <sup>−</sup> *<sup>ρ</sup>*4*ρ*6)).

> Now, it is evident that the equilibrium points *x* ∗ *i* , *y* ∗ *i* , *z* ∗ *i* , *i* = 1, 2, 3; have at least one negative term. However, although is not straightforward, the same can be concluded regarding equilibriums *x* ∗ 4 , *y* ∗ 4 , *z* ∗ 4 and (*x* ∗ 5 , *y* ∗ 5 , *z* ∗ 5 ), as these are computed by disregarding the common term in each equation as follows

$$\begin{aligned} \frac{\rho\_1 y}{\rho\_2 + y} - \rho\_3 z - \rho\_4 &= 0, \\ \rho\_5 x + \rho\_6 z + \rho\_7 &= 0, \\ \rho\_8 x + \rho\_9 y - \rho\_{10} &= 0, \end{aligned}$$

and equality *ρ*5*x* + *ρ*6*z* + *ρ*<sup>7</sup> = 0 can only be fulfilled when either *x* ∗ *j* or *z* ∗ *j* , *j* = 4, 5; are negative. Therefore, the KM mechanistic model (1)–(3) has unique biologically meaningful equilibrium given by

$$(x\_0^\*, y\_0^\*, z\_0^\*) = (0, 0, 0).$$

### **Appendix B. Logistic, Pirt, and Luedeking–Piret Equations**

This appendix presents results concerning the in silico experimentation when fitting the experimental data to the logistic, Pirt, and Luedeking–Piret Equations (17)–(19). Figures A1–A5 are aiming to qualitative compare the proposed mathematical model with the classic model of biomass-substrate-product. Further, a quantitative comparison is carried out through the coefficient of determination and the Akaike Information Criterion at Section 4, see Tables 2, 4, and 5.

] L=g[ )t( X

12

14

12

14

12

14

] L=g[ )t( P

] L=g[ )t( S

**Figure A2.** Observed data (×green marker), and approximated values (continuous red line) for strains 5–8 with the Logistic, Pirt, and Luedeking–Piret Equations (17)–(19).

*Entropy* **2023**, *25*, 497

*Entropy* **2023**, *25*, 497

**Figure A4.** Observed data (×green marker), and approximated values (continuous red line) for strains 13–16 with the Logistic, Pirt, and Luedeking–Piret Equations (17)–(19).

*Entropy* **2023**, *25*, 497

**Figure A5.** Observed data (×green marker), and approximated values (continuous red line) for strain 17 with the Logistic, Pirt, and Luedeking–Piret Equations (17)–(19).

### **Appendix C. Predictive Ability of the KM Mechanistic Model**

This appendix presents results concerning the in silico experimentation when solving the KM mechanistic model (1)–(3) for a time interval of *t* ∈ [0, 39] in order to illustrate its ability to predict the dynamics of the three variables, i.e., the concentration in g/L over time between biomass [*x*(*t*)], glucose [*y*(*t*)], and ethanol [*z*(*t*)], after the last experimental data point without further substrate addition into the batch. It should be noted that at this stage of the research there is not available data to validate if the model is be able to accurately predict the evolution of both biomass and ethanol in the system, further experimental data points could be helpful to better fit the values of biomass death rate [*ρ*4], and ethanol degradation rate [*ρ*10]. However, Figures A6–A10 allow us to illustrate results concerning Theorem 2 and Corollary 1 from the asymptotic stability analysis.

] L=g[ )t( x

] L=g[ )t( z

] L=g[ )t( y

] L=g[ )t( x

**Figure A8.** Dynamics prediction for strains 9–12 with the KM system (1)–(3). Observed data are given by the × green marker, and approximated values by the continuous blue line.

 z

] L=g[ )t( x

] L=g[ )t( z

] L=g[ )t( y

**Figure A10.** Dynamics prediction for strain 17 with the KM system (1)–(3). Observed data are given by the × green marker, and approximated values by the continuous blue line.

### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Review* **Mathematical Models of Death Signaling Networks**

**Madhumita Srinivasan <sup>1</sup> , Robert Clarke <sup>2</sup> and Pavel Kraikivski 3,\***


**Abstract:** This review provides an overview of the progress made by computational and systems biologists in characterizing different cell death regulatory mechanisms that constitute the cell death network. We define the cell death network as a comprehensive decision-making mechanism that controls multiple death execution molecular circuits. This network involves multiple feedback and feed-forward loops and crosstalk among different cell death-regulating pathways. While substantial progress has been made in characterizing individual cell death execution pathways, the cell death decision network is poorly defined and understood. Certainly, understanding the dynamic behavior of such complex regulatory mechanisms can be only achieved by applying mathematical modeling and system-oriented approaches. Here, we provide an overview of mathematical models that have been developed to characterize different cell death mechanisms and intend to identify future research directions in this field.

**Keywords:** cell death; apoptosis; necroptosis; ferroptosis; pyroptosis; immunogenic cell death; regulatory networks; death execution pathways; mathematical models; computational analysis

**Citation:** Srinivasan, M.; Clarke, R.; Kraikivski, P. Mathematical Models of Death Signaling Networks. *Entropy* **2022**, *24*, 1402. https:// doi.org/10.3390/e24101402

Academic Editor: Yoh Iwasa

Received: 7 September 2022 Accepted: 28 September 2022 Published: 1 October 2022

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

**Copyright:** © 2022 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/).

### **1. Introduction**

Mathematical modeling is a powerful tool that allows one to connect molecular biology to cell physiology by associating the qualitative and quantitative features of dynamical molecular networks with signal–response curves measured by cell biologists [1]. Mathematical and systems-oriented approaches have been successfully applied to describe the dynamics of complex molecular networks that control cell cycle [2,3], nutrient signaling [4], checkpoints [5], signaling dysregulation in cancer [6], and cell death [7–13]. Systemsoriented mathematical approaches are especially useful for analyzing complex systems that cannot be understood by intuitive reasoning. Undoubtedly, cell death regulation is one such molecular mechanism that cannot be fully understood without mathematical modeling. Here, we provide an overview of mathematical models that have been successfully applied to quantitatively characterize death signaling networks.

Cell death mechanisms are directly involved in regulations of tissue homeostasis, inflammation, immunity, development and other physiological processes [14]. Characterization of new genes and molecular components, involved in signaling pathways by regulating cell death, continues to progress. A detailed characterization of cell death regulation can help identify novel targets and develop effective therapeutic protocols to strike acquired drug resistance in cancer cells. Accurate predictive mechanistic models of complex molecular networks regulating cell death can be used to test the effects of new drugs on the system, and to search for synergistic drug combinations and effective treatment protocols. Different modeling approaches have been already successfully applied to model extensive cell death molecular networks. Ordinary differential equations (ODEs), Boolean logic, pharmacokinetic-pharmacodynamic (PK-PD), Petri nets, agentbased modeling (ABM), cellular automata and hybrid approaches are the common choices

available to model molecular mechanisms involved in cell death control, decisions and execution [6,12,13,15–19].

Cell death execution is an all-or-none, irreversible process [20]. Mathematically the activation of irreversible cell death can be described by an irreversible bistable switch with a stable survival steady state, a stable death steady state, and a third unstable steady state separating the survival and death states [21–24]. A pro-death signal can induce cell death by driving the bistable system from the survival to the death state. The transition occurs when the pro-death signal reaches a threshold that corresponds to the limit point bifurcation. Transition in the reverse direction, from death to survival, is impossible because the second limit point bifurcation, where the death steady state vanishes, occurs in the biologically irrelevant negative signal values (i.e., the concentration of a death-inducing ligand or stressor cannot be negative). Therefore, the activation of the cell death execution in such a bistable system cannot be reversed, even if the initial cell death trigger is removed. This mathematical description of the cell death activation is consistent with a threshold mechanism for cell death induction [25] and an all-or-none death decision [22,26,27]. Importantly, understanding how cells control the cell death/survival switch can help to identify targets that can force cancer cells to flip the switch to activate the irreversible cell death execution.

Complexity of cell death regulatory networks, a requirement to account for all important regulating molecular details and pathways, availability of merely small sets of sparse data for model calibration, as well as under- and over-fitting of the model are issues that must be routinely solved in order to develop a predictive model of cell death [6,15]. This review describes mathematical models that have been successfully applied to quantitatively characterize such cell death control mechanisms as apoptosis, necroptosis, ferroptosis, pyroptosis and immunogenic cell death.

### **2. Apoptosis**

Apoptosis is one of the most well-studied and characterized programmed cell death mechanisms. The detailed characterization of molecular interactions involved in apoptosis, and the growing amount of related quantitative data, has encouraged computational and systems biologists to develop mathematical models of apoptosis [12,13,17]. Over the last twenty years, several dozen mathematical models of apoptosis regulation have been described. These apoptosis models aim to explain different data or effects of different treatments on cell death. While the core molecular components regulating apoptosis are shared by all models, variations in molecular circuit designs, components, data, mathematical approaches, and study goals make each model a unique tool to study apoptosis. Most often, molecular mechanisms of apoptosis are mathematically represented using ODEs [7,21,22,25–33], Boolean logics [34–36], and Petri nets [16]; other computational approaches have also been applied [18,37,38].

The execution core of apoptosis regulation involves a family of proteases termed caspases. Caspases can be separated into the following two groups: effector or executioner caspases (caspase-3, -6, -7), and active initiator caspases (e.g., caspases-8, -9). Activation of the caspases initiates the cleavage of several important cellular proteins, such as actin and nuclear lamins, which results in cell body and nuclear shrinkage and cell death [39]. Apoptosis can be processed through mitochondria-dependent (intrinsic apoptosis) and mitochondria-independent (extrinsic apoptosis) caspase-3 activation pathways [14]. The core components involved in these two pathways are commonly included in all mathematical models of apoptosis and can be found in the earliest mathematical models of apoptosis [7].

Extrinsic apoptosis is characterized by high amounts of active caspase-8 that activates the downstream effectors caspase-3, caspase-6, and caspase-7. The activation of caspase-8 is receptor-mediated, which occurs upon receipt of a death signal that is processed by a surface death receptor such as FAS (a member of the tumor necrosis factor gene superfamily) [14]. Therefore, extrinsic apoptosis is a receptor-mediated cell death mechanism, as shown in Figure 1 (left). By contrast, intrinsic apoptosis can be executed even in cells with

lower levels of active caspase-8 but requires an additional amplification that involves activation of the pro-apoptotic functions of the mitochondria. For example, stress-related factors (e.g., DNA damage) can induce activation of the executioner caspases via a mitochondriadependent pathway in the absence of an external death signal [40] (Figure 1, right panel). The mitochondria-dependent pathway begins with the cleavage of anti-apoptotic Bcl-2 family members, which causes the aggregation of pro-apoptotic proteins such as Bax and Bak. Aggregation of pro-apoptotic proteins is followed by the release of cytochrome *c* from the mitochondria, which induces the formation of a large protein complex known as the apoptosome. The apoptosome recruits and activates caspase-9, allowing it to cleave the downstream effectors pro-caspase-3, pro-caspase-6, and pro-caspase-7. Notably, the expression of anti-apoptotic Bcl-2 family members can block the intrinsic apoptosis signaling in cells. By contrast, extrinsic apoptosis cannot be blocked by the expression of high levels of Bcl-2 proteins because large amounts of caspase-8 are already generated. involves activation of the pro-apoptotic functions of the mitochondria. For example, stress-related factors (e.g., DNA damage) can induce activation of the executioner caspases via a mitochondria-dependent pathway in the absence of an external death signal [40] (Figure 1, right panel). The mitochondria-dependent pathway begins with the cleavage of anti-apoptotic Bcl-2 family members, which causes the aggregation of pro-apoptotic proteins such as Bax and Bak. Aggregation of pro-apoptotic proteins is followed by the release of cytochrome *c* from the mitochondria, which induces the formation of a large protein complex known as the apoptosome. The apoptosome recruits and activates caspase-9, allowing it to cleave the downstream effectors pro-caspase-3, pro-caspase-6, and procaspase-7. Notably, the expression of anti-apoptotic Bcl-2 family members can block the intrinsic apoptosis signaling in cells. By contrast, extrinsic apoptosis cannot be blocked by the expression of high levels of Bcl-2 proteins because large amounts of caspase-8 are already generated.

superfamily) [14]. Therefore, extrinsic apoptosis is a receptor-mediated cell death mechanism, as shown in Figure 1 (left). By contrast, intrinsic apoptosis can be executed even in cells with lower levels of active caspase-8 but requires an additional amplification that

*Entropy* **2022**, *24*, x FOR PEER REVIEW 3 of 23

**Figure 1.** The mechanisms of receptor-induced apoptosis (left) and stress-mediated, mitochondriadependent apoptosis (right). Solid lines represent activation (arrowhead) and inhibition (bar head) influences. **Figure 1.** The mechanisms of receptor-induced apoptosis (left) and stress-mediated, mitochondriadependent apoptosis (right). Solid lines represent activation (arrowhead) and inhibition (bar head) influences.

The earliest mathematical models of apoptosis described both mitochondria-dependent and independent death activation pathways. In early 2000, Fussenegger et al. published a mechanistic ODE-based mathematical model of apoptosis that describes both receptor-mediated and stress-induced caspase activation mechanisms [7]. The receptor-mediated feature of the model describes the FAS surface receptor that activates procaspase-8. Activation of apoptosis initiator caspases involves the following reactions: the binding of an extracellular death ligand to the FAS receptor, the binding of FAS-associated death domain (FADD) protein to the FAS death domain, and the binding of caspase-8 to a domain on FADD that enables caspase-8 activation by proteolytic cleavage. Each binding process is described by a specific rate parameter in the model. Simulation results show that about 50% of procaspase-8 is activated within two hours after the death signal is received. After procaspase-8 activation, the executioner caspase is activated within minutes, The earliest mathematical models of apoptosis described both mitochondria-dependent and independent death activation pathways. In early 2000, Fussenegger et al. published a mechanistic ODE-based mathematical model of apoptosis that describes both receptormediated and stress-induced caspase activation mechanisms [7]. The receptor-mediated feature of the model describes the FAS surface receptor that activates procaspase-8. Activation of apoptosis initiator caspases involves the following reactions: the binding of an extracellular death ligand to the FAS receptor, the binding of FAS-associated death domain (FADD) protein to the FAS death domain, and the binding of caspase-8 to a domain on FADD that enables caspase-8 activation by proteolytic cleavage. Each binding process is described by a specific rate parameter in the model. Simulation results show that about 50% of procaspase-8 is activated within two hours after the death signal is received. After procaspase-8 activation, the executioner caspase is activated within minutes, and then the initiation of procaspase-9 occurs with the lag time ~20–30 min. The activation curves have a sigmoidal shape indicating, that the transition between the inactive to the active state is characterized by a threshold. If the binding between FADD and clustered FAS death

domains is disrupted, then only <0.1% of active caspase-8 is observed upon receipt of the death signal, which is consistent with experimental observations [41].

Fussenegger's model of stress-mediated apoptosis regulation describes the activation of procaspase-9 by cytosolic cytochrome *c*, and the apoptotic protease-activating factor 1 (Apaf-1) complex. Activated caspase-9 then activates apoptosis executioner caspases at some specific rate. Formation of the Apaf-1–cytochrome *c* complex is inhibited by antiapoptotic Bcl-2 family members such as Bcl-xL. Proapoptotic Bcl-2 family members (e.g., Bax, Bak) can bind to antiapoptotic family members and remove their inhibitory effect. The ratio of anti- versus pro-apoptotic Bcl-2 family members is controlled by the p53 transcription factor that is activated in cells under stress conditions. Simulation results of stress-induced caspase activation dynamics were consistent with experimental observations [42]. Specifically, the model shows that cytochrome *c* is released within 10 min after a stress death signal is received, which results in procaspase-9 activation, 35–40% of the executioner caspase being active within 1 h, and 70% of the executioner caspase being active at 2 h. In addition, simulations revealed that the active fraction of both initiator and executioner caspases is reduced in p53 mutant cells as compared to wild-type cells. Overexpression of antiapoptotic Bcl-2 family members is predicted to block the activation of procaspase-9. The model also confirms that the ratio of anti- versus pro-apoptotic Bcl-2 family members determines whether or not executioner caspases will be activated. The model was then used to predict the effects of combined therapies based on simultaneous receptor- and stress-induced caspase activation.

The model developed by Fussenegger et al. was successful in explaining qualitative experimental observations. However, more quantitative data would be required to complete the model calibration. Quantitative information on reaction rates and molecular concentrations is required to perform reliable mathematical simulations of signal transduction in the apoptosis regulatory network. In 2004, Eissing et al. developed a reduced receptor-induced apoptosis, using parameter values from the literature to evaluate the system behavior within a wide range of parameters [21]. The model revealed that caspase activity remains low for a time that is inversely proportional to the stimulus strength, followed by a steep rise in activity when the input exceeds the threshold; caspase activity then ceases at some maximum level. Bifurcation analysis of the model confirmed that the apoptosis regulation system exhibits a bistable behavior. The same year, Bentele et al. developed a data-based model of receptor-induced apoptosis with parameters estimated on the basis of quantitative experimental data [25]. The time series data for concentrations of 15 different molecules after activation of FAS receptors were used to calibrate the core model of the FAS-induced apoptosis. In addition, data from distinct apoptosis activation scenarios in response to different initial values of ligand concentration were used to improve the estimation of model parameters. The model predicted that apoptosis is not executed when a ligand–receptor concentration ratio is below a critical value, which was also confirmed by experimental observations. In conclusion, Bentele et al. proposed a threshold mechanism for induction of receptor-induced apoptosis. A year later, Hua et al. published a FAS-induced apoptosis model to investigate the effects of altering the level of Bcl-2 on the kinetics of caspase-3 activation [43]. The model predicts that Bcl-2 blocks the mitochondrial pathway by binding to proapoptotic Bax, Bak, and tBid proteins. Further, the model predicts that apoptosis signaling flow can be switched between mitochondriadependent and mitochondria-independent pathways by varying molecular component levels without changing network structure.

In 2006, Legewie et al. developed a quantitative kinetic model of intrinsic (stressinduced) apoptosis, which displays an all-or-none behavior of caspase activation in response to an apoptotic stimulus [22]. The model helped to identify the positive feedback mechanism that allows cells to achieve ultrasensitivity and bistability in cell death decision making. The pathway molecular regulators that control the apoptotic threshold stimulus and integrate multiple inputs into an all-or-none caspase output were also determined. Time-course simulation results agreed with experimental observations that the induction

of maximal caspase-3 cleavage after exogenous addition of cytochrome *c* occurs within ~15–60 min. Furthermore, cytochrome *c*-induced activation of caspase-3 was observed to be bistable and irreversible. The bistable and irreversible caspase-3 activation arises in the system due to XIAP-mediated feedback that cooperates with caspase-9 cleavage by caspase-3. X-linked inhibitor of apoptosis (XIAP) inhibits the catalytic activities of caspase-9 and caspase-3 through reversible binding. The feedback cleavage of caspase-9 by caspase-3 leads to autoamplification of the apoptotic signal. Simulation results show that XIAP-mediated feedback is observed only if caspase-9 and caspase-3 compete for binding to XIAP. Depletion and re-addition experiments using different Apaf-1, caspase-3, caspase-9, and/or XIAP concentrations were proposed to test the all-or-none caspase activation.

Also in 2006, Rehm et al. published a computational model of apoptosome-dependent caspase activation based on biochemical data from HeLa cells [26]. The model predicts that the all-or-none apoptotic response depends on caspase-3-dependent feedback signaling and XIAP, which was then verified quantitatively using single-cell experiments with a caspase fluorescence resonance energy transfer substrate. A concentration threshold of XIAP between 0.15 and 0.30 µM, controlling the substrate cleavage by effector caspases, was identified. The model suggested that high levels of XIAP may promote apoptosis resistance and sublethal caspase activation. This result agrees with a computational analysis that was performed earlier, which also suggested that the inhibitor of apoptosis plays an important role in both the induction and prevention of apoptosis [44]. Conversely, Bagci et al. proposed a mathematical model of mitochondria-dependent apoptosis to study both the role of Bax and Bcl-2 synthesis, degradation rates and the number of mitochondrial permeability transition pores involved in the cell response to a death signal [23]. The main finding was that the transition from bistable to monostable (survival) cell behavior is controlled by the synthesis and degradation rates of Bax and Bcl-2 and by the number of mitochondrial permeability transition pores. Also, the model results suggested that cooperative apoptosome formation is a much more robust mechanism to induce bistability than feedback mechanisms involving, for example, the inhibition of caspase-3 by the inhibitor of apoptosis. Later, Chen and Cui et al. analyzed the robustness of Bax and Bcl-2 apoptotic switches using both deterministic and stochastic models [38,45,46]. These mechanisms were confirmed to be bistable and robust to noise and wide ranges of parameter variation.

Albeck et al. developed a mathematical model of extrinsic, receptor-induced apoptosis to explain the molecular mechanism of the variable-delay, snap-action switch function that determines the cell choice between life and death [27]. The model was calibrated by experimental data collected from live-cell imaging, flow cytometry, and immunoblotting of cells perturbed by protein depletion and overexpression. The model was then used to reveal the mechanism by which a steady and gradual increase in caspase-8 activity is converted into a snap-action downstream signal. Permeabilization of the mitochondrial membrane and relocalization of proteins are the key factors in the extrinsic apoptosis network by which a graded signal that activates caspase-8 and promotes the formation of pores in the mitochondrial membrane is transformed into an all-or-none death decision. Importantly, such snap-action behavior at the level of the mitochondrial outer membrane permeabilization occurs independently of caspase-dependent feedback mechanisms. The formation of pores in the mitochondrial membrane involves the pore-forming proteins Bax and Bak that can self-assemble into transmembrane pores, which are antagonized by anti-apoptotic Bcl-2 proteins [47]. Cytochrome *c* is released into the cytosol when the level of active pore-forming proteins exceeds the threshold set by anti-apoptotic Bcl-2 proteins. Using experimental and modeling techniques, Spencer et al. demonstrated that cell-to-cell variability in time-to-death significantly depends on the activation rate of the tBid protein that activates the pore-forming proteins, Bax and Bak [33]. Therefore, in the case of receptor-mediated apoptosis, the timing and probability of death relies on the differences in the protein levels that can be caused, for example, by noise in gene expression. Furthermore, the stochastic protein turnover in a receptor-mediated apoptosis model can result in fractional killing [48].

Later models were developed to investigate crosstalk between apoptosis regulation and NF-κB pathways [32], the estrogen signaling network [31], endoplasmic pathways [28], and autophagy regulation [29]. Neumann et al. described a model of the crosstalk between receptor-mediated apoptosis regulation and NF-κB signaling that are activated by the same receptor in parallel to the apoptotic signaling and on a similar time scale [32]. Model and experimental analysis suggested that the balance between apoptotic and NF-κB signaling is shaped by the proteins that regulate the assembly dynamics of the death-inducing signaling complex (DISC). Therefore, the assembly of DISC acts as a signal processor, determining life/death decisions in a nonlinear manner. Tyson et al. provided a roadmap for a detailed mathematical model that would allow researchers to characterize the crosstalk among the estrogen signaling network, apoptosis, autophagy, and cell cycle regulations in breast epithelial cells [31]. Later, the same research lab published a detailed mathematical model to examine the decision process that moves a cell from autophagy to apoptosis [29]. The model was successful in explaining quantitative time-course data of autophagy and apoptosis under cisplatin treatment. Further, the model allows for characterization of the prosurvival and prodeath cell responses to cytotoxic stress. Also, in 2012, Hong et al. published a model of cisplatin-induced apoptosis that integrates the death receptor pathway, and mitochondrial and endoplasmic reticulum stress response mechanisms [28]. The model predicts the relative contribution of each signaling pathway to apoptosis. Simulation results revealed that the mitochondrial and death receptor pathways as well as crosstalk among pathways make the greatest contribution to the level of apoptosis, whereas the contribution of the endoplasmic reticulum stress pathway is negligible.

### *The Role of p53 in Apoptosis*

The tumor suppressor gene p53 (TP53) has been reported as an upregulated modulator of apoptosis and as a driver of cell fate transition from cell cycle arrest to apoptosis [49]. Mathematical models that characterize the p53 contribution to apoptosis have been developed by several groups [7,23,28,30,50]. p53 targets many genes regulating cell apoptosis, including BCL2 and BAX genes [51]. Computational study of apoptosis regulation shows that the balance between anti- and proapoptotic Bcl-2 family members is altered in p53 mutant cells [7]. Also, the active fraction of both initiator and executioner caspases is reduced in p53 mutant cells as compared with wild-type cells. The mathematical model also predicts that overexpression of the death ligand and the FAS receptor can be used to initiate executioner caspase activation in p53 mutant cells [7]. Bagci et al. have shown that apoptosis is not sensitive to caspase-3 activation when p53 expression is low, and that bistability to apoptotic stimuli is observed when p53 level is high [23]. Predictions from this apoptosis model agree with experimental data [52]. Another study reported that inhibition of p53 protects against cisplatin-induced apoptosis [28]. Cisplatin induces DNA damage that results in the phosphorylation and activation of p53. There, the activation of Bax by p53 induces mitochondrial membrane permeabilization and apoptosis [53]. Also, p53 mediates caspase-2 activation and the mitochondrial release of apoptosis-inducing factor. The model predicts time courses for p53, caspase-2, Bax activation, apoptosis-inducing factor release and apoptosis activation. Simulation results agree with experimental data that p53 inhibition prevents the mitochondrial release of apoptosis-inducing factor and cisplatin-induced apoptosis [54]. Overexpression of p53 results in caspase-2 activation and also the mitochondrial release of apoptosis-inducing factors [54].

Ballweg et al. developed a mathematical model that integrates p53 signaling, cisplatininduced events, and apoptosis regulation that was used to study the dynamics of fractional killing induced by cytotoxic drugs [30]. Many drugs activate not only apoptosis execution signaling but also expression of anti-apoptotic genes, which results only in fractional killing amongst a population of treated cells [55]. Thus, fractional killing may occur due to crosstalk between the apoptosis and survival pathways [56]. The model predicts that the probability of apoptosis depends on the dynamics of p53 and the rate of p53 activation determines the cell fate [30]. Slow activation of p53 results in cell survival, whereas fast p53 activation

induces cell death. This result also agrees with the experimental observation showing that apoptotic cells accumulate p53 much earlier than cells that survive the treatment [55]. In the model, activation of Bax and subsequent execution of apoptosis occur when the level of p53 exceeds a threshold value. However, the apoptosis initiation threshold depends on the inhibitor of apoptosis, cIAP. Cells with an elevated level of cIAP require a higher level of p53 to induce apoptosis. Because the level of apoptosis regulator cIAP increases with time, the rate of p53 activation plays an important role in the determination of cell fate. Cell-to-cell variability due to stochastic gene expression and environmental noise can also set different apoptosis initiation thresholds in different cells, resulting in fractional killing.

Up to this point, we have reviewed mathematical models of apoptosis that use ODEs to describe the mechanism of cell death (apoptosis) regulation. However, other mathematical approaches have been also used to study apoptosis regulation [16,18,34–37]. Several apoptosis models have been developed using a Boolean (logical) approach that can analyze extensive regulatory networks with many molecular components and their interactions [34–36]. Schlatter et al. developed an apoptosis regulation model that comprises 86 nodes and 125 interactions [34]. Mai et al. developed a model that describes 37 internal states of signaling molecules involved in apoptosis regulation, 2 extracellular signal inputs, and the DNA damage event as an output [35]. Calzone et al. developed a model to study crosstalk between receptor-mediated apoptosis regulation, NFκB pro-survival pathways, and RIP1-dependent necroptosis regulation [36]. These models were used to characterize feedback loops in the apoptosis regulation network structure.

While Boolean models are excellent tools to reproduce the *qualitative* behavior of a regulatory network, they are weak at addressing detailed *quantitative* questions about molecular mechanisms [19]. Petri nets have been applied to analyze and validate a qualitative model of extensive apoptosis regulation [16]. Agent-based modeling turned out to be a more appropriate approach for modeling the death-inducing signaling complex assembly than an ODE-based model that must describe a large number of intermediate products involved in DISC assembly [37]. A cellular automata approach has been applied to study apoptosis blocking in the immunological response of T cells by varying the inhibitor actions such as FLIP and IAP [18]. The model predicts that only joint suppression of both FLIP and IAP apoptosis inhibitors can effectively act to kill cancer cells through apoptosis.

In conclusion, comprehensive data and extensive experimental characterization of apoptosis allowed computational and systems biologists to develop several mathematical models of apoptosis regulation. These models not only increase our understanding of mechanisms of apoptosis execution induced by stress or signals, but also predict perturbations that can prevent or enhance apoptosis. An accurate mathematical model of apoptosis can help find novel combinations of existing therapies that can induce the death of cancer cells using low doses. Further studies that integrate apoptosis with other cell death regulations will help to understand the cell death decision mechanism that determines the execution of a specific cell death fate.

### **3. Necroptosis**

Necroptosis is a regulated cell death that can be initiated by changes in extracellular or intracellular homeostasis, detected by specific death receptors [14]. Triggering necroptosis primarily involves the receptor-interacting protein kinase 1 (RIP1), RIP3, and mixed lineage kinase domain-like protein (MLKL). Necroptosis can be induced by the binding of tumor necrosis factor (TNF) or other ligands to cell surface receptors that trigger the sequential phosphorylation of receptor-interacting protein kinases. At a cell physiology level, necroptosis results in cell volume expansion, cell membrane rupture, and intracellular material overflow that leads to a local inflammatory reaction and immune response activation. Necroptosis-inhibiting drugs can be used to treat inflammatory diseases [57]. Necroptosis-promoting drugs are potential anticancer therapies [58]. Studies of necroptosis regulation can help to identify molecular targets that can be used to reprogram the necroptosis execution in a desired direction. While many molecular components involved

in necroptosis regulation are known, the precise interaction network, signaling spread, dynamic behavior of necroptotic regulation, and the decision-making processes within the molecular network, remain poorly understood. Several mathematical models have been developed recently to characterize the dynamics of necroptosis regulation [8,59].

Xu et al. have developed a computational model of the cellular necroptosis signaling network [8], to study the necroptosis signaling dynamics that lead to cell death in the form of oscillation-induced trigger waves. The study focused on the core cellular necroptosis signaling module that includes four components: TRADD, RIP1, caspase-8, and RIP3. The activities of key components are regulated either by phosphorylation, dephosphorylation, or cleavage. The corresponding mathematical model described 4 variables and involved 10 interaction terms. Xu et al. used a Latin hypercube sampling method to randomly scan the model network parameters and evaluate the stable oscillation behavior of the cellular necroptosis signaling circuit. Bifurcation analysis and potential landscape theory were applied to explore oscillation modes in different cellular necroptosis signaling circuits. The results indicate that the cellular necroptosis signaling circuit more likely produces oscillations when the total amount of RIP1 or caspase-8 is high, while fluctuations in the value of RIP3 have no significant effects on the oscillation probability. Also, oscillations are often obtained when the activation of caspase-8 by RIP1 is fast, while RIP3 phosphorylation by RIP1 is relatively slow. Further, oscillations are more robust when the reaction rate constants that describe RIP1 activation by RIP3 are stronger than rate constants describing other interactions. Overall, oscillation robustness analysis revealed three regulatory feedback loops formed by RIP1, caspase-8, and RIP3 interactions. These loops comprise a negative feedback loop: RIP3 activates RIP1, which activates caspase-8, that inhibits RIP3; a positive feedback loop: RIP1 activates RIP3, which inhibits caspase-8, that inhibits RIP1; and an incoherent feedforward loop: RIP1 activates both caspase-8 and RIP3, and caspase-8 inhibits RIP3. Importantly, for oscillations to be robust, the reactions in the positive feedback loop must be slower than reaction rates in the negative feedback loop. Also, a stochastic parameter analysis indicated that the incoherent feedforward loop is the essential molecular mechanism that allows the necroptosis signaling system to generate oscillations.

Xu et al. classified oscillations that occur in cellular necroptosis signaling circuits into four groups according to amplitude and oscillation period. About 50% of observed oscillations had a high-amplitude (above the median value of all the counted amplitudes) and fast period (>100 min based on the oscillation period of NF-κB [60,61]), about 37% of oscillations had a low-amplitude and fast period, ~12% of oscillations had high-amplitude and slow period, and ~1% of oscillations had a slow and low-amplitude period. Further analysis revealed that the inhibition rates of RIP1 and RIP3 by caspase-8 play an important role in determining the amplitude behavior of fast oscillations. In addition, bifurcation analysis revealed that the dynamic behavior of the system can be switched from slow high-amplitude oscillations to slow low-amplitude oscillations by tuning the parameters that describe the activation of caspase-8 by RIP1. However, the transition from fast to slow oscillation behavior cannot be achieved by changing any single reaction rate constant. Also, the system changes dynamics from slow high-amplitude oscillations to fast low- or highamplitude oscillations when two parameters that describe RIP1 inhibition with caspase-8 and RIP1 phosphorylation with RIP3 are simultaneously tuned. Robustness analysis revealed that the period of fast oscillations was more robust to parameter perturbations than the period of slow oscillations. The amplitude of slow low-amplitude oscillations was robust to parameter perturbations, while the robustness of amplitude of fast high-amplitude oscillations was the weakest. Overall, the study provides a quantitative characterization for the mechanism of oscillation mode-switching behavior in the necroptosis signaling network. Xu et al. proposed that MLKL can decode the information according to the amplitude and period of RIP3, which can be an important mechanism that allows cells to generate different responses in various stressful conditions.

A more recent detailed computational model of tumor necrosis factor (TNF)-induced necroptosis has been developed by Ildefonso et al. [59]. The model was derived from

the literature-curated molecular mechanisms of necroptosis regulation, which involves 14 proteins, 37 biochemical species, and 40 reactions. The simplified molecular mechanism that shows key species involved in necroptosis execution is shown in Figure 2. Dynamics of species were described by a set of ordinary differential equations where all reactions were described by the mass action law. The model was calibrated and validated using experimental protein time-course data from a well-established necroptosis-executing cell line. Simulations then confirmed that the model is successful in explaining the dynamics of necroptosis reporter, phosphorylated mixed lineage kinase domain-like protein (pMLKL). Furthermore, four distinct necroptosis execution modes were identified by using a dynamical systems analysis and a spectral clustering algorithm. While the temporal dynamics of pMLKL were similar in each mode of necroptosis execution, the sequences of molecular events that led to MLKL phosphorylation and subsequent necroptotic cell death were different. The modes primarily differed in the values of rate constants across the necroptosis execution pathway. For example, the rate constant for binding of A20 to ubiquitinated RIP1 was significantly smaller in mode 4 than in the other modes, and also smaller in mode 2 relative to modes 1 and 3. Also, mode 4 has a significantly larger activation rate and smaller deactivation rate constant for caspase-8 in complex II. The activation/deactivation of caspase-8 in complex II is a critical step in the pathway for determining whether the cell will progress to necroptosis. Differences in rate constant values create the difference in the action of A20 and CYLD enzymes across four modes that are then able to effectively operate as an inhibitor or activator of necroptosis. necroptosis has been developed by Ildefonso et al. [59]. The model was derived from the literature-curated molecular mechanisms of necroptosis regulation, which involves 14 proteins, 37 biochemical species, and 40 reactions. The simplified molecular mechanism that shows key species involved in necroptosis execution is shown in Figure 2. Dynamics of species were described by a set of ordinary differential equations where all reactions were described by the mass action law. The model was calibrated and validated using experimental protein time-course data from a well-established necroptosis-executing cell line. Simulations then confirmed that the model is successful in explaining the dynamics of necroptosis reporter, phosphorylated mixed lineage kinase domain-like protein (pMLKL). Furthermore, four distinct necroptosis execution modes were identified by using a dynamical systems analysis and a spectral clustering algorithm. While the temporal dynamics of pMLKL were similar in each mode of necroptosis execution, the sequences of molecular events that led to MLKL phosphorylation and subsequent necroptotic cell death were different. The modes primarily differed in the values of rate constants across the necroptosis execution pathway. For example, the rate constant for binding of A20 to ubiquitinated RIP1 was significantly smaller in mode 4 than in the other modes, and also smaller in mode 2 relative to modes 1 and 3. Also, mode 4 has a significantly larger activation rate and smaller deactivation rate constant for caspase-8 in complex II. The activation/deactivation of caspase-8 in complex II is a critical step in the pathway for determining whether the cell will progress to necroptosis. Differences in rate constant values create the difference in the action of A20 and CYLD enzymes across four modes that are then able to effectively operate as an inhibitor or activator of necroptosis.

high-amplitude oscillations was the weakest. Overall, the study provides a quantitative characterization for the mechanism of oscillation mode-switching behavior in the necroptosis signaling network. Xu et al. proposed that MLKL can decode the information according to the amplitude and period of RIP3, which can be an important mechanism that al-

A more recent detailed computational model of tumor necrosis factor (TNF)-induced

lows cells to generate different responses in various stressful conditions.

*Entropy* **2022**, *24*, x FOR PEER REVIEW 9 of 23

**Figure 2.** The influence diagram shows key molecular components involved in necroptosis regulation. Lines represent activation (arrowhead) and inhibition (bar head) influences.

Taken together, the computational analysis helped to resolve the controversy in experimental observations by showing that CYLD- and A20-driven deubiquitination of RIP1 may act as pro- and anti-necroptotic in different cell types. According to Ildefonso et al.'s model, knocking out A20 decreases the probability of necroptosis execution (necroptosis sensitivity) in mode 1, and increases the sensitivity to necroptosis in mode 2 [59]. Con-

versely, knocking out CYLD increases the sensitivity to necroptosis in modes 1 and 4, and decreases the sensitivity to necroptosis in mode 2. Knocking out CYLD or A20 has no effect in mode 3. Also, A20 knockout has no effect in mode 4. These results have been compared to cell phenotype observations in A20 and CYLD knockdown experiments in different cell types. For example, it has been reported that RIP1 is deubiquitinated by both A20 and CYLD in mouse fibrosarcoma cells, but inhibition of CYLD protects cells from necroptosis, while A20 depletion can sensitize cells to death by necroptosis [62]. Thus, A20 and CYLD depletion experiments in mouse fibrosarcoma cells are consistent with the model results obtained for A20 and CYLD knockouts in mode 2.

TNF, TNFR, and MLKL are three common protein modulators of necroptosis across the four modes of necroptosis execution. Furthermore, rate constants that control the association of TNF to TNFR, ubiquitination of RIP1 by cIAP in complex I, and association of LUBAC to complex I can be used to efficiently modulate necroptosis execution across the four modes. Therefore, targeting these modulators can be used to induce or prevent necroptosis, potentially useful for both cancer therapy and treatment of inflammatory diseases.

Apoptosis and necroptosis regulation networks share common nodes and edges and may suppress each other [63]. Either apoptosis or necroptosis can be induced by TNF and the cell death decision depends on the cell state. Complex II can recruit RIP3 to form a necrosome or recruit caspase-8 to stabilize its active conformation, resulting in the release of an activated caspase 8 homodimer that then can induce apoptosis [64]. Li et al. [65] performed a quantitative study of crosstalk between the apoptosis and necroptosis pathways. Specifically, mathematical modeling was used to investigate three possible mechanisms of caspase-8 activation by (i) TRADD, (ii) RIP1, and (iii) TRADD and RIP1 together. The law of mass action was used to convert the proposed molecular mechanisms into a system of ODEs. Simulations of each mechanism were compared with data obtained using the sequential window acquisition of all theoretical fragment ion spectra mass spectrometry methods. All three mechanisms reproduced the amounts of major components in TNFR1, RIP1, and RIP3 complexes. However, only mechanism (ii) could explain a negative regulation of RIP3 phosphorylation by the increase in RIP1 levels. This result was also supported by a sensitivity analysis showing that the most robust negative regulation of RIP3 phosphorylation by RIP1 is achieved when mechanism (ii) is used in the model. To test this prediction, Li et al. experimentally knocked down RIP1 to three different expression levels by using RIP1-specific short hairpin RNA and measured the increase in TNF-induced phosphorylation of RIP3 and MLKL. Deletion of RIP1 completely blocks TNF-induced RIP3 phosphorylation [65]. In addition, simulation results show that pro-caspase-8 activity is necessary for the up-regulation of RIP3 phosphorylation by decreasing RIP1 expression. The mechanism was further refined to make it in agreement with the observation that TNF induces quick caspase-8 activation and apoptosis in RIP1 KO cells [62]. Specifically, TRADD-dependent caspase-8 activation was added to the mechanism (ii). The final model successfully explained both RIP10 s biphasic roles in necroptosis, where RIP1 promotes necroptosis within an extremely low-level range (<∼2% of wildtype) and inhibits necroptosis at higher levels, and the activation level of caspase-8 in RIP1 KO cells. Also, the response of pro-caspase-8 to RIP1 level is linear, whereas RIP3 phosphorylation is determined by the nonlinear (ultrasensitive) threshold pattern.

Overall, a quantitative approach has been applied successfully to describe the roles of RIP1 in cell death determination. In conclusion, Li et al. proposed a "speed competition" decision mechanism in which cells decide to execute apoptosis or necroptosis by the pathway that reaches the final destination first. Interestingly, simultaneous execution of necroptosis and apoptosis has been observed in some individual cells [65].

### **4. Pyroptosis**

The regulated cell death that is associated with the formation of plasma membrane pores by members of the gasdermin protein family is called pyroptosis [14]. The induction of pyroptosis may occur as a consequence of inflammatory caspase activation that can be triggered by pathogen invasion such as Gram-negative bacteria. The critical role of caspase-driven pyroptosis for innate immune responses against invading bacteria has been confirmed in experiments with mice carrying gene mutations that disrupt normal activity of caspase proteins [66]. By killing the host cell, pyroptosis removes the replication compartment of intracellular pathogens and thus prevents their spreading. Hence, pyroptosis has an important role in innate immunity against intracellular pathogens. of pyroptosis may occur as a consequence of inflammatory caspase activation that can be triggered by pathogen invasion such as Gram-negative bacteria. The critical role of caspase-driven pyroptosis for innate immune responses against invading bacteria has been confirmed in experiments with mice carrying gene mutations that disrupt normal activity of caspase proteins [66]. By killing the host cell, pyroptosis removes the replication compartment of intracellular pathogens and thus prevents their spreading. Hence, pyroptosis has an important role in innate immunity against intracellular pathogens.

Overall, a quantitative approach has been applied successfully to describe the roles of RIP1 in cell death determination. In conclusion, Li et al. proposed a "speed competition" decision mechanism in which cells decide to execute apoptosis or necroptosis by the pathway that reaches the final destination first. Interestingly, simultaneous execution of

The regulated cell death that is associated with the formation of plasma membrane pores by members of the gasdermin protein family is called pyroptosis [14]. The induction

*Entropy* **2022**, *24*, x FOR PEER REVIEW 11 of 23

**4. Pyroptosis** 

necroptosis and apoptosis has been observed in some individual cells [65].

Pyroptosis induced by inflammatory caspases is driven by the gasdermin protein GSDMD. Caspases activate GSDMD that then translocates to the plasma membrane where GSDMD induces pore formation and thus rapid plasma membrane permeabilization. The simplified molecular mechanism of the pyroptosis induced by inflammatory caspases is shown in Figure 3. In this scheme, pyroptosis relies on caspase-1 activation. Pyroptosis induced by inflammatory caspases is driven by the gasdermin protein GSDMD. Caspases activate GSDMD that then translocates to the plasma membrane where GSDMD induces pore formation and thus rapid plasma membrane permeabilization. The simplified molecular mechanism of the pyroptosis induced by inflammatory caspases is shown in Figure 3. In this scheme, pyroptosis relies on caspase-1 activation.

**Figure 3.** The influence diagram shows key molecular components involved in GSDMD- and GSDME-executed death modes called pyroptosis and secondary pyroptosis respectively. Lines represent activation (arrowhead) and inhibition (bar head) influences. **Figure 3.** The influence diagram shows key molecular components involved in GSDMD- and GSDMEexecuted death modes called pyroptosis and secondary pyroptosis respectively. Lines represent activation (arrowhead) and inhibition (bar head) influences.

Beyond inflammatory settings, pyroptotic cell death can be induced by TNF, various DNA-damaging agents, and infection with vesicular stomatitis virus [67,68]. In these cases, pyroptosis is driven by other members of the gasdermin family, specifically GSDME. This form of pyroptosis releases fewer inflammatory cytokines than is observed when pyroptosis is induced by inflammatory caspases. Pyroptotic signaling relies on the activation of caspase-3 that catalyzes proteolytic cleavage of GSDME. The identification of other gasdermin family members that execute pyroptosis in conditions that are beyond inflammatory settings has been significantly expanded [14]. Beyond inflammatory settings, pyroptotic cell death can be induced by TNF, various DNA-damaging agents, and infection with vesicular stomatitis virus [67,68]. In these cases, pyroptosis is driven by other members of the gasdermin family, specifically GSDME. This form of pyroptosis releases fewer inflammatory cytokines than is observed when pyroptosis is induced by inflammatory caspases. Pyroptotic signaling relies on the activation of caspase-3 that catalyzes proteolytic cleavage of GSDME. The identification of other gasdermin family members that execute pyroptosis in conditions that are beyond inflammatory settings has been significantly expanded [14].

A computational study of the crosstalk between caspase-1- and caspase-3-driven pyroptosis pathways was performed by Zhu et al. [9]. The molecular regulatory network that executes pyroptosis via activation of GSDMD and GSDME is shown in Figure 3. The crosstalk between caspase-1- and caspase-3-driven pyroptosis pathways is realized through tBid, caspase-9, and caspase-8 components. Zhu et al. developed a mathematical model that describes the dynamics of seven molecular components and the dynamics of the cell population governed by cell proliferation and death processes. The model consists of eight coupled ODEs and 83 parameters. Hill functions were used to describe activation A computational study of the crosstalk between caspase-1- and caspase-3-driven pyroptosis pathways was performed by Zhu et al. [9]. The molecular regulatory network that executes pyroptosis via activation of GSDMD and GSDME is shown in Figure 3. The crosstalk between caspase-1- and caspase-3-driven pyroptosis pathways is realized through tBid, caspase-9, and caspase-8 components. Zhu et al. developed a mathematical model that describes the dynamics of seven molecular components and the dynamics of the cell population governed by cell proliferation and death processes. The model consists of eight coupled ODEs and 83 parameters. Hill functions were used to describe activation and inactivation reactions for molecular components. The values of 44 parameters were estimated from sources available in the literature and 39 parameters were estimated using 138 time-course data points that were measured for eight variables (the death rate and seven molecular components) in wild-type cells and cells with single, double, and triple knockouts of the molecular components.

The pyroptosis decision mechanism was analyzed using bifurcation and sensitivity analysis methods. Bifurcation analysis revealed that the change in expression levels of caspase-1, caspase-3, and GSDMD can switch between GSDMD- and GSDME-executed pyroptosis death modes. Furthermore, the transition between pyroptosis death modes could not be efficiently controlled by varying the expression levels of caspase-8, caspase-9,

tBid, or GSDME. According to the model, GSDMD-driven pyroptosis is more likely when the caspase-1 total expression level is below ∼1.5 nM and GSDME-driven pyroptosis occurs when the caspase-1 level is above 14 nM. For caspase-1 levels ranging from 1.5–14 nm, bistability is observed when either GSDMD- or GSDME-driven pyroptosis may occur. Similarly, when GSDMD level is lower than 88 nM, GSDME-driven pyroptosis is induced, whereas cells can selectively execute either pyroptosis mode when the level of GSDMD is between 88 nM and 165 nM. GSDMD-driven pyroptosis occurs when GSDMD level is higher than 165 nM. Also, cells execute GSDMD-driven pyroptosis when caspase-3 level is lower than 250 nM, and selectively induce either GSDMD- or GSDME-executed pyroptosis with higher levels of caspase-3.

Sensitivity analysis confirmed that the expression levels of GSDMD and caspase-1 can efficiently change the pyroptosis death modes. This result agrees with experimental observations [69,70]. In addition, bifurcation analysis predicts that the expression level of caspase-3 can also change the pyroptosis death mode between caspase-1- and caspase-3 driven pyroptosis. Overall, the model predicted 3 molecular components and 12 reactions that can be targeted to control the switch between modes of pyroptosis execution. Drugs that can switch between pyroptosis death modes can help to improve treatment protocols for cancer and inflammasome-mediated diseases. For example, GSDME-induced pyroptosis can act as a tumor suppressor [71,72] and also releases fewer inflammatory cytokines when compared to pyroptosis that is executed by GSDMD.

Li et al. extended the GSDMD-induced pyroptosis model by adding apoptosis regulation [73]. The model allows one to study the crosstalk between pyroptosis and apoptosis and inflammasome-induced cell death under different perturbation conditions. Simulation results reproduce the dynamics of cell death executioners in multiple knockout cells. Pyroptosis and apoptosis events are determined by the level of cleaved GSDMD and cleaved caspase-3, respectively. Sensitivity analysis was performed to determine the molecular components that can significantly affect the occurrence of pyroptosis and apoptosis. The model predicted that caspase-1 and GSDMD are key molecular regulators directing the signal flow that can switch cell death modes between pyroptosis and apoptosis. Decreases in caspase-1 or GSDMD gradually inhibit pyroptosis and enhance apoptosis induction. These model predictions were validated by caspase-1 and GSDMD-knocked down experiments. Furthermore, the model results helped to suggest the death signal propagation pathways, resulting in pyroptosis or apoptosis in cells expressing different levels of caspase-1 or GSDMD. To understand the roles of caspase-1 and GSDMD in triggering the cell death modes, Li et al. employed a potential landscape approach. The cell death landscape was represented by potential wells corresponding to pyroptosis and apoptosis death modes. In the double-well potential landscape, the system evolved into one of the two wells from any initial condition. Caspase-1 or GSDMD could change the potential landscape from monostable to bistable. A monostable landscape corresponding to pyroptosis is obtained in cells with a high expression level of caspase-1 or GSDMD; the potential landscape changes to bistable and then to an apoptotic monostable as the expression level of caspase-1 or GSDMD decreases. Overall, the model helps to understand the inflammasome-induced cell death, crosstalk between pyroptosis and apoptosis, and may be used to determine potential molecular targets for driving cells into a desired death execution mode.

### **5. Ferroptosis**

Ferroptosis is another regulated cell death mechanism that involves iron-catalyzed lipid damage [14,74,75]. Cell death occurring by ferroptosis correlates with the accumulation of markers of lipid peroxidation and can be suppressed by iron chelators, inhibitors of lipid peroxidation, and lipophilic antioxidants [75]. Ferroptotic cell death can be modulated pharmacologically and genetically by perturbing lipid repair systems that involve glutathione and glutathione peroxidase 4 (GPX4) that convert toxic lipid hydroperoxides (L-OOH) into non-toxic lipid alcohols (L-OH) [76]. Depletion or inactivation of GPX4 results in overwhelming lipid peroxidation that causes cell death. Ferroptosis also depends on a set of

enzymatic reactions that regulate the biosynthesis of membrane polyunsaturated fatty acids (PUFA)-containing phospholipids, which are the substrates of pro-ferroptotic lipid peroxidation products [75]. Also, the formation of coenzyme-A-derivatives of PUFAs (PUFA-CoA) and their insertion into phospholipids are necessary for the induction of a ferroptotic death signal. Two enzymes, ACSL4 and LPCAT3 are involved in the biosynthesis and remodeling of PUFAs [75,77]. Depletion of PUFAs can suppress the occurrence of ferroptosis, and loss of ACSL4 and LPCAT3 gene products increases resistance to ferroptosis [75].

Iron induces the accumulation of lipid peroxides and thus is important for the execution of ferroptosis. Intracellular iron levels depend on the iron efflux pump ferroportin and the iron importer TFR1 and other proteins that regulate iron import, export, and storage [78–80]. Also, for ferroptosis to start, phospholipid molecules containing polyunsaturated fatty acids (LH-P) are formed from PUFA-CoA, which are then oxidized into lipid hydroperoxides (L-OOH) and eventually into lipid radicals (LO\*). LH-P generation is regulated by LPCAT3 and inhibited by monounsaturated fatty acids (MUFAs). Production of MUFAs depends on desaturation of the saturated fatty acids (SFAs) which is catalyzed by the desaturase SCD1 [81]. Formation of lipid radicals LO\* is promoted by reactive oxygen species (ROSs) and lipid peroxidation enzymes including ALOX15 [74,82]. The generation of endogenous lipid radicals initiates ferroptosis. In addition, ferroptotic cell-death responses can be modulated by p53 activity [83]. For example, induction of SAT1, a transcription target of p53, is correlated with the expression levels of ALOX15 [83]. The influence diagram that reflects the molecular mechanism of ferroptosis is shown in Figure 4. Overall, ferroptosis is morphologically and mechanistically different from all other types of regulated cell death. Regulation of ferroptosis has great potential for cancer therapy, and molecular targets that promote ferroptosis are being actively explored [84]. *Entropy* **2022**, *24*, x FOR PEER REVIEW 14 of 23

**Figure 4.** The influence diagram shows key molecular components involved in ferroptosis regulation. Solid arrows represent reactions of transformation, activation (arrowhead), and inhibition (bar head). Dashed lines with arrowheads or bar heads represent activation or inhibition of a reaction, respectively. **Figure 4.** The influence diagram shows key molecular components involved in ferroptosis regulation. Solid arrows represent reactions of transformation, activation (arrowhead), and inhibition (bar head). Dashed lines with arrowheads or bar heads represent activation or inhibition of a reaction, respectively.

Konstorum et al. developed a stochastic, multistate, discrete mathematical model of ferroptosis regulation [10]. The model describes states of eleven variables that represent

Konstorum et al. used a system-level analysis to study how different input conditions and parameters alter ferroptosis sensitivity. They found that ferroptosis sensitivity depends on PUFA synthesis and PUFA incorporation into the phospholipid membrane, as well as the balance between levels of pro-oxidant species (ROS, lipoxygenases) and antioxidant factors (GPX4). Ferroptosis sensitivity can be reduced by altering parameters that minimize the production of L-OOH species. High ACSL4 and low SCD1 levels result in high ferroptosis sensitivity. The model also predicted that a high level of SCD1 can inhibit ferroptotic induction even when levels of ACSL4 are high. These model predictions were validated using an in vitro experimental system of an ovarian cancer stem cell culture [10]. Overall, the model allows us to better understand the crosstalk between pathways trans-

Immunogenic cell death (ICD) is a regulated cell death mechanism that induces an immune response in the hosts [14]. Basically, ICD is an immunostimulatory form of apoptosis that is characterized by the ability of dying cells to generate robust adaptive immune

mitting signals from different inputs that induce the execution of ferroptosis.

SLC7A11. Each variable can take on three values that respectively represent low, medium, and high molecular species activity or expression level. Variables are updated using updating rules and an asynchronous update scheme at each discrete time step. Five external inputs representing ACSL4, ferroportin, p53, SCD1, and TFRC, which do not change during the course of the simulation, are used to study the sensitivity of ferroptosis induction

to different signaling and perturbation conditions.

**6. Immunogenic Cell Death** 

Konstorum et al. developed a stochastic, multistate, discrete mathematical model of ferroptosis regulation [10]. The model describes states of eleven variables that represent ALOX15, GPX4, L-HP, LIP, LO\*, L-OOH, LPCAT3, MUFA, PUFA-CoA, ROS, and SLC7A11. Each variable can take on three values that respectively represent low, medium, and high molecular species activity or expression level. Variables are updated using updating rules and an asynchronous update scheme at each discrete time step. Five external inputs representing ACSL4, ferroportin, p53, SCD1, and TFRC, which do not change during the course of the simulation, are used to study the sensitivity of ferroptosis induction to different signaling and perturbation conditions.

Konstorum et al. used a system-level analysis to study how different input conditions and parameters alter ferroptosis sensitivity. They found that ferroptosis sensitivity depends on PUFA synthesis and PUFA incorporation into the phospholipid membrane, as well as the balance between levels of pro-oxidant species (ROS, lipoxygenases) and antioxidant factors (GPX4). Ferroptosis sensitivity can be reduced by altering parameters that minimize the production of L-OOH species. High ACSL4 and low SCD1 levels result in high ferroptosis sensitivity. The model also predicted that a high level of SCD1 can inhibit ferroptotic induction even when levels of ACSL4 are high. These model predictions were validated using an in vitro experimental system of an ovarian cancer stem cell culture [10]. Overall, the model allows us to better understand the crosstalk between pathways transmitting signals from different inputs that induce the execution of ferroptosis.

### **6. Immunogenic Cell Death**

Immunogenic cell death (ICD) is a regulated cell death mechanism that induces an immune response in the hosts [14]. Basically, ICD is an immunostimulatory form of apoptosis that is characterized by the ability of dying cells to generate robust adaptive immune responses [85]. The immune response is promoted by damage-associated molecular patterns (DAMPs), which are released by dying cells [86]. DAMPS communicate a state of danger to the organism by activating pattern recognition receptors (PRRs) that are present on the surface of innate immune cells such as monocytes, macrophages, and dendritic cells (DCs) [87]. Activated macrophages and dendritic cells can migrate to the lymph node and pass the antigens to CD8<sup>+</sup> and CD4<sup>+</sup> T lymphocytes, which results in an adaptive immune response. Tumor cell systems are often used to study ICD regulation and dynamics [88]. The immune responses against cancer- or pathogen-driven antigens that induce ICD are well characterized [85]. Importantly, over the past years, developments of ICD-related cancer immunotherapy approaches are gaining great momentum [89].

To study the ICD dynamics of cancer cells, Checcoli et al. developed a mathematical model that integrates intracellular mechanisms involved in ICD and intercellular interactions among cancer cells, DCs, CD8<sup>+</sup> , and CD4<sup>+</sup> T cells [11]. The modeling approach is based on a continuous time Boolean Kinetic Monte-Carlo formalism that was successfully applied to model different complex molecular mechanisms [90]. The aim of the mathematical characterization of ICD processes was to identify the regulatory molecular targets and combinations of pharmacological compounds that can increase anticancer immunity. The model can predict the time-dependent size of different cell populations involved in ICD that is induced by a treatment exposure.

To determine the role of each of the main cell types involved in ICD, Checcoli et al. first simulated a core ICD mechanism that is merely sufficient to reproduce ICD events observed experimentally [11]. The core regulatory mechanism describes the release of CALR, ATP, and HMGB1 molecules from dying tumor cells, and inner-state activation or evolution of immature DC, activated DC, migrating DC, lymph node DC, T cell, and cytotoxic T lymphocyte cell types. As shown in Figure 5, also included are two processes: tumor cell division, which is inhibited by T cells, and death, which is promoted by cytotoxic T lymphocytes. The states of molecules and cells are described by Boolean variables that can take two values: **1** for active or present and **0** for inactive or absent. The system state is described by a vector of Boolean values that represent each molecule, process, and cell type

in the system. In the probabilistic description, the probability is assigned to each system state, such that the sum of probabilities over all possible system states is equal to **1**. Then, to determine the number of cells in a given system state, the system state probability is multiplied by the overall size of the cell population. *Entropy* **2022**, *24*, x FOR PEER REVIEW 16 of 23

**Figure 5.** The core mechanism of immunogenic cell death. Lines represent activation (arrowhead) and inhibition (bar head) influences. A-DC, M-DC, and L-DC represent activated DC, migrating DC, and lymph node DC cell types, respectively. **Figure 5.** The core mechanism of immunogenic cell death. Lines represent activation (arrowhead) and inhibition (bar head) influences. A-DC, M-DC, and L-DC represent activated DC, migrating DC, and lymph node DC cell types, respectively.

The core model can reproduce the series of events following an ICD-inducing intervention. The release of CALR, ATP, and HMGB1 molecules by dying cancer cells occurs within hours, a slow increase in T cells begins after 100 h, which peaks at 200 h, and the tumor cells are eliminated in about 220 h when a rapid increase in cytotoxic T lymphocyte cell population begins. When the clonal expansion of the cytotoxic T lymphocytes was blocked in the model, tumor cell clearance became less efficient and depended mostly on the direct cytotoxicity of the treatment. The core model can reproduce the series of events following an ICD-inducing intervention. The release of CALR, ATP, and HMGB1 molecules by dying cancer cells occurs within hours, a slow increase in T cells begins after 100 h, which peaks at 200 h, and the tumor cells are eliminated in about 220 h when a rapid increase in cytotoxic T lymphocyte cell population begins. When the clonal expansion of the cytotoxic T lymphocytes was blocked in the model, tumor cell clearance became less efficient and depended mostly on the direct cytotoxicity of the treatment.

To improve the predictive power of the model, Checcoli et al. extended their core model by including more cell types and molecular components as well as the ligand–receptor dynamics that determines intercellular communication. The extended model describes 57 entities and provides more detailed representations of the series of events that were explored by the core model. Simulation results of the extended model also reproduce the succession of events resulting in ICD. Simulations were performed starting with 80% of tumor cells, 10% of dendritic cells, and 5% of inactive CD4+ and CD8+ cells. The population of tumor cells rapidly decays starting from 250 h when cytotoxic T lymphocytes are engaged to eliminate tumor cells. To improve the predictive power of the model, Checcoli et al. extended their core model by including more cell types and molecular components as well as the ligand– receptor dynamics that determines intercellular communication. The extended model describes 57 entities and provides more detailed representations of the series of events that were explored by the core model. Simulation results of the extended model also reproduce the succession of events resulting in ICD. Simulations were performed starting with 80% of tumor cells, 10% of dendritic cells, and 5% of inactive CD4<sup>+</sup> and CD8<sup>+</sup> cells. The population of tumor cells rapidly decays starting from 250 h when cytotoxic T lymphocytes are engaged to eliminate tumor cells.

To assess the extended model robustness to parameter changes, Checcoli et al. performed a sensitivity analysis measuring the variations in sizes of tumor cell populations within the 220 h and 280 h time frame when the tumor cell population decreases in the To assess the extended model robustness to parameter changes, Checcoli et al. performed a sensitivity analysis measuring the variations in sizes of tumor cell populations within the 220 h and 280 h time frame when the tumor cell population decreases in the

standard conditions (WT) of the extended model [11]. The decrease in size of the tumor cell population was seen to be delayed only for a few parameter changes when compared with the WT condition. Changes in parameters that control the number of DCs gave the

amount enhanced the death process. Changes in parameters controlling the rate of T cell clonal expansion give a similar effect on the cell death process. Sensitivity analysis also

standard conditions (WT) of the extended model [11]. The decrease in size of the tumor cell population was seen to be delayed only for a few parameter changes when compared with the WT condition. Changes in parameters that control the number of DCs gave the strongest effect. A lower amount of DCs delayed the time of death, whereas a higher amount enhanced the death process. Changes in parameters controlling the rate of T cell clonal expansion give a similar effect on the cell death process. Sensitivity analysis also suggested the points of intervention that had the strongest effect on ICD. For example, a complete knockout of CD28 or CD80 (costimulatory molecules for T cell activation) resulted in a failure of the ICD-inducing treatment (80% of the tumor cell population persists at t = 280 h). By contrast, an external treatment that increases Interleukin-2 (IL-2) could kill the tumor cells faster, at t = 200 h.

The Boolean approach does not provide quantitative details and different regimens of drug treatments. Nevertheless, the model characterizes ICD events and dynamics in cancer cells and predicts molecular targets that could increase tumor clearance. For future directions, Checcoli et al. suggested to include specific in vitro and in vivo experiments to identify parameter values that will agree with experimentally observed timing of the different events leading to tumor clearance [11]. Further extension of the model including effects of IFNγ or TGFβ on the immune cells, and major signaling pathways inside each cell type, will allow the model to predict more feasible pharmacological interventions that can boost ICD for killing tumor cells.

### **7. Discussion**

The significant progress that has been made in the mathematical characterization of different cell death execution pathways offers quantitative insight into cell death control and mechanistically explains why and how a living cell may die. Table 3 summarizes cell death mathematical model development over a 22-year period. We include the modeled cell death mechanism, methods, a mathematical description of the cell death event used in each model, and the main modeling results obtained in each work. ODE and Boolean logic-based approaches are the most common mathematical techniques used to model cell death mechanisms. However, a physical description based on the potential landscape theory has been recently applied to study stochastic dynamics and global stability of cell death signaling pathways [8,73]. In this approach, the steady state probability distribution of a system *Pss* and a dimensionless potential function *E* are related via Boltzmann relation: *E* = −ln(*Pss*) [91]. The physical description allows one to employ thermodynamics to analyze cell death regulatory circuits. Conversely, entropy-based approaches have been applied to analyze biological networks [92] and a cell fate selection process [93,94], they have not been yet applied to characterize cell death decision mechanisms. Therefore, one promising future direction is to describe cell death networks using physical approaches that could help to reveal new functional system states and unknown properties of cell death regulatory mechanisms.


**Table 1.** Summary of cell death mechanism models.

*Entropy* **2022**, *24*, 1402


### **Table 2.** *Cont.*

*Entropy* **2022**, *24*, 1402


**Table 3.** *Cont.*

Importantly, many different cell death pathways share common molecular components, and thus all these pathways can interact together at any time to form a complex mechanism. Therefore, we hypothesize that cell death can be controlled by a singular, highly integrated cell death decision network, see Figure 6. This network enables cells to alter the signal flow through the shared nodes but with different edges and so select alternative cell death execution pathways within a single control network of cell death. A stress death signal can thus initiate multiple death mechanisms but not all reach an execution threshold. Currently, the molecular mechanism that regulates the selection of each specific death execution pathway remains elusive. In addition, mathematical models developed to study crosstalk between necroptosis and apoptosis [65], pyroptosis and apoptosis [73], autophagy and apoptosis [29] support the hypothesis that signals propagating through different cell death pathways are integrated to process the execution of specific cell death. We are developing a mathematical model of the cell death decision network to predict the molecular species and interactions that direct the signal flow towards a specific irreversible cell death fate. Such a model will provide new insights into the integrated control of cell death. Model predictions will help develop new approaches to either block or initiate irreversible cell death and identify which cell death pathways are blocked and which pathways remain accessible to execute cell death. Thus, model predictions will

Li et al., 2022 [73], crosstalk between pyroptosis and apoptosis regulations

Konstorum et al., 2020 [10], ferroptosis regulation

Checcoli et al., 2020 [11], immunogenic cell death (ICD) mechanism

*Entropy* **2022**, *24*, x FOR PEER REVIEW 19 of 23

between GSDMD- and GSDME-executed pyroptosis death modes

Caspase-1 and GSDMD are key proteins that regulate the switching between pyroptosis and apoptosis

Ferroptosis sensitivity depends on PUFA synthesis, PUFA incorporation into the phospholipid membrane, and the balance between levels of pro-oxidant species and antioxidant factors

The succession of events resulting in ICD. Points of intervention that had the strongest effect on ICD

ratio of dying cell population to the initial cell population

ODE and potential energy landscape approaches. DR: by levels of cleaved GSDMD (pyroptosis) and cleaved caspase-3 (apoptosis)

Stochastic, multistate, discrete mathematical approach. DR: intermediate and high levels of the lipid radical LO\*

Boolean Kinetic Monte-Carlo approach. DR: Death node is at **1**

suggest alternative interventions to overcome a block in cell death activation that can occur in cancer cells that acquire drug resistance. which pathways remain accessible to execute cell death. Thus, model predictions will suggest alternative interventions to overcome a block in cell death activation that can occur in cancer cells that acquire drug resistance.

Importantly, many different cell death pathways share common molecular components, and thus all these pathways can interact together at any time to form a complex mechanism. Therefore, we hypothesize that cell death can be controlled by a singular, highly integrated cell death decision network, see Figure 6. This network enables cells to alter the signal flow through the shared nodes but with different edges and so select alternative cell death execution pathways within a single control network of cell death. A stress death signal can thus initiate multiple death mechanisms but not all reach an execution threshold. Currently, the molecular mechanism that regulates the selection of each specific death execution pathway remains elusive. In addition, mathematical models developed to study crosstalk between necroptosis and apoptosis [65], pyroptosis and apoptosis [73], autophagy and apoptosis [29] support the hypothesis that signals propagating through different cell death pathways are integrated to process the execution of specific cell death. We are developing a mathematical model of the cell death decision network to predict the molecular species and interactions that direct the signal flow towards a specific irreversible cell death fate. Such a model will provide new insights into the integrated control of cell death. Model predictions will help develop new approaches to either block or initiate irreversible cell death and identify which cell death pathways are blocked and

**Figure 6.** The concept of cell death decision network that controls the irreversible execution of different cell death mechanisms. **Figure 6.** The concept of cell death decision network that controls the irreversible execution of different cell death mechanisms.

**Author Contributions:** Conceptualization, P.K. and R.C.; writing—original draft preparation, M.S. and P.K.; writing—review and editing, M.S., P.K. and R.C; visualization, M.S. and P.K.; supervision, P.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** US department of the Army, Breast Cancer Research Program W81XWH-18-1-0722; BC171885, Public Health Service NIH U01CA184902.

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

**Data Availability Statement:** Not applicable.

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

### **References**

