Parameter Estimation for a Kinetic Model of a Cellular System Using Model Order Reduction Method

Abstract

Order reduction methods are important tools for systems engineering and can be used, for example, for parameter estimation of kinetic models for systems biology applications. In particular, the Proper Orthogonal Decomposition (POD) method produces a reduced-order model of a system that is used for solving inverse problems (parameter estimation). POD is an intrusive model order reduction method that is aimed to obtain a lower-dimensional system for a high-dimensional system while preserving the main features of the original system. We use a singular value decomposition (SVD) to compute a reduced basis as it is usually numerically more robust to compute the singular values of the snapshot matrix instead of the eigenvalues of the corresponding correlation matrix. The reduced basis functions are then used to construct a data-fitting function that fits a known experimental data set of system substance concentrations. The method is applied to calibrate a kinetic model of carbon catabolite repression (CCR) in Escherichia coli, where the regulatory mechanisms on the molecular side are well understood and experimental data for a number of state variables is available. In particular, we show that the method can be used to estimate the uptake rate constants and other kinetic parameters of the CCR model.

1. Introduction

In systems biology, mathematical models of cellular systems are essential for predicting and optimizing cells behavior in culture  [1,2,3]. One category of these mathematical models is a continuous model, which describes the change of substance concentrations over time via a system of differential equations. Mathematical models are capable of understanding the full picture of biological systems, but finding the right parameter values (e.g., reaction rate constants) is a serious challenge in systems biology.

Parameter estimation aiming to determine the parameter values for a mathematical model so that the dynamical system evolves in a way such that the system characteristics obtained from experimental observations are properly matched, see [4,5,6]. Mathematically, the parameter estimation problem corresponds to an inverse problem [7], i.e., one seeks to find the unknown model inputs (e.g., parameter values) using the measured system states and other available information [8]. Different algorithms of computational optimization can be used for estimating unknown parameters of a dynamical system [9,10], e.g., Gauss–Newton Method [11], least square regression [12,13], or maximum likelihood estimation [14]. In addition, parameter estimation was recently discussed using machine learning algorithms, for instance, convolutional neural networks [15].

Proper orthogonal decomposition (POD) is a model order reduction (MOR) technique is aiming to project high-dimensional system onto a lower-dimensional subspace whilst retaining the most important features of the original system [16,17]. The POD method, also known as Karhunen–Loève Expansion (KLE) [18,19], is considered an optimal linear method since it minimizes the squared distance between the original and reduced model. Recently, the POD has received much attention for analyzing complex physical systems [20,21,22]. Moreover, it has been applied to different dynamical systems in biology, e.g., the large-scale kinetic model of the metabolic network of Escherichia coli [3,23] and oscillating biological network models [24].

The POD method generates a reduced-order model that is a linear combination of basis elements. This basis elements are computed to capture the essential features of the system [25]. Thus, the POD method is well suited in optimal control and for solving inverse problem of dynamical systems, like the parameter estimation problem. Moreover, the POD method takes an optimally large set of parameter values and then reduces the set to a few estimates to be tested, so it is considered less computationally intense [24]. Often, the POD method is used for parameter estimation of partial differential equations (PDEs), e.g., in [26] it is used to compute a reduced-order model for the bidomain equations of cardiac electrophysiology which is utilized in an inverse problem solved using an evolutionary algorithm. Also in [27], the POD method has been applied to estimate scalar parameters in an elliptic PDE. In [24], the POD method is used for parameter estimation of ordinary differential equation (ODE) systems for stable oscillating biological networks. The study shows that the POD method is more accurate than the spline method for a stable oscillating network.

The aim of this study is to discuss parameter estimation using POD method for an ODE system representing the kinetic model of a cellular system. Specifically, the method is applied to the carbon catabolite repression (CCR) mechanism in Escherichia coli [28]. The POD method provides a reduced-order model that is a linear combination of reduced basis functions. We use a singular value decomposition (SVD) of the snapshot matrix to compute the reduced basis, while in [24,25,27], the reduced basis is computed via the eigenvalues and eigenvectors of the correlation matrix. The snapshot matrix is generated from numerical simulations of the kinetic model for different parameter values and with certain initial conditions. The reduced basis functions are then used for constructing a data-fitting function that fits measurement data sets of dynamical system substance concentrations. Non-linear least square fitting is considered the most common approach to estimate parameters for non-linear systems [9,10]. Therefore, we are aiming to perform parameter estimation for the CCR kinetic model using POD and the Least square method (LSQ) for comparison. As we will show in this study, the results of the CCR kinetic model with the obtained estimated parameter values using the POD method fit well with the experimental data.

The remainder of this paper is organized in the following manner. In Section 2, we discuss the parameter estimation procedure for non-linear ordinary differential equations using the POD method, singular value decomposition, and constructing a data-fitting function for the observable components. In Section 3, the approaches are applied to the carbon catabolite repression network of E. coli. Finally, we end with some concluding remarks in Section 4.

2. Parameter Estimation for ODEs Using the POD Method

