Metabolic Reaction Network-Based Model Predictive Control of Bioprocesses

Abstract

Bioprocesses are increasingly used for the production of high added value products. Microorganisms are used in bioprocesses to mediate or catalyze the necessary reactions. This makes bioprocesses highly nonlinear and the governing mechanisms are complex. These complex governing mechanisms can be modeled by a metabolic network that comprises all interactions within the cells of the microbial population present in the bioprocess. The current state of the art in bioprocess control is model predictive control based on the use of macroscopic models, solely accounting for substrate, biomass, and product mass balances. These macroscopic models do not account for the underlying mechanisms governing the observed process behavior. Consequently, opportunities are missed to fully exploit the available process knowledge to operate the process in a more sustainable manner. In this article, a procedure is presented for metabolic network-based model predictive control. This procedure uses a combined moving horizon-model predictive control strategy to monitor the flux state and optimize the bioprocess under study. A CSTR bioreactor model has been combined with a small-scale metabolic network to illustrate the performance of the presented procedure.

1. Introduction

Biochemical processes have gained attention in the past decades in producing high added value products. Microorganisms play an important role in biochemical processes as these mediate or catalyze the reactions that are required to produce the products of interest. The mechanisms governing such biochemical processes are often complicated as a large number of reactions within the microbial cells and interactions between the microbial cells and the reactor medium take place. To summarize and model the mechanisms governing biochemical processes, metabolic network models have become an important tool. Metabolic networks comprise all interactions within the cell and between cell and environment by modeling the metabolites (i.e., the chemical compounds produced, consumed, and interacted with by the microorganisms) as the nodes and the reactions or fluxes between these different chemical compounds as the edges of a network. Mathematically, this network can be summarized in a stoichiometric matrix $\mathbf{S}$ in which an element $S_{ij}$ corresponds with the stoichiometric coefficient of the i-th metabolite in the j-th reaction or flux [1,2]. Genome-scale metabolic network model reconstructions are the result of intensive genomic studies involving numerous experiments. Hence, the result of such a study is a highly complex metabolic network with often a myriad of reactions and metabolites. Using genome-scale metabolic networks in an online control setting is to date infeasible as these models typically contain over 100 metabolites and reactions and hence a high number of differential states. Therefore, the current state of the art in model-based bioprocess optimization relies on macroscopic mass balance models.

This is one of the main reasons the current state of the art in model-based bioprocess optimization and control relies on macroscopic mass balance models, often only including biomass, substrate, and product mass balances. Macroscopic mass balance models make abstraction of the underlying mechanisms governing the studied bioprocess and therefore miss opportunities in achieving a better process monitoring and control. As there is no need to account for irrelevant parts of the metabolism for the process under study, genome-scale metabolic networks can be reduced to low complexity metabolic network models. The complexity of the metabolic network that should be used for a studied bioprocess depends on the level of detail required and the amount of experimental data or information available to identify the network. Hence, low- or medium-complexity metabolic networks can be used to account for the relevant intracellular mechanisms governing the bioprocess under study. Results from literature, e.g., [3], indicate that including such metabolic network models in bioprocess optimization results in an improved process performance compared with the use of an unstructured macroscopic bioprocess model. In addition, the intracellular fluxes can be estimated rather accurately from extracellular measurements using a low- or medium-complexity metabolic network as in e.g., [4,5,6,7,8], enabling an enhanced bioprocess monitoring by monitoring the flux state, hence the metabolic state of the microorganisms in the biochemical process.

The current state of the art in model-based bioprocess control is model predictive control (MPC) in which a process model is used to solve an open-loop dynamic optimization problem over a time interval, called the prediction horizon to optimize the future of a process. The optimal control inputs are computed until a new measurement becomes available. In biochemical processes, not all states can be measured and the states need to be estimated from the available measurements. For linear systems a Kalman filter can be used to retrieve an optimal state estimate [9]. For nonlinear systems extended Kalman filter [10] or the unscented Kalman filter [11] can be used. Alternatively, moving horizon estimation (MHE) can be used which solves a constrained optimization problem based on a set of measurements obtained at previous time points and also moves this estimation horizon when a new measurement is obtained [12,13,14]. MHE is quite similar to MPC in the sense that in both approaches a dynamic optimization problem is solved online and at every sampling instance the employed horizon is shifted. MHE looks at the past to estimate the states, while MPC looks at the future to compute controls.

