Multi-objective parameter estimation on cultivation of yeast Kluyveromyces marxianus var. lactis MC5

Abstract

In this study, the most appropriate form of the specific growth rate of the yeast Kluyweromyces marxianus var lactis MC5 was determined by a multi-criteria decision-making method (PROMETHEE II). Then, based on experimental data from a total of six batch processes, the simultaneous estimation of the kinetic parameters for four batch processes was made. The remaining two experiments were used to test the model. The parameter estimation problem was designed as a multi-objective optimisation problem. The method of scalarisation of the weighted sum was used. In this way, the multi-objective optimisation problem was converted to a single-objective one. The direct search method was used to solve the problem of obtaining the Pareto solution. In order to test the resulting solution, two batches of experimental data were additionally used. The obtained results showed that the model satisfactorily described the test data. The validated kinetic model can be used for optimal control of batch and fed-batch fermentation.

Keywords:

• Multi-objective optimisation
• parameter estimation
• batch fermentation
• kinetic model
• pareto optimal solution
• kluyweromyces marxianus var. lactis MC5

Introduction

Biotechnological processes, as an object for modelling and optimisation, differ significantly from the processes taking place in non-living nature, which is due to the complicated interaction of different influences. The growth of microorganisms and the metabolic transformations that occur in them is a multistage process with complex interrelated physiological, biochemical, genetic, physical and other factors. These features significantly complicate the process of their modelling and management. Based on these facts, the need to develop adequate models necessary for process optimisation and management emerges as an important problem. Mathematical modelling is an approach that allows studying the static and dynamic behaviour of the process without doing practical experimental studies, which are expensive and time-consuming to conduct [1–4].

In most cases, biotechnological processes (BTP) are described by unstructured models, such as a set of non-linear differential equations. Parameter estimation in non-structural models is important in the verification and subsequent use of a mathematical model for modelling, optimisation and optimal process control [5–10]. The ethanol production by a strain of Saccharomyces cerevisiae can be optimised by means of a general multiple objective optimisation framework for several metabolic responses involved in the ethanol production process [11]. Zhou et al. [12] have described a rapid method for the design of integrated bioprocesses using a combination of operating windows and a Pareto optimisation approach. The authors explored the use of a Pareto optimisation technique to find the optimal conditions for an integrated bioprocessing sequence and the benefits of first reducing the feasible space by developing a series of operating windows to provide a smaller optimisation search area.

In recent years, evolutionary methods such as hybrid differential evolution [13,14] and genetic algorithms (GAs) have attracted attention [15]. The GAs have become a popular universal tool for solving various optimisation problems and have been used for multi-objective optimisation of kinetic parameters estimations [16–18]. For example, a modified GA has been proposed for the parameter identification of an E. coli fed-batch process [19]. The results from the simulation demonstrate the high efficiency of the proposed modified GA for parameter identification of fermentation processes.

In addition to these methods, fuzzy sets theory is also widely used [20].

An interactive procedure has been applied to solve multi-objective optimisation problems [21]. A fuzzy set was used to model the decision maker’s judgment for each objective function. A fuzzy decision-making procedure has been utilised to find the optimal feed policy for the production of fuel ethanol via fed-batch fermentation [22,23]. This is done by means of an assigned membership function for each objective to convert the general multiple objective optimisation problem into a maximising decision problem.

Following the approach described in Wang et al. [22], Petrov and Ilkova [24] applied a fuzzy decision procedure to find the optimal feeding policy for a fed-batch fermentation process for L-lysine production in a stirred bioreactor from Brevibacterium flavum 22LD. The fermentation process was considered as a general multi-objective optimisation problem which was transformed into a single maximisation problem via an assigned membership function for each objective. A fuzzy set theory method was used to determine the global solution. This is a direct approach to obtain the solution of the optimisation problem. The applied multi-objective optimisation showed a significant increase in process productivity and a decrease in the concentration of glucose and threonine at the end of the process.