We consider a system of ordinary differential equations given by

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =f(x\left(t\right),\Theta ),x\left(t_0\right)=x_0\in {\mathbb{R}}^{n_m}\mathrm{on}\mathbb{I}=[t_0,t_{end}],\hfill \end{array}$$

(1)

where $x\left(t\right)\in {\mathbb{R}}^{n_m}$ is the state vector (in our case the vector of concentration of $n_m$ species) and $\Theta \in {\mathbb{R}}^{n_{\theta }}$ is a vector of parameters that contains all unknown constants determining the dynamics, e.g., kinetic parameters. While the system may depend on a number of parameters $\Theta =[{\theta }_1,\cdots ,{\theta }_{n_{\theta }}]$, at first, we estimate one parameter, $n_{\Theta }=1$.

Assume that we have some measurement data from experiments that are given in the following matrix

$$Y=\left[\begin{array}{ccc}{y}_1\left({\tau }_1\right)& \cdots & y_1\left({\tau }_{n_{\tau }}\right)\\ & ⋮& \\ y_{n_d}\left({\tau }_1\right)& \cdots & y_{n_d}\left({\tau }_{n_{\tau }}\right)\end{array}\right]\in {\mathbb{R}}^{n_d\times n_{\tau }},$$

containing measurements for $n_d$ components at $n_{\tau }$ measurement time points ${\tau }_j$ in the interval $\mathbb{I}$. Our final goal is to fit the parameters of the mathematical model to the experimental data, we assume a suitable ordering of x for the sake of implementation simplicity i.e.,  $x_i\left({\tau }_j\right)\approx y_i\left({\tau }_j\right)$ for $i=1,\cdots ,n_d<n_m$ and $j=1,\cdots ,n_{\tau }$ with sufficient high accuracy.

2.1. Creating the Snapshot Matrix

We perform $n_s$ simulations of the model using different sets of parameter values ${\Theta }_1,\cdots {\Theta }_{n_s}$ equidistantly distributed in $[lb,ub]$, where $lb$ and $ub$ are the minimum and maximum for all parameter values. A Latin hypercube sampling (LHS) is used for generating the parameter space [29,30]. In each of the $n_s$ simulations, we evaluate at $n_t$ simulation time points. Note that $n_t$ can be much larger than $n_{\tau }$. In this way, we obtain the snapshot matrix

$$X=\left[\begin{array}{ccc}\tilde{x}(t_1;{\Theta }_1)& \cdots & \tilde{x}(t_{n_t};{\Theta }_1)\\ \tilde{x}(t_1;{\Theta }_2)& \cdots & \tilde{x}(t_{n_t};{\Theta }_2)\\ & ⋮& \\ \tilde{x}(t_1;{\Theta }_{n_s})& \cdots & \tilde{x}(t_{n_t};{\Theta }_{n_s})\end{array}\right]\in {\mathbb{R}}^{n_m·n_s\times n_t},$$

where $\tilde{x}(t_i;{\Theta }_j)$ denotes a numerical approximation to the solution $x(t_i;{\Theta }_j)$ of (1) at the time point $t_i$ for parameter value ${\Theta }_j$.

We can regard $X=\left[\begin{array}{ccc}\overline{x}\left(t_1\right)& \cdots & \overline{x}\left(t_{n_t}\right)\end{array}\right]$ as the snapshot matrix taken for the enlarged system

$$\begin{array}{c}\hfill \dot{\overline{x}}=\overline{f}\left(\overline{x}\right),\overline{x}\left(t_0\right)={\overline{x}}_0\in {\mathbb{R}}^{n_m·n_s},\end{array}$$

(2)

where $\overline{x}\left(t\right):=\left[\begin{array}{c}x(t;{\Theta }_1)\\ x(t;{\Theta }_2)\\ ⋮\\ x(t;{\Theta }_{n_s})\end{array}\right]\in {\mathbb{R}}^{n_m·n_s}$, with each $x(t;{\Theta }_i)\in {\mathbb{R}}^{n_m}$ is a solution of (1), and  $\overline{f}\left(\overline{x}\right)=\left[\begin{array}{c}f(x,{\Theta }_1)\\ f(x,{\Theta }_2)\\ ⋮\\ f(x,{\Theta }_{n_s})\end{array}\right]$, as well as ${\overline{x}}_0=\left[\begin{array}{c}{x}_0\\ x_0\\ ⋮\\ x_0\end{array}\right]$.

2.2. Computing the Reduced Basis

To obtain a reduced basis, one can use the eigenvalues of the correlation matrix $C=X^{T}X$ as in [24,25]. Since the eigenvalues ${\lambda }_j$ of C, and the singular values ${\sigma }_j$ of X are related as ${\lambda }_j={\sigma }_j^2$, we use a singular value decomposition (SVD) to compute the singular values of X instead.

If the matrix X has the condition number $\kappa \left(X\right)$, then the condition number of C is $\kappa \left(C\right)=\kappa {\left(X\right)}^2$, such that for problems with high condition number the computation of eigenvalues is much less robust that the computation of singular values, see e.g., [31].

