Adaptive Monitoring of Biotechnological Processes Kinetics

: In this paper, an approach for the monitoring of biotechnological process kinetics is proposed. The kinetics of each process state variable is presented as a function of two time-varying unknown parameters. Fortheirestimation, ageneralsoftwaresensorisderivedwithon-linemeasurementsasinputs that are accessible in practice. The stability analysis with a different number of inputs shows that stability can be guaranteed for fourth- and fifth-order software sensors only. As a case study, the monitoring of the kinetics of processes carried out in stirred tank reactors is investigated. A new tuning procedure is derived that results in a choice of only one design parameter. The effectiveness of the proposed procedure is demonstrated with experimental data from Bacillus subtilis fed-batch cultivations.


Introduction
Monitoring of the main variables and parameters of biotechnological processes is of key importance for processes investigations and control, especially in industrial installations, where a limited number of measurements is available [1,2].
A widely used approach for process monitoring is model-based software sensor (SS) design [3,4]. The operational model has to be (i) as accurate as possible to mimic the main characteristics and dynamics of the process, and (ii) simple enough for a monitoring and control design. Such an operational model was proposed by Bastin and Dochain [4] as the general dynamical model (GDM) of bioprocesses in stirred tank reactors (STR). Approaches based on GDM software sensors are widely applied simultaneously with other approaches for nonlinear systems, such as extended Kalman and Luenberger filters [5][6][7][8], moving horizon [9,10] neural-network based observers [11], high-gain approach [12], multirate observers [13], sliding mode-observers [14,15], interval SS [16], cascade SS [17][18][19], and joint estimation of state variables and parameters [1,20,21], among others.
Usually, the software sensors are designed using operational models with constant yield coefficients [4,12,15,20,21]. For many industrial biotechnological processes, like wastewater treatment and processes in inhomogeneous mediums, the reproducibility is poor [18,22]. Hence, the assumption of a constant yield coefficient leads to inexact results because of considerable changes of these parameters during a process or within different production batches. This change is due to adaptations of metabolic pathways, protein expression pattern, and random mutations of organisms, as well as the occurrence and dynamics of population heterogeneities in single species, especially in multispecies bioprocesses [23,24].

Monitoring of the Unknown Kinetics
The proposed approach can be considered as an extension of the previously described method of a model-based software sensor [4].

New Formalization of the Process Kinetics
The bioprocess kinetics describe the rates (fluxes) of biochemical reactions of the process as functions of concentrations, pH, temperature, etc.
In the considered case, it is presented as product of two unknown time-varying terms as follows: where φ m (t)-the vector of known kinetics; ϕ(t)-the key kinetic parameter, which describes the dynamics of the main state variables; and Y(t)-the vector of yield coefficients comprising remaining parts of the state variables' kinetics. This presentation of the kinetics is a novelty, as both terms are unknown time-varying parameters with a clear physical meaning. The maximal dimension, m, of the vector of known kinetics, φ m (t), and the vector of unknown yield coefficients, Y(t), is three. In the case of m > 3, problems of the software sensor design arise as it is shown below.