In a series of previous works [25–28] on the cultivation of the yeast Kluyweromyces marxianus var. lactis MC5, we have investigated ten non-structural models. The models were studied independently for the two main substrates of the process, lactose and oxygen. Two forms of specific growth rate, in multiplicative and additive form, are discussed in these reports. To select the most appropriate model, two multi-criteria decision-making methods have been applied, the Intercriteria decision analysis (ICDA) method [29] and the PROMETHEE II method [30,31]. The application of ICDA showed that the best model for the specific growth rate is the Monod – Monod model in its multiplicative form. The application of PROMETHEE II, in turn, showed that the most appropriate model is the Tessier – Monod model in its multiplicative form. All investigated combinations of models for specific growth rate are compared by eleven criteria.

Without publishing the results, the following combinations of models for specific growth rate were also investigated (in multiplicative form): Moser – Monod, Han-Levenspiel – Monod, Haldane – Haldane, Moser – Mink, Monod – Tessier, Tessier – Moser, and Tessier – Mink. In additive form, the following were studied: Levenspiel – Monod, Tessier – Moser, and Tessier – Monod.

The aim of this study was to determine the most appropriate of the above-cited models for the specific growth rate. Then, based on experimental data from batch cultivation, to determine the kinetic parameters of the model, by formulating a multi-objective optimisation problem. With the kinetic parameters thus determined, the model should be tested with additional experimental data.

Materials and methods

Experimental investigation

The strain Kluyveromyces marxianus var. lactis MC5 was cultivated under the following conditions [26,27,32]:

Nutrient medium with the basic component – whey ultra filtrate with lactose concentration 44 g L⁻¹. The ultra filtrate is derived from whey separated in the production of white cheese and deproteinisation by ultra filtration on LAB 38 DDS with a GR 61 PP membrane under the following conditions: T = 40–43 °C; input pressure, P_in = 0.65 MPa; output pressure, P_out = 0.60 MPa. The nutrient medium contained 0.6% (NH)₄HPO, 5% yeast autolysate and 1% yeast extract 1%; pH = 5.0–5.2.

The air flow rate (Q_G) was 60 L L⁻¹ h⁻¹ up to hour 4 and 120 L L⁻¹ h⁻¹ up to the end of the process under continuous mixing n = 800 min⁻¹.

Temperature was 29 °C.

The changes in the microbiological process (lactose conversion in yeast cells to protein) were studied during the strain growth: lactose concentration in fermentation medium in oxidation and assimilation of lactose by Kluyweromyces marxianus var. lactis MC5 was determined by enzyme methods by UV tests (Boehringer Manheim, Germany, 1983); concentration of cell mass and protein contents were determined on the basis of the nitrogen contents (Kjeltek system 1028); dissolved oxygen concentration in the fermentation medium in the process of oxidation and assimilation of lactose was determined by an oxygen sensor. For the measurement of the oxygen concentration in the fermentation medium, an oxygen sensor (LKB) was used.

The duration of the cultivation process was t_f= 12 h.

Six fermentations where carried out in aerobic batch cultivation. The experimental investigations were carried out in the computer controlled laboratory bioreactor 2 L-M with magnetic coupling [33] and the total volume was V = 2.0 L.

The results of the experimental studies for the six fermentations are shown in Supplemental Appendix (Table S1 and Figures S1–S6).

Kinetic model

Two methods (ICDA and PROMETHEE II) were used to select the most suitable combination of models for the specific growth rate.

The ICDA method can simplify the original Multi-criteria decision making (MCDM) problem by elimination of all redundant criteria or objects, thus resulting in a simplified (almost) equivalent MCDM problem that can be solved faster and with lower computational costs. In the particular case, the objects are the growth rate models separately for lactose and oxygen. The models studied and the validation criteria used are shown in the Supplemental Appendix. Through ICDA, the Monod model for the specific growth rate depending on lactose and oxygen in its multiplicative form was chosen. In this case, the model for the specific growth rate was: $\mu (S,C)=\frac{{\mu }_{m}S}{(K_S+S)}\frac{C}{(K_C+C)}$, h⁻¹.

The second method (PROMETHEE II) showed that the most appropriate form of the specific growth rate was Tessier’s model for lactose and Monod’s model for oxygen in its multiplicative form (Supplemental Appendix). In this case, the model for the specific growth rate was: $\mu (S,C)={\mu }_m\left(1-\hspace{0.17em}\mathit{\text{exp}}\hspace{0.17em}\left(-\frac{S}{K_{\mathit{\text{SI}}}}\right)\right)\frac{C}{(K_C+C)}$, h⁻¹.

The designations in the above models are given below.