In this article, a strategy is presented for online flux estimation and bioprocess control using low complexity metabolic network models. As indicated previously, one of the typical problems in online bioprocess control is the limited number of experimental measurements that are available. Therefore, a dynamic metabolic flux analysis (DMFA) model structure is combined with three kinetic models to estimate the intracellular fluxes: a constant, a linear, and a nonlinear model. The state estimation is performed based on MHE, while MPC is used for the bioprocess control. The presented strategy is applied to a continuous stirred tank (bio)reactor (CSTR) case study which uses a a low-complexity metabolic network.

The remainder of this article is structured as follows. Firstly, the Methods section consists of the description of the used model structure for metabolic reaction network-based modeling and the used kinetic model structures, followed by the introduction of the combined MHE/MPC strategy. Next, the case study is presented together with the results obtained after implementing the presented procedure. Finally, the main conclusions of this contribution are summarized.

2. Methods

2.1. Metabolic Reaction Network-Based Modeling

For a full overview of metabolic reaction network theory, the reader is referred to [1]. Here, only the final dynamic metabolic flux analysis model structure, which is used for the simulation and control parts of this work, is given [15]:

$$\begin{array}{cc}\hfill \frac{d{\mathbf{c}}_{\mathrm{ext}}\left(t\right)}{dt}& ={\mathbf{S}}_{\mathrm{ext}}·\mathbf{K}·\mathbf{u}\left(t\right)·c_{\mathrm{bio}}\left(t\right)\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{d{c}_{\mathrm{bio}}\left(t\right)}{dt}& ={\mathbf{s}}_{\mathrm{bio}}^{\top }·\mathbf{K}·\mathbf{u}\left(t\right)·c_{\mathrm{bio}}\left(t\right)\hfill \end{array}$$

(2)

with ${\mathbf{c}}_{\mathrm{ext}}$ being the ($m_{\mathrm{ext}}\times 1$) vector of concentrations of extracellular metabolites (mol/L), $c_{\mathrm{bio}}$ the scalar biomass concentration (gDW/L), ${\mathbf{S}}_{\mathrm{ext}}$ the $m_{\mathrm{ext}}$ rows of the stoichiometric matrix which correspond to the extracellular metabolites, and ${\mathbf{s}}_{\mathrm{bio}}^T$ the row of the stoichiometric matrix which corresponds to the biomass pseudo-metabolite. It is assumed that the production of biomass can be described by a pseudo-reaction/metabolic flux (a so-called biomass flux) as illustrated in [15]. This biomass concentration is used in Equations (1) and (2) to match dimensions of the RHS and the LHS as the (free) fluxes are expressed per gDW biomass and the metabolite concentrations are expressed in mol/L. $\mathbf{u}$ is the ($d\times 1$) vector of free fluxes (mol/gDW/h), and $\mathbf{K}$ is a suitable basis for the intracellular part of the stoichiometric matrix ${\mathbf{S}}_{\mathrm{int}}$. As the fluxes are parameterized as a function of the basis and the free fluxes $\mathbf{v}=\mathbf{K}\ast \mathbf{u}$, the choice of the basis is an important one as this determines the definition of the free fluxes and the complexity of the flux profiles and the number of parameters to be estimated. The authors of [15] demonstrated that three options exist to select a basis: (i) a fixed rational basis, (ii) a fixed orthonormal basis, or (iii) an optimal orthonormal basis. In this case, a fixed rational basis based on fluxes 1, 4, and 5 has been chosen. This matrix is of dimension ($n\times d$), where n is the total number of fluxes in the network. The number of free fluxes d is equal to $n-rank\left({\mathbf{S}}_{\mathrm{int}}\right)$. This formulation stems from the pseudo steady state that is assumed at the intracellular level, with respect to the extracellular dynamics.

This model is only fully defined when all free fluxes are expressed as a function of the concentrations and possibly environmental conditions like temperature and pH:

$$\begin{array}{cc}\hfill \mathbf{u}\left(t\right)& =\mathbf{f}({\mathbf{c}}_{\mathrm{ext}},T,pH,\dots ,\mathsf{\Phi })\hfill \end{array}$$

(3)

with $\mathsf{\Phi }$ being the parameter vector specific to the model structure which is used. As the fluxes are defined as specific fluxes, i.e., on a per-biomass basis, they do not depend on the biomass concentration anymore. Furthermore, in this contribution, only the dependence on the extracellular metabolite concentrations is considered.

2.2. Moving Horizon Estimation (MHE)

Moving horizon estimation (MHE) estimates states and parameters at a point in time $t_{i+H_e}$ from the measurements at this time point and the measurements at the previous $H_e$ time points. Hence, an estimation horizon of $H_e+1$ measurements is used for the state and parameter estimation. Once a new sample is taken at $t_{H_e}+1$, the estimation horizon moves one time step, meaning that the measurements at $t_{H_e}+1$ are used for the state and parameter estimation while the measurements related to the first point of the previous horizon (i.e., $t_i$) is discarded [12,13,14]. This is also illustrated in Figure 1.