Let $d=rank\left(X\right)\le min(n_m·n_s,n_t)$. The SVD guarantees the existence of real numbers ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_d>0$ and orthogonal matrices $V\in {\mathbb{R}}^{n_m·n_s\times n_m·n_s}$ with columns ${\left\{v^i\right\}}_{i=1}^{n_m·n_s}$ and $W\in {\mathbb{R}}^{n_t\times n_t}$ with columns ${\left\{w^i\right\}}_{i=1}^{n_t}$, such that

$$\begin{array}{c}\hfill V^{T}XW=\left[\begin{array}{cc}D& 0\\ 0& 0\end{array}\right]=:\Sigma \in {\mathbb{R}}^{n_m.n_s\times n_t}\end{array}$$

(3)

where $D=diag({\sigma }_1,\cdots ,{\sigma }_d)\in {\mathbb{R}}^{d\times d}$ and the zeros denote matrices of appropriate dimensions. From (3), we obtain

$$X=V\Sigma W^T,$$

or similarly, the vectors ${\left\{v^i\right\}}_{i=1}^d$ and ${\left\{w^i\right\}}_{i=1}^d$ satisfy

$$X{w}^i={\sigma }_i{v}^i\iff v^i=\frac{1}{{\sigma }_i}X{w}^i,\mathrm{for}i=1,\cdots ,d.$$

To obtain a reduced basis we can now cut off singular values with ${\sigma }_i<\epsilon$, where $\epsilon$ is chosen such that a smaller number of basis elements $n_r<<n_m·n_s$ is sufficient to capture the main features of the solution of (2). By setting $V_r=[v^1,\cdots ,v^{n_r}]\in {\mathbb{R}}^{n_m·n_s\times n_r}$ (the so-called POD modesor POD basis vectors), then a reduced order model (ROM) for (2) is given by

$$\begin{array}{c}\hfill {\dot{\overline{x}}}_r={\overline{f}}_r(t,{\overline{x}}_r),{\overline{x}}_r\left(t_0\right)=V_r^T{\overline{x}}_0\in {\mathbb{R}}^{n_r},\end{array}$$

(4)

where $\overline{x}\approx V_r{\overline{x}}_r$, i.e., $V_r{\overline{x}}_r$ serves as approximation for $\overline{x}$, and ${\overline{f}}_r(t,{\overline{x}}_r):=V_r^T\overline{f}(t,V_r{\overline{x}}_r)$.

It is well-known that $span\{v^1,\cdots ,v^{n_r}\}$ is the best approximation of $Range\left(X\right)$ in the sense that minimizes the two-norm of the approximation error:

$$\sum _{j=1}^{n_t}{\parallel \overline{x}\left(t_j\right)-V_r{V}_r^T\overline{x}\left(t_j\right)\parallel }_2^2=\sum _{i=n_r+1}^d{\sigma }_i^2,$$

see  [32].

Usually, the goal is to choose $n_r$ small enough while the relative information content  [33] of the basis for the $n_r$-dimensional subspace, defined by

$$I\left(n_r\right)=\frac{{\sum }_{i=1}^{n_r}{\sigma }_i}{{\sum }_{i=1}^d{\sigma }_i},$$

is near to one. If the $n_r$-dimensional subspace should contain a percentage p of the information contained in the full dimensional space ${\mathbb{R}}^{n_m·n_s}$, then one should choose $n_r$ such that

$$n_r=\mathrm{argmin}\left\{I\left(n_r\right)|I\left(n_r\right)\ge \frac{p}{100}\right\}.$$

Each data vector ${\overline{x}}^i\in {\mathbb{R}}^{n_m·n_s}$, i.e., the i-th column of X, can be written as

$${\overline{x}}^i=\sum _{j=1}^d{b}_{ij}{v}^j,\mathrm{with}{b}_{ij}=\langle v^j,{\overline{x}}^i\rangle ={\left(v^j\right)}^T·{\overline{x}}^i,$$

(5)

where $\langle .,.\rangle$ denotes the canonical inner product in ${\mathbb{R}}^{n_m·n_s}$.

Thus, the state vector can be approximated in the reduced basis as

$${\overline{x}}^i\approx \sum _{j=1}^{n_r}{b}_{ij}{v}^j,i=1,\cdots ,n_t,$$

and each $n_m$-dimensional component ${\overline{x}}_{\ell }^i$, $\ell =1,\cdots ,n_s$ of ${\overline{x}}^i$ (relating to (1)) can be represented in the reduced basis by

$$x(t_i;{\Theta }_{\ell })={\overline{x}}_{\ell }^i\approx \sum _{j=1}^{n_r}{b}_{ij}{\widehat{v}}_{\ell }^j,$$

where $v^j$ is decomposed into $n_s$ components as in (2)