In order to be unambiguous in choosing a model of the specific growth rate, we applied PROMETHEE II. In it, unlike ICDA, a strict preference can be set for the highest ranked alternative. This required all criteria to be selected Type I – a common criterion from the Preference Function included in the method. We used PROMETHEE Academic Edition software [34] to solve the multi-criteria decision-making problem.

The results for different growth rate models after the application of PROMETHEE II are shown in Table 1.

Table 1. Results after application of the PROMETHEE II method.

Rank	Model	Form	ϕ	ϕ^+	ϕ^-
1	Han-Levenspiel – Monod	ADD	0.6154	0.8042	0.1888
2	Tessier – Moser	ADD	0.3636	0.6783	0.3147
3	Tessier – Monod	ADD	0.2378	0.6084	0.3706
4	Tessier – Monod	MUL	0.1469	0.5734	0.4266
5	Mink – Haldane	MUL	0.1049	0.5455	0.4406
6	Moser – Monod	MUL	0.0559	0.5245	0.4685
7	Tessier – Moser	MUL	0.0350	0.5175	0.4825
8	Han-Levenspiel – Monod	MUL	−0.0769	0.4615	0.5385
9	Monod – Monod	MUL	−0.1329	0.4336	0.5664
10	Haldane – Haldane	MUL	−0.1399	0.4266	0.5664
11	Moser – Mink	MUL	−0.1469	0.4126	0.5594
12	Monod – Tessier	MUL	−0.1888	0.4056	0.5944
13	Mink – Mink	MUL	−0.3566	0.3147	0.6713
14	Tessier – Mink	MUL	−0.5175	0.2308	0.7483

ADD: additive form; MUL: multiplicative form of specific growth rate; ϕ⁺: positive outranking flow; ϕ⁻: negative outranking flow, outranking for PROMETHEE II method: $\phi ={\phi }^{+}-{\phi }^{-}$.

The highest-ranking model (Table 1) for specific growth rate is the Han-Levensipel – Monod model in its additive form. Second in rank is the Tessier – Moser model, but with a significantly smaller value of ϕ. For modelling the specific growth rate, we choose the Han-Levenspiel – Monod model.

The specific growth rate has the additive form $\mu (S,C)=\mu (S)+\mu (C)$ of Han-Levenspiel – Monod model: (1) $$\frac{\mathit{\text{dX}}}{\mathit{\text{dt}}}=[\mu (S)+\mu (C)]X$$(1) (2) $$\frac{\mathit{\text{dS}}}{\mathit{\text{dt}}}=-\frac{1}{Y_{X/S}}\mu (S)X$$(2) (3) $$\frac{\mathit{\text{dC}}}{\mathit{\text{dt}}}=-\frac{1}{Y_{X/C}}\mu (C)X+k_{l}a(C^{*}-C)$$(3) (4) $$\mu (S)=\frac{{\mu }_{1m}S{(1-S/S_m)}^n}{S+K_S{(1-S/S_m)}^m},\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\mu (C)=\frac{{\mu }_{2m}C}{(K_C+C)}$$(4) where: t – process time, h; X – biomass (dry weight), g L⁻¹; S – lactose concentration, g L⁻¹; C – oxygen concentration, g L⁻¹; µ(S) and µ(C) – growth rate of biomass only from lactose and oxygen, respectively, h⁻¹; Y_X/S – yield coefficients of formation of the biomass from lactose, g g⁻¹; Y_X/C – yield coefficients of formation of the biomass from oxygen, g g⁻¹; k_la – mass transfer coefficient, h⁻¹; C* – maximal oxygen concentration in liquid phase C(0) = C*, g L⁻¹; ${\mu }_{1m}$ and ${\mu }_{2m}$ – maximal values of growth rate, only from lactose and oxygen, respectively, h⁻¹; K_S and K_C – Monod saturation constant, g g⁻¹;S_m – critical inhibitor concentrations, above which the reactions stops, g L⁻¹;m, n – constants, –.

The initial conditions of the six batch fermentation were as previously described [26,27,32]: $$X\left(0\right)=0.20{\text{g L}}^{-1}\left(\text{dry weight}\right);S\left(0\right)=\text{44}.0{\text{g L}}^{-1};C\left(0\right)=\text{6}.65´10^{-3}{\text{g L}}^{-1}$$