Concretely, a dynamic combined state and parameter estimation problem (in essence, a dynamic optimization problem) is solved over a time interval $\left[t_i,t_{i+H_e}\right]$. In this article, the MHE formulation of [7] is used. This formulation is presented below for completeness. For more details on this formulation, the interested reader is referred to [7].

$$\begin{array}{c}\hfill \underset{{\mathbf{x}}_{H_e}^c,{\mathbf{w}}_i,\cdots ,{\mathbf{w}}_{i+H_e-1}}{\min }\sum _{j=i}^{i+H_e}{∥{\mathbf{m}}_j-\mathbf{y}\left(t_j\right)∥}_{\mathbf{V}}^2+\sum _{j=i}^{i+H_e-1}{∥{\mathbf{w}}_j∥}_{\mathbf{W}}^2+{∥{\mathbf{x}}_i^c-{\overline{\mathbf{x}}}_i^c∥}_{\overline{\mathbf{P}}}\mathbf{i}\end{array}$$

(4)

subject to:

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}\left(t\right)& ={\mathbf{S}}_{\mathrm{e}}·\mathbf{K}·\widehat{\mathbf{u}}(\mathbf{x},\mathbf{p})·{\mathbf{q}}_{\mathrm{bio}}^{\top }·\mathbf{x}\left(t\right)+\left({\mathbf{X}}_{\mathrm{in}}·\mathbf{r}-\mathbf{x}\left(t\right)\right)·D+{\omega }_{\mathrm{x}}\left(t\right)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \dot{\mathbf{p}}\left(t\right)& ={\omega }_{\mathrm{p}}\left(t\right)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \mathbf{x}\left(t_i\right)& ={\mathbf{x}}_i\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \mathbf{p}\left(t_i\right)& ={\mathbf{p}}_i\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \mathbf{0}& =\mathbf{z}\left(t\right)-{\mathbf{I}}_{\mathrm{irr}}·\mathbf{K}·\widehat{\mathbf{u}}(\mathbf{x},\mathbf{p})\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& =\mathbf{x}\left(t\right)\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& \ge \mathbf{0},\mathbf{z}\left(t\right)\ge \mathbf{0}\hfill \end{array}$$

(11)

with $\mathbf{x}\left(t\right)={\left[{\mathbf{c}}_{\mathrm{ext}}^{\top }\left(t\right)c_{\mathrm{bio}}\right]}^{\top }\left(t\right)$ being the state vector comprising both extracellular metabolites and biomass at t, $\widehat{\mathbf{u}}$ the estimated free fluxes, $\mathbf{p}$ the flux parameters, ${\mathbf{q}}_{\mathrm{bio}}^{\top }=\left[0\cdots 01\right]$ a row vector selecting the partition of the extracellular stoichiometric ${\mathbf{S}}_{\mathrm{e}}=\left[{\mathbf{S}}_{\mathrm{ext}}^{\top }{\mathbf{s}}_{\mathrm{bio}}\right]$ corresponding with biomass, ${\mathbf{X}}_{\mathrm{in}}\in \mathbb{R}4\times 2$ the matrix of concentrations at the bioreactor inlet for the different metabolites, and $\mathbf{r}\in {\mathbb{R}}^{2\times 1}$ the vector of control variables manipulating the feed composition sent to the reactor and D the dilution rate. In addition, the process noise on the states and parameters is denoted by ${\omega }_{\mathrm{x}}\left(t\right)$ and ${\omega }_{\mathrm{p}}\left(t\right)$, respectively. Algebraic states $\mathbf{z}\left(t\right)$ have been added for the description of the irreversible fluxes as denoted in Equation (9).

Note that ${\mathbf{x}}_{H_e}^c$, i.e., the combined state and parameter vector at $t_{i+H_e}$ and the process noise estimates ${\mathbf{w}}_i,\cdots ,{\mathbf{w}}_{i+H_e-1}$ are the optimization variables for the MHE problem described in Equations (4)–(11). Note that the arrival cost parameters ${\overline{\mathbf{x}}}_i^c$ and ${\overline{\mathbf{P}}}_{\mathbf{i}}$ are updated via an unscented Kalman filter [16].

2.3. Unscented Kalman Filter (UKF)