Structure of the General Software Sensor
For linearization of relations between yield coefficients and the key parameter, the following vector of auxiliary parameters W Y and the auxiliary parameter W ϕ are introduced: where Through the differentiation of Equations (1)-(3), the following dynamic equations of φ m (t) elements are obtained: The dynamics of the vector W Y can be presented as a function of φ m (t), its time-derivative and time-derivative of W ϕ , as follows: Applying the natural logarithm to the vector φ m (t) (Equations (1)-(3)), the following relationship is obtained: Equations (4)-(6) determine the linear structure of the software sensor of parameters W Y and W ϕ , described with the following equations: ω, γ, and Ψ are a diagonal matrix and vectors, of which the elements are tuning parameters of software sensor system (7), respectively. Since Equations (8) and (9) are scalar quantities, system (7) is of fourth or fifth order depending on the number of measurements. The structure is derived by the assumption that according to Equation (9). The estimatesŶ(t) andφ(t) are obtained by applying reverse functions on Equations (2) and (3) as follows: The general structure of software sensor Equations (7)-(9) is given in Figure 1. The inputs and kinetics formalization are presented by Equation (1). The auxiliary variables are calculated using Equations (2) and (3). Linearization is carried out according to Equations (4)-(6). The software sensor is described by system (7)(8)(9). The outputs are obtained by the reverse transformations using Equations (10) and (11). ω, γ, and Ψ are a diagonal matrix and vectors, of which the elements are tuning parameters of software sensor system (7), respectively. Since Equations (8) and (9) are scalar quantities, system (7) is of fourth or fifth order depending on the number of measurements. The structure is derived by the assumption that Ŵφ = d p dt and the dynamics of the estimate of p depends on the estimation errors according to Equation (9). The estimates ( ) The general structure of software sensor Equations (7-9) is given in Figure 1. The inputs and kinetics formalization are presented by Equation (1). The auxiliary variables are calculated using Equations (2) and (3). Linearization is carried out according to Equations (4)-(6). The software sensor is described by system (7)(8)(9). The outputs are obtained by the reverse transformations using Equations (10) and (11). The application of the original idea for the linearization of kinetics, by representing it in logarithmic form, leads to the derivation of a linear structure of the software sensor, system (7)(8)(9). This allowed us, by applying the well-known linear control theory, to study the stability of both SS structures (fourth-and fifth-order), as well as to reduce the number of SS' design parameters to be tuned.

Stability Analysis
The following system of the estimation error dynamics is considered: where x-the estimation error vector with dim = 4 or 5;  The application of the original idea for the linearization of kinetics, by representing it in logarithmic form, leads to the derivation of a linear structure of the software sensor, system (7)(8)(9). This allowed us, by applying the well-known linear control theory, to study the stability of both SS structures (fourthand fifth-order), as well as to reduce the number of SS' design parameters to be tuned.

Stability Analysis
The following system of the estimation error dynamics is considered: where x-the estimation error vector with dim = 4 or 5; W Y , W ϕ , and p-'true' values of the parameters; W Y ,Ŵ ϕ , andp-estimates of W Y , W ϕ , and p; where ω i , γ i , and ψ i are the SS design parameters. The concrete presentation of matrix A is: The analysis of the observability shows that the system (7-9) (for m = 1, 2, and 3) was not completely observable, but detectable, which is equivalent to the existence of an asymptotically stable linear observer of the state [28]. Therefore, the parameter tuning of estimator (7-9) cannot guarantee exponential (arbitrarily chosen) convergence of estimates, but asymptotic convergence only.
The cases of m = 2 (fourth-order system) and m = 3 (fifth-order system) are considered below. In the case of m > 3, the polynomial degree of matrix A is greater than five, which leads to problems of the software sensor design [29].

Stability Analysis of the Fourth-Order System
The system stability is investigated based on the analysis of system (12), where matrix A and vectors x and u are presented as: The characteristic polynomial of matrix A is: where λ is the variable of the characteristic polynomial. If matrix A has four different Eigenvalues (λ 1 ÷ λ 4 ), the following system (including equalities and inequalities) has to be solved to guarantee the stability: Numerous variations of system (16)(17)(18)(19) exist that lead to a solution of a system of four algebraic equations with four Eigenvalues and six tuning parameters.

Stability Analysis of the Fifth-Order System
Matrix A and vectors x and u are presented as follows: The characteristic polynomial of matrix A is: In the case of the five different Eigenvalues λ 1 ÷ λ 5 of matrix A, the following system has to be solved to guarantee stability: The general presentation of the fifth-order polynomial cannot be solved algebraically using finite number operations [29]. Many solutions of system (22)(23)(24)(25)(26) exist, presenting five algebraic equations with five Eigenvalues and nine tuning parameters. The parameter's tuning could be different depending on the considered process and its dynamics, as it is shown in the example below.