In this study, experimental data from four batch fermentations of K. marxianus var. lactis MC5 were simultaneously used to determine the ten (n = 10) parameters in models (1)-(4). The selection of objective functions was performed in the next part of the work.

All studied models cited in our publications are presented in the Supplemental Appendix (Table S2). The statistical criteria by which the models were evaluated using the ICDA method and PROMETHEE II method are also shown there (Table S3).

Objective functions for determination of kinetic parameters

In this work, the approach proposed by Wang and Shew [5] was used. In this approach, a vector of objective functions is considered. Each of these objective functions is simultaneously used to estimate the parameters of the kinetic model, thereby transforming the parameter estimation problem into a multi-objective optimisation problem. The vector of objective functions was: (5) $$\underset{\theta }{\mathit{\text{min}}}Q=\underset{\theta }{\mathit{\text{min}}}{\left[Q_1,\dots .,Q_k,\dots ,Q_m\right]}^T$$(5) where θ is denoted as n × 1 vector of kinetics parameters in the model, m – number of batch experiments (m = 4).

The vector of objective functions (Q) consists of the criteria of the four batches of observations. Here, Q_k refers to an objective function for the k^th observation of a batch, respectively. Each objective function for the batch observation is expressed by the mean least squares error as (6) $$Q_k(\theta )=\frac{1}{N}\sum _{i=1}^N\left(\begin{array}{c}\frac{{[X_m(t_i,\theta )-X_e(t_i)]}^2}{X_{e\mathrm{\text{max}}}^2}+\frac{{[S_m(t_i,\theta )-S_e(t_i)]}^2}{S_{e\mathrm{\text{max}}}^2}+\\ +\frac{{[C_m(t_i,\theta )-C_e(t_i)]}^2}{C_{e\mathrm{\text{max}}}^2}\end{array}\right)$$(6) where: Q_k(θ) – objective function of each batch experiment; k = 1, …., m; N – number of the experiments; t_i – time partitions, h; $X_m(t_i,\theta ),S_m(t_i,\theta ),C_m(t_i,\theta )$ – simulated values with models of kinetics variables for each experiment, g L⁻¹; $X_e(t_i),S_e(t_i),C_e(t_i)$ – measurement values of kinetics variables for each experiment, g L⁻¹; $X_{e\mathrm{\text{max}}}^2,S_{e\mathrm{\text{max}}}^2,C_{e\mathrm{\text{max}}}^2$ – maximal measurement values of kinetics variables for each experiment, g L⁻¹.

Pareto optimal solution

In multi-objective optimisation, we have multiple competing functions. It is most appropriate for multi-objective optimisation to use Pareto optimal solutions [5,9–12]:

Definition 1.

The feasible region in the input space of model parameters Ω is the set of all admissible model parameters that satisfy the model Equations (5)(5) $$\underset{\theta }{\mathit{\text{min}}}Q=\underset{\theta }{\mathit{\text{min}}}{\left[Q_1,\dots .,Q_k,\dots ,Q_m\right]}^T$$(5) , i.e. $$\Omega =\left\{\theta \dot{z}=f(z,\theta ),z(0)=z_0\text{and}{\theta }_{\mathit{\text{min}}}\le \theta \le {\theta }_{\mathit{\text{max}}}\right\}$$

where θ_min and θ_max are the lower and upper bounded vectors of the model parameters, $\dot{\mathbf{z}}=\mathbf{f}(\mathbf{z},\theta )$ consists of the batch model equations.

Thus, the Pareto optimal solutions to the multi-objective parameter estimation problem can be defined.

Definition 2.

The vector of model parameters θ* is Pareto optimal for the problem Equation (5)(5) $$\underset{\theta }{\mathit{\text{min}}}Q=\underset{\theta }{\mathit{\text{min}}}{\left[Q_1,\dots .,Q_k,\dots ,Q_m\right]}^T$$(5) if and only if there does not exist $\theta \in \Omega$ such that $Q_k(\theta )\le Q_k({\theta }^{*})$, for k = 1, …, m with $Q_j(\theta )<Q_j({\theta }^{\mathbf{*}})$ at least one j [4].

The image of a Pareto optimal point is a Pareto optimal solution.