The unscented Kalman filter (UKF) is an extension of the unscented transformation to the Kalman filter. The unscented transform is a way to calculate the mean and variance of the nonlinear transformation of a random variable that follows a symmetric unimodal distribution. In this article, the states and parameters are the nonlinear transformation of the random variable (i.e., the initial conditions of the states and parameters). In UKF, the sigma points, a fixed number of deterministically chosen sampling points from the distribution of the random variable and are then propagated through the nonlinear function. The mean and covariance of this propagation are subsequently approximated. For the implementation of the UKF in this work, the sigma points are chosen based on recommendations from [11]. Different strategies can be used to select these sigma points and recently a comparative study has been presented in [17].

2.4. Model Predictive Control (MPC)

Model predictive control (MPC) solves an open-loop dynamic optimization problem over a time interval $\left[t_k,t_{k+H_p}\right]$, with $H_p$ as the prediction horizon to optimize the future of a process. The novel obtained measurements and state and parameter estimates at $t_k$ are used as initial conditions for the open-loop dynamic optimization problem that is solved at every iteration. This principle is illustrated in Figure 2: based on the state estimate or measurement at $t_k$ a dynamic optimization problem is solved over the time interval $\left[t_k,t_{k+H_p}\right]$ (a). The solution is applied to the process until a new measurement at $t_{k+1}$ is obtained and then a new dynamic optimization problem is solved over $\left[t_{k+1},t_{k+H_p+1}\right]$ (b).

MPC has become a widespread control strategy and numerous successful applications to industrial problems have been reported, e.g., [18,19]. MPC typically refers to linear MPC in which a quadratic objective function is minimized using a linear process model and subject to linear constraints. In the last decades, nonlinear MPC (NMPC) approaches using nonlinear process models, constraints, and objective functions have become more popular.

In this article the following (N)MPC formulation for $t\in \left[t_k,t_{k+H_p}\right]$ is used, based on [20]:

$$\begin{array}{c}\hfill \underset{\mathbf{x}(·),\mathbf{u}(·)}{\min }{\int }_{t_k}^{t_k+H_p}L(\mathbf{x}\left(\tau \right),\mathbf{u}\left(\tau \right),)d\tau +M(t_k+H_p)\end{array}$$

(12)

subject to:

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}\left(t\right)& ={\mathbf{S}}_{\mathrm{e}}·\mathbf{K}·\mathbf{u}(\mathbf{x},\mathbf{p})·{\mathbf{q}}_{\mathrm{bio}}^{\top }·\mathbf{x}\left(t\right)+\left({\mathbf{X}}_{\mathrm{in}}·\mathbf{r}-\mathbf{x}\left(t\right)\right)·D\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \dot{\mathbf{p}}\left(t\right)& =\mathbf{0}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \mathbf{x}\left(t_k\right)& ={\widehat{\mathbf{x}}}_k\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \mathbf{p}\left(t_k\right)& ={\widehat{\mathbf{p}}}_k\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \mathbf{0}& =\mathbf{z}\left(t\right)-{\mathbf{I}}_{\mathrm{irr}}·\mathbf{K}·\mathbf{u}(\mathbf{x},\mathbf{p})\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& =\mathbf{x}\left(t\right)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \mathbf{x}\left(t\right)& \ge \mathbf{0},\mathbf{z}\left(t\right)\ge \mathbf{0}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathbf{0}& \ge \mathbf{c}(\mathbf{x},\mathbf{r},\mathbf{p},t).\hfill \end{array}$$

(20)

with $\widehat{\mathbf{x}}\left(t_k\right)$ being the measured or estimated state values at $t_i$, $\mathbf{c}(\mathbf{x},\mathbf{r},\mathbf{p},\tau )$, constraint functions on the systemare not comprised by the previous constraints, $\tau$ is the time used in solving the dynamic optimization problems in the NMPC routine, and $\mathbf{x}$ and $\mathbf{r}$ the states and controls computed by solving the dynamic optimization problem on the interval $\left[t_k,t_k+H_p\right]$. Note that for the other symbols the same notation is used as in the MHE problem (Equations (4)–(11)). As MPC requires that the entire state vector is known, the states and parameters estimates at $t_k$ are denoted by ${\widehat{\mathbf{x}}}_k$ and ${\widehat{\mathbf{p}}}_k$.

2.5. Combined MHE/MPC Strategy

In this article, a strategy combining MHE and MPC is used for the model-based predictive control of bioprocesses exploiting metabolic network information. The moving-horizon strategy used in this work is thus different from regular implementations because of the black-box nature of the flux models which are used. The proposed strategy is outlined in Figure 3.