Structure of the General Software Sensor
As a case study, a class of aerobic fed-batch processes carried out in homogenous conditions was considered. An operational model was derived. The model describes the dynamics of the key variables chosen from the expert's point of view, as follows: where X, S 1 , S 2 , and RQ are biomass, glucose, dissolved oxygen concentrations, and the respiratory quotient, respectively; F-glucose feed rate; V-culture volume; S 1m -glucose concentration in the feed; -rates of biomass growth, substrate consumption and dissolved oxygen consumption, respectively; q RQ (t)-the change of the ration between the specific rates of carbon dioxide production and oxygen consumption q CO 2 (t)/q S 2 (t).
In the considered case, a linear relationship exists between the experimental data of the specific oxygen consumption rate and the specific carbon dioxide producti on rate. This leads to problems in software sensor design. For this reason, the respiratory quotient RQ was used, as calculated from the on-line data of the off-gas analysis. As shown below, the relationship between the yield coefficient,Y RQ/S 1 (t), related to RQ, and yield coefficients with biological meaning, Y X/S 1 and Y eth/pyr , is justified.
In practice, for the considered processes, on-line measurements of X, φ S 2 (t), and RQ were available. Using the observer-based estimators published by the authors of [16], from on-line measurements of X and RQ, we received the rates φ X (t) and q RQ (t), which together with φ S 2 (t), were accepted as inputs of the scheme from Figure 2.
where X, S1, S2, and RQ are biomass, glucose, dissolved oxygen concentrations, and the respiratory quotient, respectively; F-glucose feed rate; V-culture volume; S1m-glucose concentration in the feed; Kla-oxygen transfer coefficient; 2 * S C -oxygen saturation concentration; ( ) ( ) ( ) 1 2 , , In the considered case, a linear relationship exists between the experimental data of the specific oxygen consumption rate and the specific carbon dioxide producti on rate. This leads to problems in software sensor design. For this reason, the respiratory quotient RQ was used, as calculated from the on-line data of the off-gas analysis. As shown below, the relationship between the yield coefficient,   , , comprising the remaining parts of the corresponding kinetics: Equations (30) and (34) were derived assuming that the changes of the respiratory quotient can be expressed as: Following the general formalization, each rate, φ X , φ S 2 , and q RQ , is presented as function of two terms considered as unknown time-varying parameters: The glucose consumption φ S 1 (t) as common term and a yield coefficient Y X/S 1 (t), Y S 2 /S 1 (t), Y RQ/S 1 (t) , comprising the remaining parts of the corresponding kinetics: Equations (30) and (34) were derived assuming that the changes of the respiratory quotient can be expressed as: where Y eth/pyr is the yield coefficient for ethanol production at the pyruvic acid branch point that is specific for the considered processes. Since the respiratory quotient RQ characterizes the degree of fermentative growth, it is postulated that its changes are proportional to the split of carbon toward ethanol synthesis at the pyruvic acid branch point under substrate limited conditions. Hence, the yield Y RQ/S 1 (t) represent the term of Equation (34), including both 'true' yield coefficients, Y X/S 1 and Y eth/pyr , as follows: The monitoring of the kinetics of the processes described by system (32)-(36) was reduced to the estimation of four parameters-the yield coefficients Y X/S 1 (t), Y S 2 /S 1 (t), Y RQ/S 1 (t) and the glucose consumption rate φ S 1 (t), which are outputs of the scheme shown in Figure 2. This figure repeats the structure of the general SS from Figure 1, where the information about each component corresponds to the considered case. System (7-11), describing the general SS for the considered case, can be rewritten as follows: This represents a fifth-order SS of the auxiliary variables W Y X/S 1 , W Y S 2 /S 1 , W Y RQ/S 1 , W φ S 1 . The values of φ X and q RQ , as well as their time-derivatives at the right side of Equations (37) and (39), were obtained as explained above. The time-derivatives of φ S 2m in (38) were calculated by numerical differentiation.
The outputs of the SS,Y X/S 1 (t), Y S 2 /S 1 (t), Y RQ/S 1 (t), and φ S 1 (t), are obtained by the reverse functions:Ŷ