$$v^j=\left[\begin{array}{c}{v}_1^j\\ ⋮\\ v_{n_m}^j\\ v_{n_m+1}^j\\ ⋮\\ v_{2n_m}^j\\ ⋮\\ v_{(n_s-1)n_m+1}^j\\ ⋮\\ v_{n_m·n_s}^j\end{array}\right]=:\left[\begin{array}{c}{\widehat{v}}_1^j\\ {\widehat{v}}_2^j\\ ⋮\\ {\widehat{v}}_{n_s}^j\end{array}\right].$$

Or analogously (in a continuous version of POD) as

$$\begin{array}{c}\hfill x(t;{\Theta }_{\ell })\approx \sum _{j=1}^{n_r}{b}_j\left(t\right){\widehat{v}}_{\ell }^j,\end{array}$$

(6)

for coefficient functions $b_j\left(t\right)=\langle x(t;{\Theta }_{\ell }),{\widehat{v}}^j\rangle$.

2.3. Construction of Data-Fitting Function

Next, we want to construct a data-fitting function $\Psi \left(t\right)\in {\mathbb{R}}^{n_m}$ that fits the known experimental data for the component $k=1,\cdots ,n_d$. We assume $\Psi \left(t\right)$ to be of the form

$$\Psi \left(t\right):=\sum _{\ell =1}^{n_s}{c}_{\ell }{\overline{x}}_{\ell }^k\left(t\right)=\sum _{\ell =1}^{n_s}{c}_{\ell }x(t;{\Theta }_{\ell }).$$

With (6), we obtain

$$\Psi \left(t\right)=\sum _{\ell =1}^{n_s}{c}_{\ell }\sum _{j=1}^{n_r}{b}_j\left(t\right){\widehat{v}}_{\ell }^j=\sum _{\ell =1}^{n_s}\sum _{j=1}^{n_r}{c}_{\ell }\langle x(t;{\Theta }_{\ell }),{\widehat{v}}_{\ell }^j\rangle {\widehat{v}}_{\ell }^j.$$

Usually only a few specific states of the dynamical system can be measured in experiments. Here, we assume that the first $n_d$ components of the state vector can be measured and we want to determine the coefficients $c_{\ell }$ in such a way that the first $n_d$ components of $\Psi ={\left[{\Psi }_k\right]}_{k=1}^{n_d}$ are an acceptable fit of the experimental data points, i.e.,  ${\Psi }_k\left({\tau }_s\right)\approx y_k\left({\tau }_s\right)$ for $s=1,\cdots ,n_{\tau }$ and $k=1,\cdots ,n_d$.

Thus, we construct a linear system to be solved for $c_{\ell }$:

$$\begin{array}{cc}\hfill y_k\left({\tau }_s\right)& \stackrel{!}{=}{\Psi }_k\left({\tau }_s\right)={\left[\sum _{\ell =1}^{n_s}\sum _{j=1}^{n_r}{c}_{\ell }{b}_j\left({\tau }_s\right){\widehat{v}}_{\ell }^j\right]}_k={\left[\sum _{\ell =1}^{n_s}\sum _{j=1}^{n_r}{c}_{\ell }\langle x({\tau }_s;{\Theta }_{\ell }),{\widehat{v}}_{\ell }^j\rangle {\widehat{v}}_{\ell }^j\right]}_k\hfill \end{array}$$

for all $k=1,\cdots ,n_d$ and $s=1,\cdots ,n_{\tau }$.

We obtain $n_d·n_{\tau }$ equations for the $n_s$ unknowns $c_{\ell }$, $\ell =1,\cdots ,n_s$. For each $k=1,\cdots ,n_d$, we obtain a linear system

$$A^{k}c=y^k$$

where $A^k=\left[a_{s\ell }^k\right]$, and  $a_{s\ell }^k:={\sum }_{j=1}^{n_r}{b}_j\left({\tau }_s\right)V((\ell -1)n_m+k,j),$ for $\ell =1,\cdots ,n_s$, $s=1,\cdots ,n_{\tau }$ and $y^k=\left[\begin{array}{c}{y}_k\left({\tau }_1\right)\\ ⋮\\ y_k\left({\tau }_{n_{\tau }}\right)\end{array}\right]$. Summarizing all $A^k$ and $y^k$ into one large system yields $Ac=y$ with $A\in {\mathbb{R}}^{n_d{n}_{\tau }\times n_s}$ and $y\in {\mathbb{R}}^{n_d{n}_{\tau }}$.

2.4. Estimating the Parameters

Once we have a representation for $\Psi \left(t\right)$ we can substitute $\Psi$ and

$${\dot{\Psi }}_i\approx \frac{\Psi \left(t_{i+1}\right)-\Psi \left(t_i\right)}{h_s}$$

as approximation for ${\dot{\Psi }}_i$ into the ODE (1) at the simulation time points $t_i$ to obtain

$$\frac{\Psi \left(t_{i+1}\right)-\Psi \left(t_i\right)}{h_s}\approx f(t_i,\Psi \left(t_i\right),\theta ),i=1,\cdots ,n_t.$$

Here, $\Psi$ can be assumed to be smooth since it is composed as linear combination of solution of an ODE. Inserting all information that we have gives