The algorithm starts with taking a sample and use an unscented Kalman filter (UKF) to estimate the combined state and parameter vector and the combined state and parameter variance-covariance matrix. Subsequently, it is checked whether a sufficient number of measurements is available to perform MHE, i.e., $H_e+1$ measurements need to be available or equivalently: $t_k\ge t_{H_e}$. If $t_k<t_{H_e}$, the state and parameter estimates from UKF are used for the MPC step. Otherwise, the state and parameter estimates from the MHE procedure are used for the MPC step. The control action computed during the MPC step is then applied to the process until the next sample is obtained. Once the end time $t_f$ is reached, the algorithm is stopped.

This combined MHE-MPC strategy has been implemented in the Pomodoro toolkit [21], which is available at https://cit.kuleuven.be/biotec/software/pomodoro (accessed on 11 October 2021). Pomodoro is a toolkit that can be used for solving dynamic (multi-objective) optimization problems and for model-based control and estimation. Pomodoro uses CasADi [22] as a symbolic framework for the problem formulation and derivative computation. Orthogonal collocation is used as default to discretize the dynamic optimization problems. A piecewise constant discretization is used for the controls and cubic Lagrange polynomials with collocation points at the Radau roots are used for the states [23]. The resulting NLPs are solved with IPOPT [24].

3. Results and Discussion

The proposed combined MHE-MPC strategy for metabolic network model-based predictive bioprocess control has been implemented on a CSTR bioreactor using a low-complexity metabolic network. The general setup for the implementation of this case study is shown in Figure 4. Two different models are used: (i) the simulation model, which gives the process dynamics, and to which output noise is added to get the process measurements, and (ii) the estimation/control model, which is an approximation since the process model is not known. It is important to note that only the flux kinetics are assumed to be unknown in the estimation/control model. The metabolic reaction network in the estimation/control model is for this case study the same as the one used in the simulation model. Additionally, the general model structure, i.e., the transport terms, etc., are the same for both the simulation and the estimation/control model.

3.1. Case Study Description

The bioprocess in this case study is a continuous bioreactor with two substrates and two products. The composition of the medium that is fed to the bioreactor is controlled, but the flow rate of this feed is always kept constant, as well as the flow rate of the spent medium that is removed from the bioreactor. Since both these flow rates are equal, the volume of the bioreactor is kept constant. Mathematically, the dynamic model consists of (i) the reaction term as defined in Equations (21) and (22), and (ii) transport terms because of the supply of substrate to the reactor and removal of medium from the reactor. The resulting model is the following:

$$\begin{array}{cc}\hfill \frac{d{\mathbf{c}}_{\mathrm{ext}}\left(t\right)}{dt}& ={\mathbf{S}}_{\mathrm{ext}}·\mathbf{K}·\mathbf{u}\left(t\right)·c_{\mathrm{bio}}\left(t\right)+({\mathbf{C}}_{\mathrm{in}}·\mathbf{w}-{\mathbf{c}}_{\mathrm{ext}})·D\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill \frac{d{c}_{\mathrm{bio}}\left(t\right)}{dt}& ={\mathbf{s}}_{\mathrm{bio}}^T·\mathbf{K}·\mathbf{u}\left(t\right)·c_{\mathrm{bio}}\left(t\right)-c_{\mathrm{bio}}·D\hfill \end{array}$$

(22)

with F as the flow rate (L/h) in and out of the bioreactor and V the volume (L) of the bioreactor, making $D=\frac{F}{V}$ the dilution rate (1/h). In this article, the dilution rate is fixed to 0.1 h${}^{-1}$. The ($3\times 3$) substrate composition matrix ${\mathbf{C}}_{\mathrm{in}}$ is made up of three columns that express the composition of each of the three feed mixtures, and the ($3\times 1$) control vector $\mathbf{w}$ contains the weights of these feeds in the final mixture which is sent to the bioreactor. To make sure the total flow rate to the bioreactor is equal to the outgoing flow rate, these controls should always sum to one:

$$\begin{array}{c}\hfill \sum _{i=1}^3{w}_i=1\end{array}$$

(23)

This way, there are only two independent controls in the system. The numerical values for the model parameters are given in

Table 1. Model parameter values.

Model Parameter	Numerical Value
F	$0.1$	L/h
V	$1.0$	L
${\mathbf{C}}_{\mathrm{in}}$	$\left[\begin{array}{ccc}6.0& 0.0& 0.0\\ 0.0& 9.0& 0.0\\ 0.0& 0.0& 0.0\end{array}\right]$	mol/L

.

The metabolic reaction network, simulation model, and the estimation/control model are briefly discussed below.

3.1.1. Metabolic Reaction Network