Tuning Procedure
The structure of the derived software sensor systems (37)-(45) was simplified in comparison to the general structure (7)- (11). An optimal dimension of vectors γ and Ψ was proposed with the aim to facilitate the tuning procedure. It was based on stability analysis of the system's error.
Taking into account the dynamics of processes in homogenous conditions, it is accepted that: The value of λ can be obtained by an optimization procedure (minimization of the square of sum of estimation errors) using experimental data of each of the investigated processes of the class described by system (27)-(31): where φ X (t) andq RQ (t) are estimates of φ X and q RQ , obtained similarly to [16]. EstimatesŶ X/S 1 (t),Ŷ S 2 /S 1 (t),Ŷ RQ/S 1 (t), andφ S 1 (t) were received by the system (37)-(45).
The calculation of the design parameters is based on the following stability conditions: Based on (54) and (55), the following relation between the design parameter ω 3 and the Eigenvalue is obtained: The following auxiliary parameters D, B, and C are introduced: Substituting the terms ω 2 ω 3 and (ω 2 + ω 3 ) from (52) with their equivalents from the Equations (53) and (51), respectively, the following third order polynomial of parameter ω 1 is derived: One of the roots of (60) is denoted as: The design parameter ω 2 is calculated from (51): The values of the parameters ω 3 and ω 1 are obtained from Equation (56) and Equation (61), respectively. The value of the parameter γ 2 is calculated from (55) by: To satisfy the stability condition (55), the value of ψ 1 has to be negative. For the parameter γ 1 , we assumed the relationship γ 1 = γ 2 > 0.
Realizing the steps from Equation (51) to Equation (63), the values of all tuning parameters were received. The proposed tuning procedure reduced the SS design to choose of one tuning parameter only.
A fed-batch cultivation process of B. subtilis, carried out in the laboratory bioreactor, was considered as being representative. Two experiments (Experiment I and Experiment II), carried out in STR, were used for the demonstration of theoretical results [30].
The calculation of the specific oxygen uptake rate (q S2m ) and carbon dioxide production rate (q CO2m ) was performed on the basis of a gas balance of the input and output flow using O 2 and CO 2 measurements in the off-gas. The specific oxygen uptake rate, cell dry weight, and respiratory quotient RQ = q CO2m /q S2m were available as on-line measurements, which were used in the calculation of the SS inputs.
The two sets of experimental data are shown in Figure 3. On-line measurements are shown in Figure 3a,c,d. In Figure 3b, off-line measurements of glucose concentration are shown. Both experiments were performed under the same conditions, but as it can be seen, they varied due to a lack of reproducibility of experiments. In order to demonstrate adaptive performance of the proposed SS shown in Figure 2, the tuning of software sensors was conducted using the data of Experiment I, while the data of the Experiment II were used for validation.
To satisfy the stability condition (55), the value of ψ1 has to be negative. For the parameter γ1, we assumed the relationship γ1 = γ2 > 0.
Realizing the steps from Equation (51) to Equation (63), the values of all tuning parameters were received. The proposed tuning procedure reduced the SS design to choose of one tuning parameter only.
A fed-batch cultivation process of B. subtilis, carried out in the laboratory bioreactor, was considered as being representative. Two experiments (Experiment I and Experiment II), carried out in STR, were used for the demonstration of theoretical results [30].
The calculation of the specific oxygen uptake rate (qS2m) and carbon dioxide production rate (qCO2m) was performed on the basis of a gas balance of the input and output flow using O2 and CO2 measurements in the off-gas. The specific oxygen uptake rate, cell dry weight, and respiratory quotient RQ = qCO2m / qS2m were available as on-line measurements, which were used in the calculation of the SS inputs.
The two sets of experimental data are shown in Figure 3. On-line measurements are shown in Figure 3a,c,d. In Figure 3b, off-line measurements of glucose concentration are shown. Both experiments were performed under the same conditions, but as it can be seen, they varied due to a lack of reproducibility of experiments. In order to demonstrate adaptive performance of the proposed SS shown in Figure 2, the tuning of software sensors was conducted using the data of Experiment I, while the data of the Experiment II were used for validation.
For investigations, a program package was prepared in a MATLAB environment (MATLAB R2019a, the MathWorks, Inc., Natick, Massachusetts, United States, 2019). The package included software sensor systems (37)-(45) and both sets of experimental data from Figure 3. The value λ = −42, used in the tuning procedure, was obtained based on the optimization procedure with the criteria given in Equation (47-50). , , Estimates of Y X/S 1 (t), Y S 2 /S 1 (t), Y RQ/S 1 (t), and φ S 1 (t) are shown in Figures 4 and 5. The estimation results under the Experiment I are shown in Figure 4. The verification results are presented in Figure 5. To show the accuracy of the outputs, their reference values were derived using laboratory measurements of glucose concentration, S 1 . A discrete observer-based estimator allowed us to obtain discrete estimates of glucose consumption rate, R S − re f = φ S 1− re f , for both experiments. The reference values of Y X/S 1 (t), Y S 2 /S 1 (t), Y RQ/S 1 (t) were calculated using relationships (32)-(34) and φ S 1− re f . All reference values are shown in Figures 4 and 5 with points.
in Figure 5. To show the accuracy of the outputs, their reference values were derived using laboratory measurements of glucose concentration, S1. A discrete observer-based estimator allowed us to obtain discrete estimates of glucose consumption rate,   As it can be seen from the results in Figures 4 and 5, the dynamics of the yield coefficients and the substrate consumption rate had different profiles for the two experiments, although the experiments were realized under homogeneous conditions. This can be explained with the dynamic Y X/S1 Y S2/S1 Y RQ/S1 Glucose consumption rate (g/lh) Y X/S1 Y S2/S1 Y RQ/S1 Glucose consumption rate (g/lh) measurements of glucose concentration, S1. A discrete observer-based estimator allowed us to obtain discrete estimates of glucose consumption rate,   As it can be seen from the results in Figures 4 and 5, the dynamics of the yield coefficients and the substrate consumption rate had different profiles for the two experiments, although the experiments were realized under homogeneous conditions. This can be explained with the dynamic Y X/S1 Y S2/S1 Y RQ/S1 Glucose consumption rate (g/lh) Y X/S1 Y S2/S1 Y RQ/S1 Glucose consumption rate (g/lh) As it can be seen from the results in Figures 4 and 5, the dynamics of the yield coefficients and the substrate consumption rate had different profiles for the two experiments, although the experiments were realized under homogeneous conditions. This can be explained with the dynamic nature of bioprocesses and fluctuations in the strain performance. Comparing the two experiments, the values of the yield coefficients retained the trend of change to some extent, while the rate of consumption of the substrate slightly decreased in Experiment 1 (Figure 4d), but stayed more or less constant in Experiment 2 (Figure 5d).
The results, shown in Figures 4 and 5, demonstrate the good performance and adaptive properties of the proposed SS. As can be expected, the accuracy of the estimation and their convergence were better for the experiment used for the tuning of the estimator. The estimates converged to the corresponding reference values with a mean relative error below 3%.
The simulations proved the reliability of the proposed methodology for real-time assessment. By setting only one design parameter, a good accuracy of the estimates of each term of the kinetics with a clear physical meaning was achieved. By verifying with data from another experiment, the adaptability of the proposed approach was demonstrated.