$$\begin{array}{cc}\hfill \frac{\Psi \left(t_{i+1}\right)-\Psi \left(t_i\right)}{h_s}& =f(t_i,\Psi \left(t_i\right),\theta )\hfill \\ \hfill \frac{{\sum }_{\ell =1}^{n_s}{c}_{\ell }x(t_{i+1};{\Theta }_{\ell })-{\sum }_{\ell =1}^{n_s}{c}_{\ell }x(t_i;{\Theta }_{\ell })}{h_s}& =f(t_i,\sum _{\ell =1}^{n_s}{c}_{\ell }x(t_i;{\Theta }_{\ell }),\theta ).\hfill \end{array}$$

This system can be solved for the parameter $\theta$ by a nonlinear system solver (e.g., fsolve in Matlab). We can sum up all these steps in the following Algorithm 1.

Algorithm 1 Parameter estimation for a kinetic model.

1. Compute $n_s$ solutions of a kinetic model given by (1) for different parameter values ${\Theta }_1,\cdots ,{\Theta }_{n_s}$ and form the snapshot matrix X. 2. Compute the SVD of X. Truncate and remain only the most dominant singular values $n_r$ to obtain the POD basis. 3. Express every vector in the snapshot matrix X as linear combination of the POD basis ${\overline{x}}^i={\sum }_{j=1}^{n_r}{b}_{ij}{v}^j$ with $b_{ij}={\left(v^j\right)}^T·{\overline{x}}^i$ for $i=1,\cdots ,n_t$. 4. Construct a data-fitting function that fits the known experimental data ( $y\left({\tau }_s\right)$) as a linear combination of POD basis $\Psi \left(t\right):={\sum }_{\ell =1}^{n_s}{c}_{\ell }{\overline{x}}_{\ell }^i\left(t\right)$. 5. Solve the linear system $\Psi \left({\tau }_s\right)=y\left({\tau }_s\right),s=1,\cdots ,n_{\tau }$ to obtain the value of $c_l,l=1,\cdots ,n_s$. 6. Substitute the value of $\psi$ into the original model $\frac{\Psi \left(t_{i+1}\right)-\Psi \left(t_i\right)}{h_s}\approx f(t_i,\Psi \left(t_i\right),\theta ),i=1,\cdots ,n_t$ and solve the system for $\theta$ parameter.

3. Application of the Parameter Estimation Method to a Kinetic Model of CCR in E. coli

In this section, we apply the parameter estimation procedure introduced in the previous section to the kinetic model of carbon catabolite repression (CCR) in E. coli [28]. CCR is the main regulatory mechanism in E. coli for the control of carbohydrate uptake. The regulatory network is strongly hierarchical with a regulator Crp (cyclic AMP receptor protein) on the top. Crp serves as accelerator for gene expression in case of low C source availability and in this way coordinates different subsystems of the cell, responsible for the uptake of carbon sources, their breakdown for the production of energy and precursors, and the conversion of the latter to biomass. In the second level, substrate specific regulator proteins (main repressor proteins) are inactivated in case the specific substrate is available in the growth medium. The interplay between regulators and the transcription apparatus results for example in the observation that only one substrate is taken up although two substrates are provided in the beginning. Two distinct growth phases are observed, and a sequential uptake of the substrates take place (diauxic growth).

3.1. The Kinetic Model of the Carbon Catabolite Repression

The main structure of the reaction scheme is outlined in Figure 1. Two extracellular substrates ($S_1$ and $S_2$) are taken up in enzyme catalysed reactions by their respective specific transporters ($E_1$ and $E_2$). The substrates have representatives inside the cell ($X_1$ and $X_2$). For example, in the case of glucose, $X_1$ represents the phosphorylated form glucose 6-phospahte while in the case of lactose, the substrate is not modified during the transport step. In general, transport systems for carbohydrate are inducible, that is, metabolites in the cell activate gene expression of the respective enzymes. In the model, these mechanisms are represented by a positive feedback from the $X_i$ to the synthesis of the respective enzymes $E_i$ (dashed blue lines). A special case for the interaction between two uptake systems in case of glucose and lactose is inducer exclusion (red dashed line); high levels of $X_1$ inhibit the activity of the second transporter for lactose uptake. The inhibition is so strong that the level of $X_2$ is minimal and, this way, the transporter $E_2$ cannot be synthesized. Only after they run out of glucose, the level of $X_1$ decreases and the inhibition is degraded. Both uptake branches converge in the central pathways, represented by metabolite M. From M, the synthesis of the main biomass compartment $B^{\prime }$ is described with a single reaction.

The system of differential equations associated with the chemical reactions network is given by

$$\dot{x}\left(t\right)=N·r\left(x\left(t\right)\right),$$

where $x={[S_1,S_2,X_1,X_2,E_1,E_2,M,B^{\prime }]}^T\in {\mathbb{R}}^8$ is the vector of metabolite concentrations, $r={[r_{s1},r_{s2},r_{d1},r_{d2},r_{e1},r_{e2},r_b]}^T$ is the vector of reaction rates and N is the stoichiometric matrix given as follows:

$N=\begin{array}{ll}& \begin{array}{lllllll}{r}_{s1}& r_{s2}& r_{d1}& r_{d2}& r_{e1}& r_{e2}& r_b\end{array}\\ \begin{array}{c}{S}_1\\ S_2\\ X_1\\ X_2\\ E_1\\ E_2\\ M\\ B^{\prime }\end{array}& \left[\begin{array}{ccccccc}-1& 0& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0& 0\\ 0& 0& 0& -1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 1& 1& 0& 0& -1\\ 0& 0& 0& 0& 0& 0& 1\end{array}\right]\in {\mathbb{R}}^{8\times 7}.\end{array}$

The kinetic rate laws $r_{s1},r_{s2}$ are defined by Michaelis–Menten kinetics [34] and $r_{d1},r_{d2},r_{e1},r_{e2},r_b$ by mass action kinetics [35] as follows:

$$\begin{array}{cc}\hfill r_{s1}& =k_{s1}{E}_1{\displaystyle \frac{S_1}{K_1+S_1}},r_{s2}=k_{s2}{E}_2{\displaystyle \frac{S_2}{K_2+S_2}},\hfill \\ \hfill r_{e1}& =k_{e1}{f}_1,r_{e2}=k_{e2}{f}_2,\hfill \\ \hfill r_{d1}& =k_{x1}{X}_1,r_{d2}=k_{x2}{X}_2,r_b=k_{m}M,\hfill \end{array}$$

where $k_{x1}$, $k_{x2}$, and $k_m$ are the constant rates all in 1/h unit, $K_1,K_2$ are the Michaelis–Menten constant in g/L unit and $k_{e1}$, $k_{e2}$ are the constant rates in the mol/gDW.h. The turnover number of enzymes $k_{s1},k_{s2}$ is in the unit 1/h.

The system of ordinary differential equations of the carbon catabolite repression is given by

$$\begin{array}{cc}\hfill \dot{S_i}& =-r_{si}{w}_{i}B,i=1,2\hfill \\ \hfill \dot{X_i}& =r_{si}-r_{di},\hfill \\ \hfill \dot{E_i}& =r_{ei}-\mu E_i,\hfill \\ \hfill \dot{M}& =r_{d1}+r_{d2}-r_b,\hfill \\ \hfill \dot{B}& =\mu B,\hfill \\ \hfill \dot{B^{\prime }}& =r_b-\mu B^{\prime },\hfill \end{array}$$

where $w_i$ is the molecular weight of the substrates (g/mol). The main biomass compartment $B^{\prime }$ consists of macromolecular species like protein, DNA, RNA, lipid, etc., and variable B is the entire biomass, growing with the specific growth rate $\mu$. The growth rate $\mu =Y_1{r}_{s1}+Y_2{r}_{s2}$ has the unit 1/h and $Y_1,Y_2$ are yield coefficients with unit gDW/mol. The process of induction is described with rates $f_1,f_2$ that are expressed in the following form

$$\begin{array}{cc}\hfill f_1& ={\displaystyle \frac{X_1}{{\beta }_1+X_1}},f_2={\displaystyle \frac{X_2}{{\beta }_2+X_2}},\hfill \end{array}$$

where ${\beta }_1,{\beta }_2$ are in the unit mol/gDW and $k_1$, $k_2$ are constant rates in mol/gDW.h unit. The intracellular components are assumed to be at the steady state i.e., $\dot{X_1}=0,\dot{X_2}=0$, and $\dot{M}=0$. This case study considers inducer exclusion where metabolite $X1$ works as a regulatory metabolite and inhibits enzyme $E2$ which is responsible for substrate $S2$ uptake. In a scheme, this is represented by an extension of the uptake rate $r_{s2}$ to $r_{s2}·{\displaystyle \frac{\alpha }{\alpha +X_1}}$, where $\alpha$ is a constant. Then the kinetic model of the above network is given as follows

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{S_1}& =-{\displaystyle \frac{k_{s1}·E_1·S_1}{K_1+S_1}}·w_1·B,\hfill \\ \hfill \dot{S_2}& =-{\displaystyle \frac{k_{s2}·E_2·S_2}{K_2+S_2}}·w_2·B·f_{inhibit},\hfill \\ \hfill \dot{E_1}& =k_{e1}{\displaystyle \frac{X_1}{{\beta }_1+X_1}}-\mu E_1,\hfill \\ \hfill \dot{E_2}& =k_{e2}{\displaystyle \frac{X_2}{{\beta }_2+X_2}}-\mu E_2,\hfill \\ \hfill \dot{B}& =\mu B,\hfill \\ \hfill \dot{B^{\prime }}& =k_m·M-\mu B^{\prime },\hfill \end{array}\end{array}$$

(7)

where