The metabolic reaction network that is used in the case study is shown in Figure 5 on the left, along with the corresponding ${\mathbf{S}}_{\mathrm{int}}$, ${\mathbf{S}}_{\mathrm{e}}$, ${\mathbf{I}}_{\mathrm{irr}}$, and the selected null basis $\mathbf{K}$ (note that ${\mathbf{K}}^{\top }$ is depicted in Figure 5) matrices. It consists of threee extracellular metabolites (A${}_{\mathrm{ext}}$, E${}_{\mathrm{ext}}$, F${}_{\mathrm{ext}}$) and biomass (X), four intracellular metabolites (A, B, C, D), and seven fluxes ($v_1$, $v_2$, $v_3$, $v_4$, $v_5$, $v_6$ and $v_{\mathrm{bio}}$). Thus, the number of free fluxes is three.

3.1.2. Simulation Flux Model

For the simulation of the measurements, these were chosen as flux 1, 4 and 5. Measurements were simulated by choosing reference profiles for these three fluxes, and simulating the states using the dynamic system.

$$\begin{array}{cc}\hfill u_{1,\mathrm{ref}}& =\frac{c_{\mathrm{Aext}}}{1.5+c_{\mathrm{Aext}}}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill u_{4,\mathrm{ref}}& =0.2·\frac{c_{\mathrm{Eext}}}{3+c_{\mathrm{Eext}}}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill u_{5,\mathrm{ref}}& =\frac{1}{1+c_{\mathrm{Fext}}}\hfill \end{array}$$

(26)

This simulation flux model is used to generate measurements over a time range of 140 h, measuring ${\mathrm{A}}_{\mathrm{ext}}$, ${\mathrm{E}}_{\mathrm{ext}}$, ${\mathrm{F}}_{\mathrm{ext}}$. Additive, independent, Gaussian measurement noise was added to the output of the simulation model with 0.01 as standard deviation. The initial conditions were set to ${\mathbf{x}}_0={\left[0.5760,3.5527,2.4976,0.8736\right]}^{\top }$. The process noise was also considered to be of the same type but with a standard deviation of 0.001. The standard deviations to be used in the arrival cost parameters ${\overline{\mathbf{P}}}_0$ is taken equal to the measurement variance for the states and set to 1.0 for the parameters. Different sampling frequencies are studied and the same holds for the prediction and estimation horizons $H_p$ and $H_e$, which are chosen to be equal in this article. This is not necessarily the case in practice.

3.1.3. Estimation/Control Flux Model

Three different kinetic flux model structures are studied as estimation/control flux model, a constant, a linear, and a non-linear model.

The constant model considers the fluxes to be constant:

$$\begin{array}{cc}\hfill \mathbf{u}\left(t\right)& ={\mathbf{p}}_{\mathrm{u}}\hfill \end{array}$$

(27)

with ${\mathbf{p}}_{\mathrm{u}}\in {\mathbb{R}}^{d\times 1}$, thus three parameters to be identified.

In the linear model, the fluxes are just a linear function of the extracellular concentrations:

$$\begin{array}{cc}\hfill \mathbf{u}\left(t\right)& =\mathbf{C}·{\mathbf{c}}_{\mathrm{ext}}\left(t\right)\hfill \end{array}$$

(28)

with $\mathbf{C}$ as the ($d\times m_{\mathrm{ext}}$) matrix of parameters that define the linear relationship between fluxes and metabolite concentrations. This results in nine parameters to be identified.

The non-linear kinetics model that is used in this work is the general biological reaction kinetics model of [25]. This model can describe both positive and negative effects of components on reaction rates, and uses only one parameter per component to describe both effects, and one maximum rate parameter per flux, i.e., $\left(m_{\mathrm{ext}}+1\right)\times d$ parameters:

$$\begin{array}{cc}\hfill u_i\left(t\right)& =u_{\max ,i}·\prod _{j=1}^{m_{\mathrm{ext}}}{\alpha }_{ij}\left(c_{\mathrm{ext},j}\left(t\right)\right)\hfill \end{array}$$

(29)

with