Conclusions
In this paper, a general approach for kinetics monitoring of a class of biotechnological processes was proposed. This class of processes is characterized by a relationship between measured and estimated parameters that guarantees the synthesis of asymptotically stable SS.
The originality of this approach consists of its (i) presentation of the process kinetics with two unknown time-varying parameters with clear physical meaning, (ii) derivation of a linear structure of the SS using logarithmic transformations of these parameters that facilitate stability analysis, and (iii) derivation of stable fourth-and fifth-order structures of the SS, satisfying conditions for asymptotical stability.
The proposed approach was applied for the monitoring of process kinetics realized in homogeneous conditions. The tuning procedure was reduced to one design parameter, which led to an automatic calculation of the others. The adaptability of the proposed method was proven by its application to two B. subtilis fed-batch cultivations. The tuning of the SS was realized with one set of experimental data, and its validation was demonstrated by the other dataset. The accuracy of the obtained estimates is good enough for monitoring of the process dynamics and can be applied for process investigations and control design.
In comparison to our previous investigations, the presented methodology provides the possibility to reduce the tuning procedure to only one tuning parameter.
The obtained results are a good basis to extend both the theoretical investigations and their applications in the future. One of the possible directions of research is the integration of physiological parameters that can be measured on-line. This will expand the applicability of the presented approach, as there is a strong relationship between yield coefficient dynamics and the physiological state of cells.