$$\begin{array}{cc}\hfill X_1& ={\displaystyle \frac{k_{s1}·E_1·S_1}{k_{x1}(K_1+S_1)}},X_2={\displaystyle \frac{k_{s2}·E_2·S2}{k_{x2}(K_2+S_2)}}·f_{inhibit}\hfill \\ \hfill M& ={\displaystyle \frac{k_{x1}·X_1}{k_m}}+{\displaystyle \frac{k_{x2}·X_2}{k_m}},f_{inhibit}={\displaystyle \frac{\alpha }{\alpha +X_1}}.\hfill \end{array}$$

3.2. Numerical Results

In the following we present some examples in which we estimate certain parameter values of the kinetic model (7). All results presented in the following sections were computed with MATLAB (R2022a).

3.2.1. Estimation of the Parameter $k_{s1}$

In this first example we estimate one parameter value, namely, the uptake rate constant $k_{s1}$. First, we compute the snapshot matrix X from the simulation of the kinetic model (7) for different values of the parameter $k_{s1}$ over the time interval $[0,7]$ h. The different values of parameter $k_{s1}$ are equidistantly distributed in the interval $[{10}^{-5},{10}^{-2}]$, which are the lower and upper bound of parameter value, respectively. In each of $n_s$ simulation, we evaluate at $n_t$ simulation time points. In this example, we set $n_s=50$. The remaining constant reaction rates are given in

Table 1. The set of all parameters of the kinetic model that are taken from literature [28].

Constant Rates	Value	Unit
$K_1$	$1\times {10}^{-4}$	g/L
$K_2$	$0.2$	g/L
$k_{e1}$	6	mol/gDWh
$k_{e2}$	$8.69$	mol/gDWh
$k_{x1}$	10	1/h
$k_{x2}$	10	1/h
$k_m$	10	1/h
${\beta }_1$	1	mol/gDW
${\beta }_2$	$0.6$	mol/gDW
$\alpha$	$11.05$	mol/gDW
$Y_1$	90	gDW/mol
$Y_2$	$102.6$	gDW/mol
$w_1$	180	gDW/mol
$w_2$	342	gDW/mol

and the parameter value of $k_{s2}=0.0033$ is assumed. We use the MATLAB function ode15s with tolerances $RTOL=ATOL={10}^{-6}$ to compute a numerical solution of the ODE system (7) with 50 different values for the parameter $k_{s1}$. The initial concentrations of extracellular and intracellular metabolites are assumed to be $(S_1,S_2,E_1,E_2,B,B^{\prime })=(0.22,1.18,0.1,0.1,0.032,0.5)\in {\mathbb{R}}^6$.

We used the MATLAB function svd to calculate the singular value decomposition of X. The singular values are shown in Figure 2.

The behavior of the model order reduction method strongly depends on the decay of the singular values of the snapshot matrix. Figure 2 shows gradually decaying singular values with a strong decay for the smallest singular values indicating that ignoring these will not lead to any considerable loss of information. We truncate the most dominant singular values such that the relative information coefficient I is close to one. In this example, we use $n_r=60$ where $I=0.99$. The date-fitting function $\Psi \left(t\right)$ is constructed as in Section 2.3 to fit the experimental data set for the three measurement components $S_1,S_2$, and B at the measurement time points $n_{\tau }=15$. The trajectories of the measurement data set and the data-fitting function are depicted in Figure 3.

We can observe that the data-fitting function fits well with the measurement data. Once the data-fitting function $\Psi \left(t\right)$ is obtained, we can substitute the function $\Psi \left(t\right)$ in the kinetic model (7) and solve the system for the parameter $k_{s1}$. We use the MATLAB function fsolve with function tolerance $TOL={10}^{-6}$ using the Levenberg–Marquardt algorithm with initial guess $k_{s1}=0.0167$. The solution gives an estimated parameter value of $k_{best}=0.018$. For the purpose of comparison, we also provide results from the least square method (LSQ). We perform the LSQ method using MATLAB function lsqnonlin and lhsdesign with the same assumptions (e.g., number of fits $n_s$ and lower and upper bounds) as in the POD method. The solution gives an estimated value $k_{lsq}=0.01$ of the parameter $k_{s1}$. The simulation results of the kinetic model with $k_{best}$ in comparison with the simulation results of the kinetic model with $k_{lsq}$, and the measurement data are depicted in Figure 4.

We can observe that the dynamical behavior of the kinetic model states with parameter value $k_{best}$ matches well with the experimental data for the substrates $S_1$, $S_2$, and the Biomass B. While the simulation of the solution of the kinetic model with parameter value $k_{lsq}$ for the substrate $S_1$ is unable to match well with the data. The substrate $S_1$ spends more time to consume compared to the measured data. In contrast, the trajectories of $S_2$ and Biomass still fit good with the measured data. From the above figure, the POD method gives better results than the LSQ method. Note, that with a change of the lower and upper bounds in the LSQ method, we could obtain a better result. Both methods are sensitive to assumptions, which means that the lower and upper bound of parameter value $lb,ub$, the mesh size M, and the reduced basis $n_r$ can be adjusted to obtain a good fitting.