A weighted min-max method [4] was used in this work. This method for characterising Pareto optimal solutions is for solving the following weighted min-max problem: (7) $$\underset{\theta \in \Omega }{\mathrm{\text{min}}}\underset{k}{\mathrm{\text{max}}}\left[w_k{Q}_k,k=1,\dots ,m\right]$$(7)

The weighted min-max problem means that the model parameters are determined in terms of the worst fit result. In Wang and Shew [5], theorems are proved which are used to ensure that the optimal solution of the weighted min-max problem is a Pareto optimal solution of the multi-objective optimisation problem Equations (5)(5) $$\underset{\theta }{\mathit{\text{min}}}Q=\underset{\theta }{\mathit{\text{min}}}{\left[Q_1,\dots .,Q_k,\dots ,Q_m\right]}^T$$(5) .

Results

In this work, experimental data from four batches were used to estimate the parameters of the kinetic model, and data from two batch experiments were used for testing.

As the approach that was used here was that of Wang and Sheu [5], the multi-objective parameter estimation problem in Equations (5)(5) $$\underset{\theta }{\mathit{\text{min}}}Q=\underset{\theta }{\mathit{\text{min}}}{\left[Q_1,\dots .,Q_k,\dots ,Q_m\right]}^T$$(5) was first converted to the weighted min-max problem in Equations (7)(7) $$\underset{\theta \in \Omega }{\mathrm{\text{min}}}\underset{k}{\mathrm{\text{max}}}\left[w_k{Q}_k,k=1,\dots ,m\right]$$(7) . This approach assumes that we are minimising the worst-case objective function in the parameter estimation problem.

In this work, to deal with the multi-objective optimisation problem, we used the weighted sum method (WSM) [10], as the simplest and most widely used approach, to minimise the composite function Q_WSM, which is the weighted sum of the objectives: (8) $$\underset{\theta }{\mathrm{\text{min}}}{Q}_{\mathit{\text{WSM}}}=\sum _{k=1}^m{v}_k{Q}_k(\theta ),v_k>0;k=1,\dots ,m$$(8) where v_k are weights which represent the relative importance of each objective.

If the weights are positive, minimising provides a Pareto optimal solution.

In our work, we assumed that all objective functions $Q_k(\theta )$ have equal weights, i.e. $v_k=(1/m)=0.25$, k = 1, …, 4; m = 4.

The min-max problem Equations (7)(7) $$\underset{\theta \in \Omega }{\mathrm{\text{min}}}\underset{k}{\mathrm{\text{max}}}\left[w_k{Q}_k,k=1,\dots ,m\right]$$(7) and the determination of the optimal solution by Equations (8)(8) $$\underset{\theta }{\mathrm{\text{min}}}{Q}_{\mathit{\text{WSM}}}=\sum _{k=1}^m{v}_k{Q}_k(\theta ),v_k>0;k=1,\dots ,m$$(8) were solved by the direct search method. A double-precision BCPOL subroutine in IMSL Math/Library [35] was used to estimate the kinetic parameters, with a stopping tolerance 10⁻⁷ and a maximum number of iterations 10,000. Increasing the tolerance to 10⁻¹⁰, did not significantly improve the results.

The algorithm was implemented by the CОMPAQ Visual FORTRAN 90 [36]. All computations were performed on an AMD Phenom II X6 1100 T (3.3 GHz) computer using Microsoft Windows XP Pro Edition (32 bit) operating system.

The weighted min-max problem Equations (7)(7) $$\underset{\theta \in \Omega }{\mathrm{\text{min}}}\underset{k}{\mathrm{\text{max}}}\left[w_k{Q}_k,k=1,\dots ,m\right]$$(7) is applied to estimate the kinetic parameters in models (1)–(4) for different weights shown in Table 2. These weights are used to illustrate the expected effects.

Table 2. Various weights for 4^th batch process.

Weight					Objective functions, Qk(θ)×10^-3
Case	1^st	2^nd	3^rd	4^th	1^st	2^nd	3^rd	4^th
1	1	1	0	0	5.941	5.736	8.462	9.496
2	0	0	1	1	5.735	3.914	7.537	8.314
3	1	1	1	1	5.941	5.736	8.462	9.496
4	1	1	1	1	3.685	5.698	3.830	3.292
5	1	0	0	0	5.941	5.736	8.462	9.496
6	0	1	0	0	5.869	3.929	7.475	8.171