$${\alpha }_{ij}\left(c_{\mathrm{ext},j}\left(t\right)\right)=\left\{\begin{array}{cc}\frac{c_{\mathrm{ext},j}}{c_{\mathrm{ext},j}+D_{ij}^2}\hfill & \mathrm{if}{D}_{ij}\ge 0\hfill \\ \frac{1}{1+c_{\mathrm{ext},j}·D_{ij}^2}\hfill & \mathrm{if}{D}_{ij}<0\hfill \end{array}\right.$$

(30)

In these equations, ${\mathbf{u}}_{\max }$ is the ($d\times 1$) vector of maximum rate parameters, and $\mathbf{D}$ is the ($d\times m_{\mathrm{ext}}$) matrix of modulation parameters. If these parameters are positive, they describe a positive activation/saturation effect of the concentrations on the fluxes, if they are negative, an inhibition effect is represented. Hence, 12 parameters need to be identified i.e., 4 parameters per free flux.

3.2. Numerical Results

Initially, the CSTR bioreactor is studied over a period of 140 h with a sampling rate of 1 h, resulting in 140 simulation intervals. During the first 10 h, a constant ratio of 50% A${}_{\mathrm{ext}}$ and E${}_{\mathrm{ext}}$ is fed to the CSTR reactor. A horizon length of 4 h has been selected for both the estimation horizon as the prediction horizon (i.e., $H=H_e=H_p=4$ h, corresponding to 20 intervals). This means that after 4 h, the first MPC control action is applied to the process. Moreover, the set point ${\mathbf{x}}_{\mathrm{SP}}$ for the MPC is defined as follows:

$$\begin{array}{c}\hfill {\mathbf{x}}_{\mathrm{SP}}=\left\{\begin{array}{c}{\left[0.5448,1.3355,3.3852,1.4845\right]}^{\top }\mathrm{if}40\mathrm{h}\le t_k\le 80\mathrm{h}\hfill \\ {\left[0.5760,3.5527,2.4976,0.8736\right]}^{\top }\mathrm{else}\hfill \end{array}\right.\end{array}$$

(31)

3.2.1. State Estimates

In Figure 6 the simulated state profiles are depicted together with the estimates obtained when using the estimation/control model. The full lines indicate the simulations while the markers indicate the estimates. For all different models the proposed strategy results in accurate state estimates. This complies with the results presented in [7]. For an extensive discussion focusing only on state estimation the interested reader is referred to [7]. Recall that the arrival cost matrices are computed with an unscented Kalman filter in this article, contrary to the arrival cost computation in [7]. In [7] the arrival cost is approximated by a quadratic cost and linearizing the nonlinear functions to obtain a linear least squares problem that can be solved analytically and results in an analytical expression for the arrival cost [14].

A notable difference can be observed when comparing the state profiles for the different flux models. The cause for this difference can be twofold: (i) the linear and Haag models are more similar to the true flux model than the constant flux model hence resulting in more similar flux profiles, and (ii) the control profiles for the true, linear, and nonlinear flux models are more similar, compared with the constant flux model. This is investigated in more depth below.

3.2.2. MPC Tracking Performance

In Figure 7 the extracellular metabolite concentration profiles (i.e., the state profiles) are depicted for the different methods, together with the set point that is used for the different concentrations during the MPC procedure. For the extracellular metabolites that are consumed, i.e., A${}_{\mathrm{ext}}$ and E${}_{\mathrm{ext}}$, the set points are better tracked by the true flux model and the constant flux model than by the linear and Haag models. The extracellular metabolites that are excreted (or produced) by the cells are better approximated with the true, linear, and Haag flux models. These observations are made from Figure 7.

In order to quantify the difference in tracking performance between the MPC control strategies obtained with the different flux models, the unscaled absolute values of the absolute tracking errors for each state are depicted in

Table 2. Absolute values of the absolute tracking errors, i.e., difference between state and set point.

	True	Constant	Linear	Haag
A${}_{\mathrm{ext}}$	13.44	5.69	19.48	24.38
E${}_{\mathrm{ext}}$	148.52	191.61	172.80	234.86
F${}_{\mathrm{ext}}$	14.99	57.87	14.66	15.06
Biomass	6.19	41.64	6.96	6.62

. From this table it becomes clear that the constant model has the lowest tracking error for the first state A${}_{\mathrm{ext}}$. For E${}_{\mathrm{ext}}$ the tracking error is the lowest for the true flux model, while the constant and linear models perform slightly worse. The Haag model performs the worst for E${}_{\mathrm{ext}}$. F${}_{\mathrm{ext}}$ is comparably well tracked by the true model, linear model, and Haag model. The constant flux model has a poor tracking for F${}_{\mathrm{ext}}$. The biomass concentration is well tracked by the true model, linear model, and Haag model, while the constant model is again performing poorly. These observations confirm the observations made from Figure 7.

For completeness, the control profiles are depicted in Figure 8. The control profiles computed with the constant flux model are flatter than the control profiles computed with the other flux models. The difference between these control profiles confirms the tracking performance discussed above.

3.2.3. Influence of Horizon Length

Same Horizon and Estimation Horizon

Firstly, the influence of the horizon length on the tracking performance is studied by selecting the same horizon for estimation and prediction (control): $H=H_e=H_p$. Four different horizons have been simulated, assuming that a sample is taken every hour: $H=2\mathrm{h}$, $H=4\mathrm{h}$, $H=12\mathrm{h}$, $H=24\mathrm{h}$. The absolute values of the absolute tracking error for all states are depicted in Figure 9.

Globally, it can be observed that the MPC control performance obtained with the true flux model and the linear flux model are more similar for a small horizon length, i.e., H = 1 h. However, no clear global trends can be observed for the increase of the horizon length. As both the estimation and prediction horizon have been modified, it is difficult to interpret these graphs. Therefore the effect of only changing the prediction horizon or solely changing the estimation horizon is addressed in the next paragraphs.

Influence of Only Changing the Prediction Horizon

In this paragraph the estimation horizon is set to $H_e=4\mathrm{h}$ and the prediction horizon has been varied as in the previous paragraph. Four different horizons have been simulated, assuming that a sample is taken every hour: $H_p=2\mathrm{h}$, $H_p=4\mathrm{h}$ (4 h), $H_p=12\mathrm{h}$, $H_p=24\mathrm{h}$. In Figure 10, the absolute values of the absolute tracking error are displayed.

Increasing the prediction horizon results in a slight reduction of the tracking error of all states using the true flux model. For the other flux models, such clear observations cannot be made.

Influence of Only Changing the Estimation Horizon

Contrary to the previous paragraph, the estimation horizon is varied and the prediction horizon is kept constant at 4 h. Following estimation horizons have been studied for a sampling rate of 1 h: $H_e=2\mathrm{h}$, $H_e=4\mathrm{h}$ (4 h), $H_e=12\mathrm{h}$, $H_e=24\mathrm{h}$. In Figure 11, the absolute values of the absolute tracking error are displayed.

The estimation horizon has a small effect on the MPC tracking performance, except for the linear flux model, in which the estimation horizon of 24 h clearly results in a higher on the MPC tracking error for all extracellular metabolites except F${}_{\mathrm{ext}}$.

3.2.4. Influence of Sampling Frequency

The influence of the sampling frequency on the tracking performance has been studied for an estimation and horizon length of 4 h ($H=H_e=H_p=4$ h). Four sampling frequency scenarios have been studied: 12 min, 30 min, 1 h, and 2 h. The MPC tracking error divided by the number of sampling points has been depicted in Figure 12.

It is expected that a higher sampling frequency will result in better tracking. This corresponds with lower sampling times $dt$. In Figure 12 this is only observed for the true model. For the linear model this is the case for all states, except F${}_{\mathrm{ext}}$. In the Haag model, this is not observed for E${}_{\mathrm{ext}}$. A lower sampling frequency only results in a better tracking for the constant flux model for A${}_{\mathrm{ext}}$.

4. Conclusions

One of the current limitations in online bioprocess control is the use of macroscopic models (that only account for substrate, biomass, and product mass balances), as the available process knowledge on the relevant underlying intracellular mechanisms is not fully exploited and as such opportunities are missed towards an optimal bioprocess operation. In this article, a procedure for metabolic network-based model predictive control has been presented, by combining online flux estimation and bioprocess control using metabolic network models. This procedure firstly uses moving horizon estimation (MHE) to estimate the states and fluxes from the available measurements and subsequently solves a model predictive control (MPC) problem to determine the optimal control strategy for the studied bioprocess. To illustrate the performance of this procedure, a case study combining a CSTR bioreactor model with a metabolic network model has been implemented. Three different model structures (a constant, a linear, and a nonlinear model) have been compared, clearly indicating that the linear flux model outperforms the nonlinear Haag flux model and the constant flux model. The true model is, as expected, outperforming the other models. In addition, the influence of the horizon lengths and sampling frequency on the MPC tracking performance has been studied. While an increase of the prediction horizon length was expected to reduce the tracking error, such observation could only be made for the true model. An increase in sampling frequency (or reduction in sampling time) was expected to reduce the tracking error, but this observation could only be made for the true model and for the linear model.

In future work, the presented procedure will be applied to the online control of bioprocesses involving large-scale metabolic networks. In addition, uncertainty will be accounted for during the MPC routine. Furthermore, the presented metabolic network-based model predictive control procedure can also be combined with high-throughput biological datasets [26,27] and advanced metabolomics techniques as e.g., low-frequency NMR, to improve the metabolic model predictions and process monitoring capabilities [8,27]. One of the current limitations of metabolic network models and methods as metabolic flux analysis is the inherent pseudo-steady state assumption. This assumption is not always valid, certainly not at timescales of multiple scales due to shifting from exponential growth phase to stationary phase. The authors of [28] proposed a hybrid modeling approach, combining machine learning with metabolic network modeling which provided accurate predictions for amino acid consumption rates in Chinese Hamster Ovary Cells. Such hybrid models could also be used in the proposed metabolic network-based MPC procedure.