3.2.2. Estimation of the Parameters $k_{s1}$ and $k_{s2}$

In this example, we consider the two parameters for the uptake rate constants $k_{s1}$ and $k_{s2}$. The snapshot matrices X is computed from $n_s=40$ simulation of the kinetic model equations (7) over the time interval. We use a Latin hypercube sampling (LHS) to obtain different values of the parameters via the MATLAB function lhsdesign. To generate parameter sampling for the parameters $k_{s1}$ and $k_{s2}$, we set the lower and upper bound for the parameters as $lb={10}^{-5}$, $ub={10}^{-2}$. We use the MATLAB function svd to calculate the singular value decomposition of X. The singular values are depicted in Figure 5.

Here, we truncate at $n_r=60$, with relative information coefficient $I=0.9$. Substituting $\psi$ in the kinetic model and solving the system for the two parameters using MATLAB function fsolve with function tolerance ${10}^{-6}$ and initial guesses $k_{s1}=0.0167,k_{s2}=0.0033$, we obtain the estimated values of the parameters $k_{best}=[0.0194,0.0032]$. In this example, we perform the LSQ method with different assumptions of bounds. We set the lower and upper bounds in the LSQ method to be $lb={10}^{-3}$, $ub={10}^{-1}$, we obtain the estimated values $k_{lsq}=[0.0163,0.0029]$.

From Figure 6, we can observe that the trajectories for the experimental data and the simulation results of kinetic model with the estimated parameter values using POD method match good for the substrates trajectories $S_1,S_2$ and Biomass curve. While the simulation of the kinetic model with parameters value using the LSQ method is slightly different from the measured data.

3.2.3. Estimation of the Parameters $k_{s1},k_{e1},k_{x1}$ and $k_{x2}$

In this example, we apply the method to estimate four parameters in the kinetic model (7). We use the MATLAB function lhsdesign to generate parameter sampling for the parameters $k_{s1},k_{s2},k_{e2}$, and ${\beta }_2$, where the lower and upper bound for the parameters are $[{10}^{-2},{10}^{+1}]$. We set the number of simulation to $n_s=20$. The snapshot matrix X is computed and the SVD is applied; Figure 7 shows the singular values. We choose $n_r=60$ and the system is solved for the four parameters using MATLAB function lsqnonlin with options algorithm levenberg-marquardt and initial guess of parameters $k_{s1}=0.0167$, $k_{e1}=6$, $k_{x1}=10$, and $k_{x2}=10$. The estimated values of the parameters are $kbest=[0.01,125.7,103.6,9.97]$. We perform the LSQ method using the same assumptions as in the POD method, the solution gives an estimated values $k_{lsq}=[0.059,0.084,0.333,9.99]$.

At the beginning of the simulation, we can observe from Figure 8 that the trajectory of $S_1$ using the estimated values from the POD method consumes slowly compared to the trajectories by the LSQ method and the experimental data. Around the time 2 h, it is exhausted a bit earlier than other curves. The trajectories of $S_2$ and B from the POD method are qualitatively similar to the measured data while their trajectories of the LSQ method are slightly different.

4. Conclusions and Discussion

In this paper, we have discussed the parameter estimation procedure using model order reduction by the POD method for a kinetic model of a cellular system. The process of finding parameter values is considered to be a main challenge in mathematical modelling of biological system, so the possibility to use the POD method in the parameter estimation problem is of crucial importance. We have applied the POD method and the singular value decomposition to obtain the reduced basis elements, which are used for constructing a function that fits the trajectory of the kinetic model to the observable data set of the biological system. The parameter estimation using the POD method is capable to estimate the uptake rate constants and other parameters of the CCR kinetic model. The solution of the kinetic model with the estimated parameter values matches good with the experimental data set. We performed the least square method to estimate parameter values of the CRR kinetic model for the sake of comparison with the POD method and show to what extent the POD method could estimate the values of parameters of the CRR kinetic model. We could observe that the number of reduced basis, the lower and upper bounds, and the number of simulations play an important role in the parameter estimation problem. For every parameter estimation case, these numbers can be adjusted until the desired results are obtained. In addition, the presented algorithm of parameter estimation can be used as a reference for different kinetic models of biological systems. The presented approach has some limitations as it strongly depends on the behavior of the singular value spectrum. If there is only a small number of dominant eigenmodes, then the singular values will rapidly decay and only a small number of them is enough to capture the characteristic behavior of the system. However, in some ill-posed cases, this might not be the case and the heuristic approach presented in the algorithm might not lead to satisfactory results. The choice of a suitable threshold for discarding certain singular values also depends on the quality of the data, in particular, rounding and approximation errors. It has been shown in [36] that the truncated SVD is a suitable method for the regularization of an ill-posed problem when the coefficient matrix is ill-conditioned with a well-determined numerical rank. In this case, the solution obtained by the truncated SVD with a truncation threshold equal to the numerical rank of the matrix is guaranteed to be similar to the regularized solution where the regularization parameter is chosen near its intuitive optimum value.