In the first case in Table 2, weights are assigned only to the 1st and 2nd batch experiments and zero to the 3rd and 4th. Therefore, only the first two pieces of experimental data are used to estimate the parameters of the kinetic model. The BCPOL settings are shown above. They are the same for all cases. The objective function values for each experiment are shown in Table 2 (row 1).

In the second case of Table 2, the weights of the 3rd and 4th batch of experiments are ones and zero for the 1st and 2nd batches. Therefore, only the 3rd and 4th batch of experimental data are used to estimate the parameters of the kinetic model. The objective function values for each experiment are shown Table 2 (row 2). From the obtained results, it is observed that the objective functions of the 3rd and 4th experiments are smaller than those of Case 1.

In the third case of Table 2, all weights are equal to one. Therefore, all four batches of experimental data are used to estimate the parameters of the kinetic model. The objective function values for each experiment are shown in Table 2 (row 3).

From the results shown in Table 2, it can be seen that the objective function values for case 1 and case 3 (rows 1 and 3) are equal. This should show that the result is not significantly affected by the 3rd and 4th batch of experiments. This means that the first and second batch experiments have top priority. To clarify this, Case 5 and Case 6 were also investigated (Table 2). The obtained results show that the main priority in solving the min-max problem Equation (7)(7) $$\underset{\theta \in \Omega }{\mathrm{\text{min}}}\underset{k}{\mathrm{\text{max}}}\left[w_k{Q}_k,k=1,\dots ,m\right]$$(7) is the first batch experiment.

The results obtained from the third case in Table 2 were used to determine the Pareto optimal solution Equations (8)(8) $$\underset{\theta }{\mathrm{\text{min}}}{Q}_{\mathit{\text{WSM}}}=\sum _{k=1}^m{v}_k{Q}_k(\theta ),v_k>0;k=1,\dots ,m$$(8) . The solution is performed by BCPOL with the same tuning parameters. This is shown in Case 4 (Table 2, row 4).

The calculated optimal kinetic parameters are shown in Table 3.

Table 3. Optimal estimated parameters of the fourth case.

Parameter	Search bound	Optimal values
${\mu }_{1\mathrm{\text{max}}}$	[0.1, 1.0]	0.367
KS	[5.0, 150.0]	143.080
Sm	[45, 150]	126.730
n	[0.1, 3.0]	1.090
m	[0.1, 3.0]	2.125
${\mu }_{2\mathrm{\text{max}}}$	[0.1, 1.0]	0.633
KC	[0.1, 15.0] × 10^−3	0.715 × 10^−3
YX/S	[0.01, 2.5]	0.056
YX/C	[1.0, 25.0]	3.419
kla	[80.0, 180.0]	145.769

The calculated profiles, with the determined optimal parameters of the kinetic model, for biomass, lactose and oxygen are shown in Figures 1–4.

Figure 1. Concentration of biomass, lactose and oxygen for the 1st batch experiment. Xe_j, Se_j and Ce_j – the experimental data for the four batch experiments, j = 1, …, 4; X_k,j – biomass concentration results obtained by using the model parameters for the first (k = 1) and the second case (k = 2), j = 1, …, 4. The indices ‘_k,j’ for the concentration of lactose (S_k,j) and oxygen (C_k,j) are the same as those for X_k,j; Xm_j, Sm_j and Cm_j – simulated results obtained by the optimal solutions (case 4 in Table 2) of the model.

Figure 2. Concentration of biomass, lactose and oxygen for the 2nd batch experiment.

Figure 3. Concentration of biomass, lactose and oxygen for the 3rd batch experiment.

Figure 4. Concentration of biomass, lactose and oxygen for the 4th batch experiment.

The experimental data from the fifth and sixth fermentations were used to test the developed mathematical model. The obtained results are shown in Figures 5 and 6.

Figure 5. Concentration of biomass, lactose and oxygen for the tested results of the 5th batch experiment. Xe_j, Se_j and Ce_j – experimental data of biomass, lactose and oxygen concentration for 5th and 6th test batch processes; Xm_j, Sm_j and Cm_j (j = 5, 6) – simulated data, respectively.

Figure 6. Concentration of biomass, lactose and oxygen for the tested results of the 6^th batch experiment.

From the obtained results, it can be seen that for all batch processes we have a decrease in the value of the objective functions.

From Figure 1 (first batch experiment), it is observed that for lactose, the optimal solution is worse (solid line) than the solutions for Case 1 and Case 2. For the second batch experiment, the results for the optimal solution are better than are the solutions obtained by Case 1 and Case 2.

In the third and fourth sets (Figures 3 and 4) of experimental data, the situation is the same as in the first batch experiment. The results are worse for the optimal solution than for the solutions obtained from Case 1 and Case 2. In general, it can be seen that the calculated profiles (Figures 1–4) agree satisfactorily with the batch experiments.

Very good results (Figures 5 and 6) were obtained for the test experiments. The computational results in the batch model fit the experimental data satisfactorily.

Discussion

In all our previous publications dealing with the process of cultivation of the yeast Kluyveromyces marxianus var. lactis MC5, different optimisation and optimal control methods were used [37–42]. In these works on the specific growth rate, the same dependence was used – the Mink model for lactose and the Haldane model for oxygen, in its multiplicative form.

For the first time in 2015 [33], the modelling of the process was done separately for the two main substrates – lactose and oxygen. Nine growth rate models were examined. From these models, based on the first five criteria (Supplemental Appendix) for assessment of adequacy, five combinations of the specific growth rate in its multiplicative form were selected and investigated. In the same year, the ICDA method was used for model selection. Two sets of models emerged from this research. The details of the model selection can be found in [25]. The research continued, adding a tenth model to the researched models. Again, two groups of best-performing models formed. The details of the model selection can be found in [26]. The completion of these studies, using the ICDA method, was done in 2021 [32], where the ten growth rate models were again examined. In this publication, the comparison between the multiplicative and the additive form of the specific growth rate was also made for the first time. The obtained results showed that the additive form is the most suitable for modelling the specific growth rate. In the same year, another multi-criteria decision-making method was used [28]. The conflicting results obtained by the ICDA and PROMETHEE II method led to the present work.

In this work identification of parameters in the model is based on the multi-criteria identification of parameters in kinetic models of biotechnological processes described in Chen and Wang [22]. The approach was applied in the multi-criteria parametric identification of batch processes in the cultivation of K. marxianus var. lactis MC5. The identification of the parameters in the model is based on the experimental data from four batch processes with simultaneous evaluation of the kinetic parameters of the model. Tests with additional experimental data from two batch process models showed very good agreement between the experimental data and the data simulated by the model.

The PROMETHEE II method developed by Brans, Mareschal and Vincke [30, 31] for multi-criteria decision making was used to select the most appropriate model of specific growth rate from 14 combinations of models in batch cultivation of K. marxianus var. lactis MC5. The application of PROMETHEE II showed that the Han-Levenspiel–Monod model in its additive form has the highest rank. This form of specific growth rate was then applied to identify the parameters in the kinetic model of the process.

Solving the task of multi-criteria parametric identification, testing the model with the additional fifth and sixth fermentations, solving the system of differential equations, was carried out using the purpose-developed algorithm and COMPAQ Visual FORTRAN 90 program [36]. The developed algorithm and program can be used to solve similar tasks in the field of biotechnology.

In future studies of the cultivation process of the yeast Kluyweromyces marxianus var. lactis MC5 it is useful to use methods that allow the determination of a global extremum (evolutionary algorithms, such as GAs, etc.) since the direct search method is highly dependent on the default search values. It is also recommended to explore other methods for the determination of the Pareto optimal solution, as well as to analyse the sensitivity of the kinetic parameters in the model. In this way, an additional assessment of the suitability of the proposed model will be made.

Conclusions

In this study, by using multi-criteria decision making, based on experimental data from batch cultivation, a simultaneous estimation of the kinetic parameters of the model was made. Based on statistical criteria, the most appropriate model was determined for the specific growth rate of K. marxianus var. lactis MC5. Additionally, two more batch processes with different experimental data were used to evaluate and test the model. The obtained results show that the model satisfactorily describes all six batches of experimental observations. The model can be successfully used for optimal feed rate control in fed-batch fermentation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data that support the findings reported in this study are available from the corresponding authors [MP] and are shown in the Supplemental Appendix, which is an integral part of the article.

Additional information

Funding

This work was supported by TBEQ-2024-C16214.